1# Copyright (C) 2011 Atsushi Togo
2# All rights reserved.
3#
4# This file is part of phonopy.
5#
6# Redistribution and use in source and binary forms, with or without
7# modification, are permitted provided that the following conditions
8# are met:
9#
10# * Redistributions of source code must retain the above copyright
11#   notice, this list of conditions and the following disclaimer.
12#
13# * Redistributions in binary form must reproduce the above copyright
14#   notice, this list of conditions and the following disclaimer in
15#   the documentation and/or other materials provided with the
16#   distribution.
17#
18# * Neither the name of the phonopy project nor the names of its
19#   contributors may be used to endorse or promote products derived
20#   from this software without specific prior written permission.
21#
22# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
23# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
24# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
25# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
26# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
27# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
28# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
29# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
30# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
32# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
33# POSSIBILITY OF SUCH DAMAGE.
34
35import numpy as np
36from phonopy.units import VaspToTHz
37from phonopy.structure.grid_points import GridPoints
38
39
40class MeshBase(object):
41    """Base class of Mesh and IterMesh classes
42
43    Attributes
44    ----------
45    mesh_numbers: ndarray
46        Mesh numbers along a, b, c axes.
47        dtype='intc'
48        shape=(3,)
49    qpoints: ndarray
50        q-points in reduced coordinates of reciprocal lattice
51        dtype='double'
52        shape=(ir-grid points, 3)
53    weights: ndarray
54        Geometric q-point weights. Its sum is the number of grid points.
55        dtype='intc'
56        shape=(ir-grid points,)
57    grid_address: ndarray
58        Addresses of all grid points represented by integers.
59        dtype='intc'
60        shape=(prod(mesh_numbers), 3)
61    ir_grid_points: ndarray
62        Indices of irreducibple grid points in grid_address.
63        dtype='intc'
64        shape=(ir-grid points,)
65    grid_mapping_table: ndarray
66        Index mapping table from all grid points to ir-grid points.
67        dtype='intc'
68        shape=(prod(mesh_numbers),)
69    dynamical_matrix: DynamicalMatrix
70        Dynamical matrix instance to compute dynamical matrix at q-points.
71
72    """
73
74    def __init__(self,
75                 dynamical_matrix,
76                 mesh,
77                 shift=None,
78                 is_time_reversal=True,
79                 is_mesh_symmetry=True,
80                 with_eigenvectors=False,
81                 is_gamma_center=False,
82                 rotations=None,  # Point group operations in real space
83                 factor=VaspToTHz):
84        self._mesh = np.array(mesh, dtype='intc')
85        self._with_eigenvectors = with_eigenvectors
86        self._factor = factor
87        self._cell = dynamical_matrix.primitive
88        self._dynamical_matrix = dynamical_matrix
89
90        self._gp = GridPoints(self._mesh,
91                              np.linalg.inv(self._cell.cell),
92                              q_mesh_shift=shift,
93                              is_gamma_center=is_gamma_center,
94                              is_time_reversal=(is_time_reversal and
95                                                is_mesh_symmetry),
96                              rotations=rotations,
97                              is_mesh_symmetry=is_mesh_symmetry)
98
99        self._qpoints = self._gp.qpoints
100        self._weights = self._gp.weights
101
102        self._frequencies = None
103        self._eigenvectors = None
104
105        self._q_count = 0
106
107    @property
108    def mesh_numbers(self):
109        return self._mesh
110
111    def get_mesh_numbers(self):
112        return self.mesh_numbers
113
114    @property
115    def qpoints(self):
116        return self._qpoints
117
118    def get_qpoints(self):
119        return self.qpoints
120
121    @property
122    def weights(self):
123        return self._weights
124
125    def get_weights(self):
126        return self.weights
127
128    @property
129    def grid_address(self):
130        return self._gp.grid_address
131
132    def get_grid_address(self):
133        return self.grid_address
134
135    @property
136    def ir_grid_points(self):
137        return self._gp.ir_grid_points
138
139    def get_ir_grid_points(self):
140        return self.ir_grid_points
141
142    @property
143    def grid_mapping_table(self):
144        return self._gp.grid_mapping_table
145
146    def get_grid_mapping_table(self):
147        return self.grid_mapping_table
148
149    @property
150    def dynamical_matrix(self):
151        return self._dynamical_matrix
152
153    def get_dynamical_matrix(self):
154        return self.dynamical_matrix
155
156    @property
157    def with_eigenvectors(self):
158        return self._with_eigenvectors
159
160
161class Mesh(MeshBase):
162    """Class for phonons on mesh grid
163
164    Frequencies and eigenvectors can be also accessible by iterator
165    representation to be compatible with IterMesh.
166
167    Attributes
168    ----------
169    frequencies: ndarray
170        Phonon frequencies at ir-grid points. Imaginary frequenies are
171        represented by negative real numbers.
172        dtype='double'
173        shape=(ir-grid points, bands)
174    eigenvectors: ndarray
175        Phonon eigenvectors at ir-grid points. See the data structure at
176        np.linalg.eigh.
177        shape=(ir-grid points, bands, bands)
178        dtype=complex of "c%d" % (np.dtype('double').itemsize * 2)
179        order='C'
180    group_velocities: ndarray
181        Phonon group velocities at ir-grid points.
182        shape=(ir-grid points, bands, 3)
183        dtype='double'
184    More attributes from MeshBase should be watched.
185
186    """
187    def __init__(self,
188                 dynamical_matrix,
189                 mesh,
190                 shift=None,
191                 is_time_reversal=True,
192                 is_mesh_symmetry=True,
193                 with_eigenvectors=False,
194                 is_gamma_center=False,
195                 group_velocity=None,
196                 rotations=None,  # Point group operations in real space
197                 factor=VaspToTHz,
198                 use_lapack_solver=False):
199        MeshBase.__init__(self,
200                          dynamical_matrix,
201                          mesh,
202                          shift=shift,
203                          is_time_reversal=is_time_reversal,
204                          is_mesh_symmetry=is_mesh_symmetry,
205                          with_eigenvectors=with_eigenvectors,
206                          is_gamma_center=is_gamma_center,
207                          rotations=rotations,
208                          factor=factor)
209
210        self._group_velocity = group_velocity
211        self._group_velocities = None
212        self._use_lapack_solver = use_lapack_solver
213
214    def __iter__(self):
215        if self._frequencies is None:
216            self.run()
217        return self
218
219    def next(self):
220        return self.__next__()
221
222    def __next__(self):
223        if self._q_count == len(self._qpoints):
224            self._q_count = 0
225            raise StopIteration
226        else:
227            i = self._q_count
228            self._q_count += 1
229            if self._eigenvectors is None:
230                return self._frequencies[i], None
231            else:
232                return self._frequencies[i], self._eigenvectors[i]
233
234    def run(self):
235        self._set_phonon()
236        if self._group_velocity is not None:
237            self._set_group_velocities(self._group_velocity)
238
239    @property
240    def frequencies(self):
241        if self._frequencies is None:
242            self.run()
243        return self._frequencies
244
245    def get_frequencies(self):
246        return self.frequencies
247
248    @property
249    def eigenvectors(self):
250        """
251        Eigenvectors is a numpy array of three dimension.
252        The first index runs through q-points.
253        In the second and third indices, eigenvectors obtained
254        using numpy.linalg.eigh are stored.
255
256        The third index corresponds to the eigenvalue's index.
257        The second index is for atoms [x1, y1, z1, x2, y2, z2, ...].
258        """
259        if self._frequencies is None:
260            self.run()
261        return self._eigenvectors
262
263    def get_eigenvectors(self):
264        return self.eigenvectors
265
266    @property
267    def group_velocities(self):
268        if self._frequencies is None:
269            self.run()
270        return self._group_velocities
271
272    def get_group_velocities(self):
273        return self.group_velocities
274
275    def write_hdf5(self):
276        import h5py
277        with h5py.File('mesh.hdf5', 'w') as w:
278            w.create_dataset('mesh', data=self._mesh)
279            w.create_dataset('qpoint', data=self._qpoints)
280            w.create_dataset('weight', data=self._weights)
281            w.create_dataset('frequency', data=self._frequencies)
282            if self._eigenvectors is not None:
283                w.create_dataset('eigenvector', data=self._eigenvectors)
284            if self._group_velocities is not None:
285                w.create_dataset('group_velocity', data=self._group_velocities)
286
287    def write_yaml(self):
288        natom = self._cell.get_number_of_atoms()
289        rec_lattice = np.linalg.inv(self._cell.get_cell())  # column vectors
290        distances = np.sqrt(
291            np.sum(np.dot(self._qpoints, rec_lattice.T) ** 2, axis=1))
292
293        lines = []
294
295        lines.append("mesh: [ %5d, %5d, %5d ]" % tuple(self._mesh))
296        lines.append("nqpoint: %-7d" % self._qpoints.shape[0])
297        lines.append("reciprocal_lattice:")
298        for vec, axis in zip(rec_lattice.T, ('a*', 'b*', 'c*')):
299            lines.append("- [ %12.8f, %12.8f, %12.8f ] # %2s" %
300                         (tuple(vec) + (axis,)))
301        lines.append("natom:   %-7d" % natom)
302        lines.append(str(self._cell))
303        lines.append("")
304        lines.append("phonon:")
305
306        for i, (q, d) in enumerate(zip(self._qpoints, distances)):
307            lines.append("- q-position: [ %12.7f, %12.7f, %12.7f ]"
308                         % tuple(q))
309            lines.append("  distance_from_gamma: %12.9f" % d)
310            lines.append("  weight: %-5d" % self._weights[i])
311            lines.append("  band:")
312
313            for j, freq in enumerate(self._frequencies[i]):
314                lines.append("  - # %d" % (j + 1))
315                lines.append("    frequency:  %15.10f" % freq)
316
317                if self._group_velocities is not None:
318                    lines.append("    group_velocity: "
319                                 "[ %13.7f, %13.7f, %13.7f ]" %
320                                 tuple(self._group_velocities[i, j]))
321
322                if self._with_eigenvectors:
323                    lines.append("    eigenvector:")
324                    for k in range(natom):
325                        lines.append("    - # atom %d" % (k+1))
326                        for l in (0, 1, 2):
327                            lines.append(
328                                "      - [ %17.14f, %17.14f ]"
329                                % (self._eigenvectors[i, k*3+l, j].real,
330                                   self._eigenvectors[i, k*3+l, j].imag))
331            lines.append("")
332
333        with open('mesh.yaml', 'w') as w:
334            w.write("\n".join(lines))
335
336    def _set_phonon(self):
337        num_band = self._cell.get_number_of_atoms() * 3
338        num_qpoints = len(self._qpoints)
339
340        self._frequencies = np.zeros((num_qpoints, num_band), dtype='double')
341        if self._with_eigenvectors or self._use_lapack_solver:
342            dtype = "c%d" % (np.dtype('double').itemsize * 2)
343            self._eigenvectors = np.zeros(
344                (num_qpoints, num_band, num_band,), dtype=dtype, order='C')
345
346        if self._use_lapack_solver:
347            from phono3py.phonon.solver import get_phonons_at_qpoints
348            get_phonons_at_qpoints(self._frequencies,
349                                   self._eigenvectors,
350                                   self._dynamical_matrix,
351                                   self._qpoints,
352                                   self._factor,
353                                   nac_q_direction=None,
354                                   lapack_zheev_uplo='L')
355        else:
356            for i, q in enumerate(self._qpoints):
357                self._dynamical_matrix.run(q)
358                dm = self._dynamical_matrix.dynamical_matrix
359                if self._with_eigenvectors:
360                    eigvals, self._eigenvectors[i] = np.linalg.eigh(dm)
361                    eigenvalues = eigvals.real
362                else:
363                    eigenvalues = np.linalg.eigvalsh(dm).real
364                self._frequencies[i] = np.array(np.sqrt(abs(eigenvalues)) *
365                                                np.sign(eigenvalues),
366                                                dtype='double',
367                                                order='C') * self._factor
368
369    def _set_group_velocities(self, group_velocity):
370        group_velocity.run(self._qpoints)
371        self._group_velocities = group_velocity.group_velocities
372
373
374class IterMesh(MeshBase):
375    """Generator class for phonons on mesh grid
376
377    Not like as Mesh class, frequencies and eigenvectors are not
378    stored, instead generated by iterator. This may be used for
379    saving memory space even with very dense samplig mesh.
380
381    Attributes
382    ----------
383    Attributes from MeshBase should be watched.
384
385    """
386    def __init__(self,
387                 dynamical_matrix,
388                 mesh,
389                 shift=None,
390                 is_time_reversal=True,
391                 is_mesh_symmetry=True,
392                 with_eigenvectors=False,
393                 is_gamma_center=False,
394                 rotations=None,  # Point group operations in real space
395                 factor=VaspToTHz):
396        MeshBase.__init__(self,
397                          dynamical_matrix,
398                          mesh,
399                          shift=shift,
400                          is_time_reversal=is_time_reversal,
401                          is_mesh_symmetry=is_mesh_symmetry,
402                          with_eigenvectors=with_eigenvectors,
403                          is_gamma_center=is_gamma_center,
404                          rotations=rotations,
405                          factor=factor)
406
407    def __iter__(self):
408        return self
409
410    def next(self):
411        return self.__next__()
412
413    def __next__(self):
414        if self._q_count == len(self._qpoints):
415            self._q_count = 0
416            raise StopIteration
417        else:
418            q = self._qpoints[self._q_count]
419            self._dynamical_matrix.run(q)
420            dm = self._dynamical_matrix.dynamical_matrix
421            if self._with_eigenvectors:
422                eigvals, eigenvectors = np.linalg.eigh(dm)
423                eigenvalues = eigvals.real
424            else:
425                eigenvalues = np.linalg.eigvalsh(dm).real
426            frequencies = np.array(np.sqrt(abs(eigenvalues)) *
427                                   np.sign(eigenvalues),
428                                   dtype='double',
429                                   order='C') * self._factor
430            self._q_count += 1
431            return frequencies, eigenvectors
432