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