1"""Phonopy class.""" 2 3# Copyright (C) 2015 Atsushi Togo 4# All rights reserved. 5# 6# This file is part of phonopy. 7# 8# Redistribution and use in source and binary forms, with or without 9# modification, are permitted provided that the following conditions 10# are met: 11# 12# * Redistributions of source code must retain the above copyright 13# notice, this list of conditions and the following disclaimer. 14# 15# * Redistributions in binary form must reproduce the above copyright 16# notice, this list of conditions and the following disclaimer in 17# the documentation and/or other materials provided with the 18# distribution. 19# 20# * Neither the name of the phonopy project nor the names of its 21# contributors may be used to endorse or promote products derived 22# from this software without specific prior written permission. 23# 24# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 25# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 26# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS 27# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE 28# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, 29# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, 30# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 31# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 32# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 33# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN 34# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 35# POSSIBILITY OF SUCH DAMAGE. 36 37import sys 38import warnings 39import textwrap 40import numpy as np 41from phonopy.version import __version__ 42from phonopy.interface.phonopy_yaml import PhonopyYaml 43from phonopy.structure.atoms import PhonopyAtoms 44from phonopy.structure.symmetry import Symmetry, symmetrize_borns_and_epsilon 45from phonopy.structure.grid_points import length2mesh 46from phonopy.structure.cells import ( 47 get_supercell, get_primitive, guess_primitive_matrix, 48 get_primitive_matrix, shape_supercell_matrix) 49from phonopy.structure.dataset import ( 50 get_displacements_and_forces, forces_in_dataset) 51from phonopy.harmonic.displacement import ( 52 get_least_displacements, directions_to_displacement_dataset, 53 get_random_displacements_dataset) 54from phonopy.harmonic.force_constants import ( 55 symmetrize_force_constants, symmetrize_compact_force_constants, 56 show_drift_force_constants, cutoff_force_constants, 57 set_tensor_symmetry_PJ) 58from phonopy.harmonic.force_constants import get_fc2 as get_phonopy_fc2 59from phonopy.interface.calculator import get_default_physical_units 60from phonopy.interface.fc_calculator import get_fc2 61from phonopy.harmonic.dynamical_matrix import get_dynamical_matrix 62from phonopy.phonon.band_structure import ( 63 BandStructure, get_band_qpoints_by_seekpath) 64from phonopy.phonon.thermal_properties import ThermalProperties 65from phonopy.phonon.mesh import Mesh, IterMesh 66from phonopy.units import VaspToTHz 67from phonopy.phonon.dos import TotalDos, PartialDos 68from phonopy.phonon.thermal_displacement import ( 69 ThermalDisplacements, ThermalDisplacementMatrices) 70from phonopy.phonon.random_displacements import RandomDisplacements 71from phonopy.phonon.animation import write_animation 72from phonopy.phonon.modulation import Modulation 73from phonopy.phonon.qpoints import QpointsPhonon 74from phonopy.phonon.irreps import IrReps 75from phonopy.phonon.group_velocity import GroupVelocity 76from phonopy.phonon.moment import PhononMoment 77from phonopy.spectrum.dynamic_structure_factor import DynamicStructureFactor 78 79 80class Phonopy(object): 81 """Phonopy main API given as a class. 82 83 Attributes 84 ---------- 85 version : str 86 unitcell : PhonopyAtoms 87 primitive : Primitive 88 supercell : Supercell 89 symmetry : Symmetry 90 Symmetry of supercell. 91 primitive_symmetry : Symmetry 92 Symmetry of primitive cell. 93 supercell_matrix : ndarray 94 shape=(3, 3), dtype='intc', order='C'. 95 primitive_matrix : ndarray 96 shape=(3, 3), dtype='double', order='C'. 97 unit_conversion_factor : float 98 Phonon frequency unit conversion factor. 99 calculator : str 100 dataset : dict 101 displacements : ndarray or list of list (getter) and array-like (setter). 102 forces : ndarray (getter) or array_like (setter). 103 force_constants : ndarray (getter) and array_like (setter). 104 nac_params : dict 105 supercells_with_displacements : list of PhonopyAtoms. 106 dynamical_matrix : DynamicalMatrix 107 108 qpoints : QpointsPhonon 109 band_structure : BandStructure 110 mesh : Mesh or IterMesh 111 thermal_properties : ThermalProperties 112 thermal_displacements : ThermalDisplacements 113 thermal_displacement_matrix : ThermalDisplacementMatrices 114 random_displacements : RandomDisplacements 115 dynamic_structure_factor : DynamicStructureFactor. 116 irreps : IrReps 117 moment : PhononMoment 118 total_dos : TotalDos 119 120 """ 121 122 def __init__(self, 123 unitcell, 124 supercell_matrix=None, 125 primitive_matrix=None, 126 nac_params=None, 127 factor=VaspToTHz, 128 frequency_scale_factor=None, 129 dynamical_matrix_decimals=None, 130 force_constants_decimals=None, 131 group_velocity_delta_q=None, 132 symprec=1e-5, 133 is_symmetry=True, 134 store_dense_svecs=False, 135 calculator=None, 136 use_lapack_solver=False, 137 log_level=0): 138 """Init Phonopy API.""" 139 self._symprec = symprec 140 self._factor = factor 141 self._frequency_scale_factor = frequency_scale_factor 142 self._is_symmetry = is_symmetry 143 self._calculator = calculator 144 self._store_dense_svecs = store_dense_svecs 145 self._use_lapack_solver = use_lapack_solver 146 self._log_level = log_level 147 148 # Create supercell and primitive cell 149 self._unitcell = PhonopyAtoms(atoms=unitcell) 150 self._supercell_matrix = self._shape_supercell_matrix(supercell_matrix) 151 if isinstance(primitive_matrix, str): 152 self._primitive_matrix = self._set_primitive_matrix( 153 primitive_matrix) 154 elif primitive_matrix is not None: 155 self._primitive_matrix = np.array(primitive_matrix, 156 dtype='double', order='c') 157 else: 158 self._primitive_matrix = None 159 self._supercell = None 160 self._primitive = None 161 self._build_supercell() 162 self._build_primitive_cell() 163 164 # Set supercell and primitive symmetry 165 self._symmetry = None 166 self._primitive_symmetry = None 167 self._search_symmetry() 168 self._search_primitive_symmetry() 169 170 # displacements 171 self._displacement_dataset = {'natom': 172 self._supercell.get_number_of_atoms()} 173 self._supercells_with_displacements = None 174 175 # set_force_constants or set_forces 176 self._force_constants = None 177 self._force_constants_decimals = force_constants_decimals 178 179 # set_dynamical_matrix 180 self._dynamical_matrix = None 181 self._nac_params = nac_params 182 self._dynamical_matrix_decimals = dynamical_matrix_decimals 183 184 # set_band_structure 185 self._band_structure = None 186 187 # set_mesh 188 self._mesh = None 189 190 # set_tetrahedron_method 191 self._tetrahedron_method = None 192 193 # set_thermal_properties 194 self._thermal_properties = None 195 196 # set_thermal_displacements 197 self._thermal_displacements = None 198 199 # set_thermal_displacement_matrices 200 self._thermal_displacement_matrices = None 201 202 # set_dynamic_structure_factor 203 self._dynamic_structure_factor = None 204 205 # set_partial_DOS 206 self._pdos = None 207 208 # set_total_DOS 209 self._total_dos = None 210 211 # set_modulation 212 self._modulation = None 213 214 # set_character_table 215 self._irreps = None 216 217 # set_group_velocity 218 self._group_velocity = None 219 self._gv_delta_q = group_velocity_delta_q 220 221 @property 222 def version(self): 223 """Return phonopy release version number. 224 225 str 226 Phonopy release version number 227 228 """ 229 return __version__ 230 231 def get_version(self): 232 """Return phonopy release version number.""" 233 warnings.warn("Phonopy.get_version() is deprecated." 234 "Use Phonopy.version attribute.", 235 DeprecationWarning) 236 return self.version 237 238 @property 239 def primitive(self): 240 """Return primitive cell. 241 242 Primitive 243 Primitive cell. 244 245 """ 246 return self._primitive 247 248 def get_primitive(self): 249 """Return primitive cell.""" 250 warnings.warn("Phonopy.get_primitive() is deprecated." 251 "Use Phonopy.primitive attribute.", 252 DeprecationWarning) 253 return self.primitive 254 255 @property 256 def unitcell(self): 257 """Return input unit cell. 258 259 PhonopyAtoms 260 Input unit cell. 261 262 """ 263 return self._unitcell 264 265 def get_unitcell(self): 266 """Return input unit cell.""" 267 warnings.warn("Phonopy.get_unitcell() is deprecated." 268 "Use Phonopy.unitcell attribute.", 269 DeprecationWarning) 270 return self.unitcell 271 272 def set_unitcell(self, unitcell): 273 """Set input unit cell.""" 274 warnings.warn("Phonopy.set_unitcell() is deprecated." 275 "No way to set unit cell will be provided.", 276 DeprecationWarning) 277 self._unitcell = unitcell 278 self._build_supercell() 279 self._build_primitive_cell() 280 self._search_symmetry() 281 self._search_primitive_symmetry() 282 self._displacement_dataset = None 283 284 @property 285 def supercell(self): 286 """Return supercell. 287 288 Supercell 289 Supercell. 290 291 """ 292 return self._supercell 293 294 def get_supercell(self): 295 """Return supercell.""" 296 warnings.warn("Phonopy.get_supercell() is deprecated." 297 "Use Phonopy.supercell attribute.", 298 DeprecationWarning) 299 return self.supercell 300 301 @property 302 def symmetry(self): 303 """Return symmetry of supercell. 304 305 Symmetry 306 Symmetry of supercell. 307 308 """ 309 return self._symmetry 310 311 def get_symmetry(self): 312 """Return symmetry of supercell.""" 313 warnings.warn("Phonopy.get_symmetry() is deprecated." 314 "Use Phonopy.symmetry attribute.", 315 DeprecationWarning) 316 return self.symmetry 317 318 @property 319 def primitive_symmetry(self): 320 """Return symmetry of primitive cell. 321 322 Symmetry 323 Symmetry of primitive cell. 324 325 """ 326 return self._primitive_symmetry 327 328 def get_primitive_symmetry(self): 329 """Return symmetry of primitive cell.""" 330 warnings.warn("Phonopy.get_primitive_symmetry() is deprecated." 331 "Use Phonopy.primitive_symmetry attribute.", 332 DeprecationWarning) 333 return self.primitive_symmetry 334 335 @property 336 def supercell_matrix(self): 337 """Return transformation matrix to supercell cell from unit cell. 338 339 ndarray 340 Supercell matrix with respect to unit cell. 341 shape=(3, 3), dtype='intc', order='C'. 342 343 """ 344 return self._supercell_matrix 345 346 def get_supercell_matrix(self): 347 """Return transformation matrix to supercell cell from unit cell.""" 348 warnings.warn("Phonopy.get_supercell_matrix() is deprecated." 349 "Use Phonopy.supercell_matrix attribute.", 350 DeprecationWarning) 351 return self.supercell_matrix 352 353 @property 354 def primitive_matrix(self): 355 """Return transformation matrix to primitive cell from unit cell. 356 357 ndarray 358 Primitive matrix with respect to unit cell. 359 shape=(3, 3), dtype='double', order='C'. 360 361 """ 362 return self._primitive_matrix 363 364 def get_primitive_matrix(self): 365 """Return transformation matrix to primitive cell from unit cell.""" 366 warnings.warn("Phonopy.get_primitive_matrix() is deprecated." 367 "Use Phonopy.primitive_matrix attribute.", 368 DeprecationWarning) 369 return self.primitive_matrix 370 371 @property 372 def unit_conversion_factor(self): 373 """Return phonon frequency unit conversion factor. 374 375 float 376 Phonon frequency unit conversion factor. This factor 377 converts sqrt(<force>/<distance>/<AMU>)/2pi/1e12 to the 378 other favorite phonon frequency unit. Normally this factor 379 is recommended to be that converts to THz (ordinary 380 frequency) to calculate a variety of phonon properties 381 that assumes that input phonon frequencies have THz unit. 382 383 """ 384 return self._factor 385 386 def get_unit_conversion_factor(self): 387 """Return phonon frequency unit conversion factor.""" 388 warnings.warn("Phonopy.get_unit_conversion_factor() is deprecated." 389 "Use Phonopy.unit_conversion_factor attribute.", 390 DeprecationWarning) 391 return self.unit_conversion_factor 392 393 @property 394 def calculator(self): 395 """Return calculator name. 396 397 str 398 Calculator name such as 'vasp', 'qe', etc. 399 400 """ 401 return self._calculator 402 403 @property 404 def dataset(self): 405 """Return dataset to store displacements and forces. 406 407 Dataset containing information of displacements in supercells. 408 This optionally contains forces of respective supercells. 409 410 dataset : dict 411 The format can be either one of two types 412 413 Type 1. One atomic displacement in each supercell: 414 {'natom': number of atoms in supercell, 415 'first_atoms': [ 416 {'number': atom index of displaced atom, 417 'displacement': displacement in Cartesian coordinates, 418 'forces': forces on atoms in supercell}, 419 {...}, ...]} 420 Elements of the list accessed by 'first_atoms' corresponds to each 421 displaced supercell. Each displaced supercell contains only one 422 displacement. dict['first_atoms']['forces'] gives atomic forces in 423 each displaced supercell. 424 425 Type 2. All atomic displacements in each supercell: 426 {'displacements': ndarray, dtype='double', order='C', 427 shape=(supercells, natom, 3) 428 'forces': ndarray, dtype='double',, order='C', 429 shape=(supercells, natom, 3)} 430 431 To set in type 2, displacements and forces can be given by numpy 432 array with different shape but that can be reshaped to 433 (supercells, natom, 3). 434 435 """ 436 return self._displacement_dataset 437 438 @property 439 def displacement_dataset(self): 440 """Return dataset to store displacements and forces.""" 441 warnings.warn("Phonopy.displacement_dataset attribute is deprecated." 442 "Use Phonopy.dataset attribute.", 443 DeprecationWarning) 444 return self.dataset 445 446 def get_displacement_dataset(self): 447 """Return dataset to store displacements and forces.""" 448 warnings.warn("Phonopy.get_displacement_dataset() is deprecated." 449 "Use Phonopy.dataset attribute.", 450 DeprecationWarning) 451 return self.dataset 452 453 @dataset.setter 454 def dataset(self, dataset): 455 if dataset is None: 456 self._displacement_dataset = None 457 elif 'first_atoms' in dataset: 458 self._displacement_dataset = dataset 459 elif 'displacements' in dataset: 460 self._displacement_dataset = {} 461 self.displacements = dataset['displacements'] 462 if 'forces' in dataset: 463 self.forces = dataset['forces'] 464 else: 465 raise RuntimeError("Data format of dataset is wrong.") 466 467 self._supercells_with_displacements = None 468 469 def set_displacement_dataset(self, displacement_dataset): 470 """Set displacements.""" 471 warnings.warn("Phonopy.set_displacement_dataset() is deprecated." 472 "Use Phonopy.dataset attribute.", 473 DeprecationWarning) 474 self.dataset = displacement_dataset 475 476 @property 477 def displacements(self): 478 """Getter and setter of displacements in supercells. 479 480 There are two types of displacement dataset. See the docstring 481 of dataset about types 1 and 2 for the displacement dataset formats. 482 Displacements set returned depends on either type-1 or type-2 as 483 follows: 484 485 Type-1, List of list 486 The internal list has 4 elements such as [32, 0.01, 0.0, 0.0]]. 487 The first element is the supercell atom index starting with 0. 488 The remaining three elements give the displacement in Cartesian 489 coordinates. 490 Type-2, array_like 491 Displacements of all atoms of all supercells in Cartesian 492 coordinates. 493 shape=(supercells, natom, 3) 494 dtype='double' 495 496 497 To set displacements set, only type-2 datast case is allowed. 498 499 displacemens : array_like 500 Atomic displacements of all atoms of all supercells. 501 Only all displacements in each supercell case (type-2) is 502 supported. 503 shape=(supercells, natom, 3), dtype='double', order='C' 504 505 """ 506 disps = [] 507 if 'first_atoms' in self._displacement_dataset: 508 for disp in self._displacement_dataset['first_atoms']: 509 x = disp['displacement'] 510 disps.append([disp['number'], x[0], x[1], x[2]]) 511 elif 'displacements' in self._displacement_dataset: 512 disps = self._displacement_dataset['displacements'] 513 514 return disps 515 516 def get_displacements(self): 517 """Return displacements in supercells.""" 518 warnings.warn("Phonopy.get_displacements() is deprecated." 519 "Use Phonopy.displacements attribute.", 520 DeprecationWarning) 521 return self.displacements 522 523 @displacements.setter 524 def displacements(self, displacements): 525 disp = np.array(displacements, dtype='double', order='C') 526 if disp.ndim != 3 or disp.shape[1:] != (len(self._supercell), 3): 527 raise RuntimeError("Array shape of displacements is incorrect.") 528 529 if 'first_atoms' in self._displacement_dataset: 530 raise RuntimeError("This displacement format is not supported.") 531 532 self._displacement_dataset['displacements'] = disp 533 534 @property 535 def force_constants(self): 536 """Getter and setter of supercell force constants. 537 538 Force constants matrix. 539 540 ndarray to get 541 There are two shapes: 542 full: 543 shape=(atoms in supercell, atoms in supercell, 3, 3) 544 compact: 545 shape=(atoms in primitive cell, atoms in supercell, 3, 3) 546 dtype='double', order='C' 547 548 array_like to set 549 If this is given in own condiguous ndarray with order='C' and 550 dtype='double', internal copy of data is avoided. Therefore 551 some computational resources are saved. 552 shape=(atoms in supercell, atoms in supercell, 3, 3), 553 dtype='double' 554 555 """ 556 return self._force_constants 557 558 def get_force_constants(self): 559 """Return supercell force constants.""" 560 warnings.warn("Phonopy.get_force_constants() is deprecated." 561 "Use Phonopy.force_constants attribute.", 562 DeprecationWarning) 563 return self.force_constants 564 565 @force_constants.setter 566 def force_constants(self, force_constants): 567 if type(force_constants) is np.ndarray: 568 fc_shape = force_constants.shape 569 if fc_shape[0] != fc_shape[1]: 570 if len(self._primitive) != fc_shape[0]: 571 msg = ("Force constants shape disagrees with crystal " 572 "structure setting. This may be due to " 573 "PRIMITIVE_AXIS.") 574 raise RuntimeError(msg) 575 576 self._force_constants = force_constants 577 if self._primitive.masses is not None: 578 self._set_dynamical_matrix() 579 580 def set_force_constants(self, force_constants, show_drift=True): 581 """Set force constants.""" 582 warnings.warn("Phonopy.set_force_constants() is deprecated." 583 "Use Phonopy.force_constants attribute.", 584 DeprecationWarning) 585 self.force_constants = force_constants 586 if show_drift and self._log_level: 587 show_drift_force_constants(self._force_constants, 588 primitive=self._primitive) 589 590 def set_force_constants_zero_with_radius(self, cutoff_radius): 591 """Set zero to force constants within cutoff radius.""" 592 cutoff_force_constants(self._force_constants, 593 self._supercell, 594 self._primitive, 595 cutoff_radius, 596 symprec=self._symprec) 597 if self._primitive.masses is not None: 598 self._set_dynamical_matrix() 599 600 @property 601 def forces(self): 602 """Return forces of supercells. 603 604 ndarray to get and array_like to set 605 A set of atomic forces in displaced supercells. The order of 606 displaced supercells has to match with that in displacement 607 dataset. 608 shape=(supercells with displacements, atoms in supercell, 3) 609 dtype='double', order='C' 610 611 [[[f_1x, f_1y, f_1z], [f_2x, f_2y, f_2z], ...], # first supercell 612 [[f_1x, f_1y, f_1z], [f_2x, f_2y, f_2z], ...], # second supercell 613 ... 614 ] 615 616 """ 617 if 'forces' in self._displacement_dataset: 618 return self._displacement_dataset['forces'] 619 elif 'first_atoms' in self._displacement_dataset: 620 forces = [] 621 for disp in self._displacement_dataset['first_atoms']: 622 if 'forces' in disp: 623 forces.append(disp['forces']) 624 if forces: 625 return np.array(forces, dtype='double', order='C') 626 else: 627 None 628 else: 629 return None 630 631 @forces.setter 632 def forces(self, sets_of_forces): 633 if 'first_atoms' in self._displacement_dataset: 634 for disp, forces in zip(self._displacement_dataset['first_atoms'], 635 sets_of_forces): 636 disp['forces'] = forces 637 elif 'displacements' in self._displacement_dataset: 638 forces = np.array(sets_of_forces, dtype='double', order='C') 639 natom = len(self._supercell) 640 if forces.ndim != 3 or forces.shape[1:] != (natom, 3): 641 raise RuntimeError("Array shape of input forces is incorrect.") 642 643 self._displacement_dataset['forces'] = forces 644 645 def set_forces(self, sets_of_forces): 646 """Set forces of supercells.""" 647 warnings.warn("Phonopy.set_forces() is deprecated." 648 "Use Phonopy.forces attribute.", 649 DeprecationWarning) 650 self.forces = sets_of_forces 651 652 @property 653 def dynamical_matrix(self): 654 """Return DynamicalMatrix instance. 655 656 This is not dynamical matrices but the instance of DynamicalMatrix 657 class. 658 659 """ 660 return self._dynamical_matrix 661 662 def get_dynamical_matrix(self): 663 """Return DynamicalMatrix instance.""" 664 warnings.warn("Phonopy.get_dynamical_matrix() is deprecated." 665 "Use Phonopy.dynamical_matrix attribute.", 666 DeprecationWarning) 667 return self.dynamical_matrix 668 669 @property 670 def nac_params(self): 671 """Getter and setter of parameters for non-analytical term correction. 672 673 dict 674 Parameters used for non-analytical term correction 675 'born': ndarray 676 Born effective charges 677 shape=(primitive cell atoms, 3, 3), dtype='double', order='C' 678 'factor': float 679 Unit conversion factor 680 'dielectric': ndarray 681 Dielectric constant tensor 682 shape=(3, 3), dtype='double', order='C' 683 684 """ 685 return self._nac_params 686 687 def get_nac_params(self): 688 """Return parameters for non-analytical term correction.""" 689 warnings.warn("Phonopy.get_nac_params() is deprecated." 690 "Use Phonopy.nac_params attribute.", 691 DeprecationWarning) 692 return self.nac_params 693 694 @nac_params.setter 695 def nac_params(self, nac_params): 696 self._nac_params = nac_params 697 if self._force_constants is not None: 698 self._set_dynamical_matrix() 699 700 def set_nac_params(self, nac_params): 701 """Set parameters for non-analytical term correction.""" 702 warnings.warn("Phonopy.set_nac_params() is deprecated." 703 "Use Phonopy.nac_params attribute.", 704 DeprecationWarning) 705 self.nac_params = nac_params 706 707 @property 708 def supercells_with_displacements(self): 709 """Return supercells with displacements. 710 711 list of PhonopyAtoms 712 Supercells with displacements generated by 713 Phonopy.generate_displacements. 714 715 """ 716 if self._displacement_dataset is None: 717 return None 718 else: 719 if self._supercells_with_displacements is None: 720 self._build_supercells_with_displacements() 721 return self._supercells_with_displacements 722 723 def get_supercells_with_displacements(self): 724 """Return supercells with displacements.""" 725 warnings.warn( 726 "Phonopy.get_supercells_with_displacements() is deprecated." 727 "Use Phonopy.supercells_with_displacements attribute.", 728 DeprecationWarning) 729 return self.supercells_with_displacements 730 731 @property 732 def mesh_numbers(self): 733 """Return sampling mesh numbers in reciprocal space.""" 734 if self._mesh is None: 735 return None 736 else: 737 return self._mesh.mesh_numbers 738 739 @property 740 def qpoints(self): 741 """Return QpointsPhonon instance.""" 742 return self._qpoints 743 744 @property 745 def band_structure(self): 746 """Return BandStructure instance.""" 747 return self._band_structure 748 749 @property 750 def mesh(self): 751 """Return Mesh or IterMesh instance.""" 752 return self._mesh 753 754 @property 755 def random_displacements(self): 756 """Return RandomDisplacements instance.""" 757 return self._random_displacements 758 759 @property 760 def dynamic_structure_factor(self): 761 """Return DynamicStructureFactor instance.""" 762 return self._dynamic_structure_factor 763 764 @property 765 def thermal_properties(self): 766 """Return ThermalProperties instance.""" 767 return self._thermal_properties 768 769 @property 770 def thermal_displacements(self): 771 """Return ThermalDisplacements instance.""" 772 return self._thermal_displacements 773 774 @property 775 def thermal_displacement_matrices(self): 776 """Return ThermalDisplacementMatrices instance.""" 777 return self._thermal_displacement_matrices 778 779 @property 780 def irreps(self): 781 """Return IrReps instance.""" 782 return self._irreps 783 784 @property 785 def moment(self): 786 """Return PhononMoment instance.""" 787 return self._moment 788 789 @property 790 def total_dos(self): 791 """Return TotalDos instance.""" 792 return self._total_dos 793 794 @property 795 def partial_dos(self): 796 """Return PartialDos instance.""" 797 warnings.warn("Phonopy.partial_dos is deprecated." 798 "Use Phonopy.projected_dos.", 799 DeprecationWarning) 800 return self.projected_dos 801 802 @property 803 def projected_dos(self): 804 """Return PartialDos instance.""" 805 return self._pdos 806 807 @property 808 def masses(self): 809 """Getter and setter of masses of primitive cell atoms.""" 810 return self._primitive.masses 811 812 @masses.setter 813 def masses(self, masses): 814 p_masses = np.array(masses) 815 self._primitive.set_masses(p_masses) 816 p2p_map = self._primitive.p2p_map 817 s_masses = p_masses[[p2p_map[x] for x in self._primitive.s2p_map]] 818 self._supercell.set_masses(s_masses) 819 u2s_map = self._supercell.u2s_map 820 u_masses = s_masses[u2s_map] 821 self._unitcell.set_masses(u_masses) 822 if self._force_constants is not None: 823 self._set_dynamical_matrix() 824 825 def set_masses(self, masses): 826 """Set masses of primitive cell atoms.""" 827 self.masses = masses 828 829 def generate_displacements(self, 830 distance=0.01, 831 is_plusminus='auto', 832 is_diagonal=True, 833 is_trigonal=False, 834 number_of_snapshots=None, 835 random_seed=None, 836 temperature=None, 837 cutoff_frequency=None): 838 """Generate displacement dataset. 839 840 There are two modes, finite difference method with systematic 841 displacements and fitting approach between arbitrary displacements and 842 their forces. The default approach is the finite difference method that 843 is built-in phonopy. The fitting approach requires external force 844 constant calculator. 845 846 The random displacement supercells are created by setting positive 847 integer values 'number_of_snapshots' keyword argument. Unless 848 this is specified, systematic displacements are created for the finite 849 difference method as the default behaviour. 850 851 Parameters 852 ---------- 853 distance : float, optional 854 Displacement distance. Unit is the same as that used for crystal 855 structure. Default is 0.01. 856 is_plusminus : 'auto', True, or False, optional 857 For each atom, displacement of one direction (False), both 858 direction, i.e., one directiona and its opposite direction 859 (True), and both direction if symmetry requires ('auto'). 860 Default is 'auto'. 861 is_diagonal : bool, optional 862 Displacements are made only along basis vectors (False) and 863 can be made not being along basis vectors if the number of 864 displacements can be reduced by symmetry (True). Default is 865 True. 866 is_trigonal : bool, optional 867 Existing only testing purpose. 868 number_of_snapshots : int or None, optional 869 Number of snapshots of supercells with random displacements. 870 Random displacements are generated displacing all atoms in 871 random directions with a fixed displacement distance specified 872 by 'distance' parameter, i.e., all atoms in supercell are 873 displaced with the same displacement distance in direct space. 874 random_seed : 32bit unsigned int or None, optional 875 Random seed for random displacements generation. 876 temperature : float 877 With given temperature, random displacements at temperature is 878 generated by sampling probability distribution from canonical 879 ensemble of harmonic oscillators (harmonic phonons). 880 cutoff_frequency : float 881 In random displacements generation from canonical ensemble 882 of harmonic phonons, phonon occupation number is used to 883 determine the deviation of the distribution function. 884 To avoid too large deviation, this value is used to exclude 885 the phonon modes whose absolute frequency are smaller than 886 this value. 887 888 """ 889 if (np.issubdtype(type(number_of_snapshots), np.integer) and 890 number_of_snapshots > 0): # noqa: E129 891 if temperature is None: 892 displacement_dataset = get_random_displacements_dataset( 893 number_of_snapshots, 894 distance, 895 len(self._supercell), 896 random_seed=random_seed) 897 else: 898 self.run_random_displacements( 899 temperature, 900 number_of_snapshots=number_of_snapshots, 901 random_seed=random_seed, 902 cutoff_frequency=cutoff_frequency) 903 units = get_default_physical_units(self._calculator) 904 d = np.array( 905 self._random_displacements.u / units['distance_to_A'], 906 dtype='double', order='C') 907 displacement_dataset = {'displacements': d} 908 else: 909 displacement_directions = get_least_displacements( 910 self._symmetry, 911 is_plusminus=is_plusminus, 912 is_diagonal=is_diagonal, 913 is_trigonal=is_trigonal, 914 log_level=self._log_level) 915 displacement_dataset = directions_to_displacement_dataset( 916 displacement_directions, 917 distance, 918 self._supercell) 919 self.dataset = displacement_dataset 920 921 def produce_force_constants(self, 922 forces=None, 923 calculate_full_force_constants=True, 924 fc_calculator=None, 925 fc_calculator_options=None, 926 show_drift=True): 927 """Compute supercell force constants from forces-displacements dataset. 928 929 Supercell force constants are computed from forces and displacements. 930 As the default behaviour, those stored in dataset are used. But 931 with setting ``forces``, this set of forces and the set of 932 displacements stored in the dataset are used for the computation. 933 934 Parameters 935 ---------- 936 forces : array_like, optional 937 See docstring of Phonopy.forces. Default is None. 938 calculate_full_force_constants : Bool, optional 939 With setting True, full force constants matrix is stored. 940 With setting False, compact force constants matrix is stored. 941 For more detail, see docstring of Phonopy.force_constants. 942 Default is True. 943 fc_calculator : str, optional 944 fc_calculator_options : str, optional 945 External force constants calculator is used. Currently, 946 'alm' is supported. See more detail at the docstring of 947 phonopy.interface.fc_calculator.get_fc2. Default is None. 948 show_drift : Bool, optional 949 With setting 950 951 """ 952 if forces is not None: 953 self.forces = forces 954 955 # A primitive check if 'forces' key is in displacement_dataset. 956 if 'first_atoms' in self._displacement_dataset: 957 for disp in self._displacement_dataset['first_atoms']: 958 if 'forces' not in disp: 959 raise RuntimeError("Forces are not yet set.") 960 elif 'forces' not in self._displacement_dataset: 961 raise RuntimeError("Forces are not yet set.") 962 963 if calculate_full_force_constants: 964 self._run_force_constants_from_forces( 965 fc_calculator=fc_calculator, 966 fc_calculator_options=fc_calculator_options, 967 decimals=self._force_constants_decimals) 968 else: 969 self._run_force_constants_from_forces( 970 distributed_atom_list=self._primitive.p2s_map, 971 fc_calculator=fc_calculator, 972 fc_calculator_options=fc_calculator_options, 973 decimals=self._force_constants_decimals) 974 975 if show_drift and self._log_level: 976 show_drift_force_constants(self._force_constants, 977 primitive=self._primitive) 978 979 if self._primitive.masses is not None: 980 self._set_dynamical_matrix() 981 982 def symmetrize_force_constants(self, level=1, show_drift=True): 983 """Symmetrize force constants. 984 985 This applies translational and permutation symmetries successfully, 986 but not simultaneously. 987 988 Parameters 989 ---------- 990 level : int, optional 991 Application of translational and permulation symmetries is 992 repeated by this number. Default is 1. 993 show_drift : bool, optioanl 994 Drift forces are displayed when True. Default is True. 995 996 """ 997 if self._force_constants is None: 998 raise RuntimeError("Force constants have not been produced yet.") 999 1000 if self._force_constants.shape[0] == self._force_constants.shape[1]: 1001 symmetrize_force_constants(self._force_constants, level=level) 1002 else: 1003 symmetrize_compact_force_constants(self._force_constants, 1004 self._primitive, 1005 level=level) 1006 if show_drift and self._log_level: 1007 sys.stdout.write("Max drift after symmetrization by translation: ") 1008 show_drift_force_constants(self._force_constants, 1009 primitive=self._primitive, 1010 values_only=True) 1011 1012 if self._primitive.masses is not None: 1013 self._set_dynamical_matrix() 1014 1015 def symmetrize_force_constants_by_space_group(self, show_drift=True): 1016 """Symmetrize force constants using space group operations. 1017 1018 Space group operations except for pure translations are applied 1019 to force constants. 1020 1021 Parameters 1022 ---------- 1023 show_drift : bool, optioanl 1024 Drift forces are displayed when True. Default is True. 1025 1026 """ 1027 set_tensor_symmetry_PJ(self._force_constants, 1028 self._supercell.cell.T, 1029 self._supercell.scaled_positions, 1030 self._symmetry) 1031 1032 if show_drift and self._log_level: 1033 sys.stdout.write("Max drift after symmetrization by space group: ") 1034 show_drift_force_constants(self._force_constants, 1035 primitive=self._primitive, 1036 values_only=True) 1037 1038 if self._primitive.masses is not None: 1039 self._set_dynamical_matrix() 1040 1041 ##################### 1042 # Phonon properties # 1043 ##################### 1044 1045 # Single q-point 1046 def get_dynamical_matrix_at_q(self, q): 1047 """Calculate dynamical matrix at a given q-point. 1048 1049 Parameters 1050 ---------- 1051 q: array_like 1052 A q-vector. 1053 shape=(3,), dtype='double' 1054 1055 Returns 1056 ------- 1057 dynamical_matrix: ndarray 1058 Dynamical matrix. 1059 shape=(bands, bands) 1060 dtype=complex of "c%d" % (np.dtype('double').itemsize * 2) 1061 order='C' 1062 1063 """ 1064 self._set_dynamical_matrix() 1065 if self._dynamical_matrix is None: 1066 msg = ("Dynamical matrix has not yet built.") 1067 raise RuntimeError(msg) 1068 1069 self._dynamical_matrix.run(q) 1070 return self._dynamical_matrix.get_dynamical_matrix() 1071 1072 def get_frequencies(self, q): 1073 """Calculate phonon frequencies at a given q-point. 1074 1075 Parameters 1076 ---------- 1077 q: array_like 1078 A q-vector. 1079 shape=(3,), dtype='double' 1080 1081 Returns 1082 ------- 1083 frequencies: ndarray 1084 Phonon frequencies. Imaginary frequenies are represented by 1085 negative real numbers. 1086 shape=(bands, ), dtype='double' 1087 1088 """ 1089 self._set_dynamical_matrix() 1090 if self._dynamical_matrix is None: 1091 msg = ("Dynamical matrix has not yet built.") 1092 raise RuntimeError(msg) 1093 1094 self._dynamical_matrix.run(q) 1095 dm = self._dynamical_matrix.get_dynamical_matrix() 1096 frequencies = [] 1097 for eig in np.linalg.eigvalsh(dm).real: 1098 if eig < 0: 1099 frequencies.append(-np.sqrt(-eig)) 1100 else: 1101 frequencies.append(np.sqrt(eig)) 1102 1103 return np.array(frequencies) * self._factor 1104 1105 def get_frequencies_with_eigenvectors(self, q): 1106 """Calculate phonon frequencies and eigenvectors at a given q-point. 1107 1108 Parameters 1109 ---------- 1110 q: array_like 1111 A q-vector. 1112 shape=(3,) 1113 1114 Returns 1115 ------- 1116 (frequencies, eigenvectors) 1117 1118 frequencies: ndarray 1119 Phonon frequencies. Imaginary frequenies are represented by 1120 negative real numbers. 1121 shape=(bands, ), dtype='double', order='C' 1122 eigenvectors: ndarray 1123 Phonon eigenvectors. 1124 shape=(bands, bands) 1125 dtype=complex of "c%d" % (np.dtype('double').itemsize * 2) 1126 order='C' 1127 1128 """ 1129 self._set_dynamical_matrix() 1130 if self._dynamical_matrix is None: 1131 msg = ("Dynamical matrix has not yet built.") 1132 raise RuntimeError(msg) 1133 1134 self._dynamical_matrix.run(q) 1135 dm = self._dynamical_matrix.get_dynamical_matrix() 1136 frequencies = [] 1137 eigvals, eigenvectors = np.linalg.eigh(dm) 1138 frequencies = [] 1139 for eig in eigvals: 1140 if eig < 0: 1141 frequencies.append(-np.sqrt(-eig)) 1142 else: 1143 frequencies.append(np.sqrt(eig)) 1144 1145 return np.array(frequencies) * self._factor, eigenvectors 1146 1147 # Band structure 1148 def run_band_structure(self, 1149 paths, 1150 with_eigenvectors=False, 1151 with_group_velocities=False, 1152 is_band_connection=False, 1153 path_connections=None, 1154 labels=None, 1155 is_legacy_plot=False): 1156 """Run phonon band structure calculation. 1157 1158 Parameters 1159 ---------- 1160 paths : List of array_like 1161 Sets of qpoints that can be passed to phonopy.set_band_structure(). 1162 Numbers of qpoints can be different. 1163 shape of each array_like : (qpoints, 3) 1164 with_eigenvectors : bool, optional 1165 Flag whether eigenvectors are calculated or not. Default is False. 1166 with_group_velocities : bool, optional 1167 Flag whether group velocities are calculated or not. Default is 1168 False. 1169 is_band_connection : bool, optional 1170 Flag whether each band is connected or not. This is achieved by 1171 comparing similarity of eigenvectors of neghboring poins. Sometimes 1172 this fails. Default is False. 1173 path_connections : List of bool, optional 1174 This is only used in graphical plot of band structure and gives 1175 whether each path is connected to the next path or not, 1176 i.e., if False, there is a jump of q-points. Number of elements is 1177 the same at that of paths. Default is None. 1178 labels : List of str, optional 1179 This is only used in graphical plot of band structure and gives 1180 labels of end points of each path. The number of labels is equal 1181 to (2 - np.array(path_connections)).sum(). 1182 is_legacy_plot: bool, optional 1183 This makes the old style band structure plot. Default is False. 1184 1185 """ 1186 if self._dynamical_matrix is None: 1187 msg = ("Dynamical matrix has not yet built.") 1188 raise RuntimeError(msg) 1189 1190 if with_group_velocities: 1191 if self._group_velocity is None: 1192 self._set_group_velocity() 1193 group_velocity = self._group_velocity 1194 else: 1195 group_velocity = None 1196 1197 self._band_structure = BandStructure( 1198 paths, 1199 self._dynamical_matrix, 1200 with_eigenvectors=with_eigenvectors, 1201 is_band_connection=is_band_connection, 1202 group_velocity=group_velocity, 1203 path_connections=path_connections, 1204 labels=labels, 1205 is_legacy_plot=is_legacy_plot, 1206 factor=self._factor) 1207 1208 def set_band_structure(self, 1209 bands, 1210 is_eigenvectors=False, 1211 is_band_connection=False, 1212 path_connections=None, 1213 labels=None, 1214 is_legacy_plot=False): 1215 """Calculate phonon band structure.""" 1216 warnings.warn("Phonopy.set_band_structure() is deprecated. " 1217 "Use Phonopy.run_band_structure().", DeprecationWarning) 1218 1219 if self._group_velocity is None: 1220 with_group_velocities = False 1221 else: 1222 with_group_velocities = True 1223 self.run_band_structure(bands, 1224 with_eigenvectors=is_eigenvectors, 1225 with_group_velocities=with_group_velocities, 1226 is_band_connection=is_band_connection, 1227 path_connections=path_connections, 1228 labels=labels, 1229 is_legacy_plot=is_legacy_plot) 1230 1231 def get_band_structure_dict(self): 1232 """Return calculated band structures. 1233 1234 Returns 1235 ------- 1236 dict 1237 keys: qpoints, distances, frequencies, eigenvectors, and 1238 group_velocities 1239 Each dict value is a list containing properties on number of paths. 1240 The number of q-points on one path can be different from that on 1241 the other path. Each set of properties on a path is ndarray and is 1242 explained as below: 1243 1244 qpoints[i]: ndarray 1245 q-points in reduced coordinates of reciprocal space without 1246 2pi. 1247 shape=(q-points, 3), dtype='double' 1248 distances[i]: ndarray 1249 Distances in reciprocal space along paths. 1250 shape=(q-points,), dtype='double' 1251 frequencies[i]: ndarray 1252 Phonon frequencies. Imaginary frequenies are represented by 1253 negative real numbers. 1254 shape=(q-points, bands), dtype='double' 1255 eigenvectors[i]: ndarray 1256 Phonon eigenvectors. None if eigenvectors are not stored. 1257 shape=(q-points, bands, bands) 1258 dtype=complex of "c%d" % (np.dtype('double').itemsize * 2) 1259 order='C' 1260 group_velocities[i]: ndarray 1261 Phonon group velocities. None if group velocities are not 1262 calculated. 1263 shape=(q-points, bands, 3), dtype='double' 1264 1265 """ 1266 if self._band_structure is None: 1267 msg = ("Phonopy.run_band_structure() has to be done.") 1268 raise RuntimeError(msg) 1269 1270 retdict = {'qpoints': self._band_structure.qpoints, 1271 'distances': self._band_structure.distances, 1272 'frequencies': self._band_structure.frequencies, 1273 'eigenvectors': self._band_structure.eigenvectors, 1274 'group_velocities': self._band_structure.group_velocities} 1275 1276 return retdict 1277 1278 def get_band_structure(self): 1279 """Return calculated band structures. 1280 1281 Returns 1282 ------- 1283 (q-points, distances, frequencies, eigenvectors) 1284 1285 Each tuple element is a list containing properties on number of paths. 1286 The number of q-points on one path can be different from that on the 1287 other path. Each set of properties on a path is ndarray and is 1288 explained as below: 1289 1290 q-points[i]: ndarray 1291 q-points in reduced coordinates of reciprocal space without 2pi. 1292 shape=(q-points, 3), dtype='double' 1293 distances[i]: ndarray 1294 Distances in reciprocal space along paths. 1295 shape=(q-points,), dtype='double' 1296 frequencies[i]: ndarray 1297 Phonon frequencies. Imaginary frequenies are represented by 1298 negative real numbers. 1299 shape=(q-points, bands), dtype='double' 1300 eigenvectors[i]: ndarray 1301 Phonon eigenvectors. None if eigenvectors are not stored. 1302 shape=(q-points, bands, bands) 1303 dtype=complex of "c%d" % (np.dtype('double').itemsize * 2) 1304 order='C' 1305 group_velocities[i]: ndarray 1306 Phonon group velocities. None if group velocities are not 1307 calculated. 1308 shape=(q-points, bands, 3), dtype='double' 1309 1310 """ 1311 warnings.warn("Phonopy.get_band_structure() is deprecated. " 1312 "Use Phonopy.get_band_structure_dict().", 1313 DeprecationWarning) 1314 1315 if self._band_structure is None: 1316 msg = ("run_band_structure has to be done.") 1317 raise RuntimeError(msg) 1318 1319 retvals = (self._band_structure.qpoints, 1320 self._band_structure.distances, 1321 self._band_structure.frequencies, 1322 self._band_structure.eigenvectors) 1323 return retvals 1324 1325 def auto_band_structure(self, 1326 npoints=101, 1327 with_eigenvectors=False, 1328 with_group_velocities=False, 1329 plot=False, 1330 write_yaml=False, 1331 filename="band.yaml"): 1332 """Conveniently calculate and draw band structure. 1333 1334 Parameters 1335 ---------- 1336 See docstring of ``Phonopy.run_band_structure`` for the parameters of 1337 ``with_eigenvectors`` (default is False) and ``with_group_velocities`` 1338 (default is False). 1339 1340 npoints : int, optional 1341 Number of q-points in each segment of band struture paths. 1342 The number includes end points. Default is 101. 1343 plot : Bool, optional 1344 With setting True, band structure is plotted using matplotlib and 1345 the matplotlib module (plt) is returned. To watch the result, 1346 usually ``show()`` has to be called. Default is False. 1347 write_yaml : Bool 1348 With setting True, ``band.yaml`` like file is written out. The 1349 file name can be specified with the ``filename`` parameter. 1350 Default is False. 1351 filename : str, optional 1352 File name used to write ``band.yaml`` like file. Default is 1353 ``band.yaml``. 1354 1355 """ 1356 bands, labels, path_connections = get_band_qpoints_by_seekpath( 1357 self._primitive, npoints, is_const_interval=True) 1358 self.run_band_structure(bands, 1359 with_eigenvectors=with_eigenvectors, 1360 with_group_velocities=with_group_velocities, 1361 path_connections=path_connections, 1362 labels=labels, 1363 is_legacy_plot=False) 1364 if write_yaml: 1365 self.write_yaml_band_structure(filename=filename) 1366 if plot: 1367 return self.plot_band_structure() 1368 1369 def plot_band_structure(self): 1370 """Plot calculated band structure. 1371 1372 Returns 1373 ------- 1374 matplotlib.pyplot. 1375 1376 """ 1377 import matplotlib.pyplot as plt 1378 1379 if self._band_structure.labels: 1380 from matplotlib import rc 1381 rc('text', usetex=True) 1382 1383 if self._band_structure.is_legacy_plot: 1384 fig, axs = plt.subplots(1, 1) 1385 else: 1386 from mpl_toolkits.axes_grid1 import ImageGrid 1387 n = len([x for x in self._band_structure.path_connections 1388 if not x]) 1389 fig = plt.figure() 1390 axs = ImageGrid(fig, 111, # similar to subplot(111) 1391 nrows_ncols=(1, n), 1392 axes_pad=0.11, 1393 label_mode="L") 1394 self._band_structure.plot(axs) 1395 return plt 1396 1397 def write_hdf5_band_structure(self, 1398 comment=None, 1399 filename="band.hdf5"): 1400 """Write band structure in hdf5 format. 1401 1402 Parameters 1403 ---------- 1404 comment : dict, optional 1405 Items are stored in hdf5 file in the way of key-value pair. 1406 filename : str, optional 1407 Default is ``band.hdf5``. 1408 1409 """ 1410 self._band_structure.write_hdf5(comment=comment, filename=filename) 1411 1412 def write_yaml_band_structure(self, 1413 comment=None, 1414 filename=None, 1415 compression=None): 1416 """Write band structure in yaml. 1417 1418 Parameters 1419 ---------- 1420 comment : dict 1421 Data structure dumped in YAML and the dumped YAML text is put 1422 at the beggining of the file. 1423 filename : str 1424 Default filename is 'band.yaml' when compression=None. 1425 With compression, an extention of filename is added such as 1426 'band.yaml.xz'. 1427 compression : None, 'gzip', or 'lzma' 1428 None gives usual text file. 'gzip and 'lzma' compresse yaml 1429 text in respective compression methods. 1430 1431 """ 1432 self._band_structure.write_yaml(comment=comment, 1433 filename=filename, 1434 compression=compression) 1435 1436 def init_mesh(self, 1437 mesh=100.0, 1438 shift=None, 1439 is_time_reversal=True, 1440 is_mesh_symmetry=True, 1441 with_eigenvectors=False, 1442 with_group_velocities=False, 1443 is_gamma_center=False, 1444 use_iter_mesh=False): 1445 """Initialize mesh sampling phonon calculation without starting to run. 1446 1447 Phonon calculation starts explicitly with calling Mesh.run() or 1448 implicitly with accessing getters of Mesh instance, e.g., 1449 Mesh.frequencies. 1450 1451 Parameters 1452 ---------- 1453 mesh: array_like or float, optional 1454 Mesh numbers along a, b, c axes when array_like object is given. 1455 dtype='intc', shape=(3,) 1456 When float value is given, uniform mesh is generated following 1457 VASP convention by 1458 N = max(1, nint(l * |a|^*)) 1459 where 'nint' is the function to return the nearest integer. In this 1460 case, it is forced to set is_gamma_center=True. 1461 Default value is 100.0. 1462 shift: array_like, optional 1463 Mesh shifts along a*, b*, c* axes with respect to neighboring grid 1464 points from the original mesh (Monkhorst-Pack or Gamma center). 1465 0.5 gives half grid shift. Normally 0 or 0.5 is given. 1466 Otherwise q-points symmetry search is not performed. 1467 Default is None (no additional shift). 1468 dtype='double', shape=(3, ) 1469 is_time_reversal: bool, optional 1470 Time reversal symmetry is considered in symmetry search. By this, 1471 inversion symmetry is always included. Default is True. 1472 is_mesh_symmetry: bool, optional 1473 Wheather symmetry search is done or not. Default is True 1474 with_eigenvectors: bool, optional 1475 Eigenvectors are stored by setting True. Default False. 1476 with_group_velocities : bool, optional 1477 Group velocities are calculated by setting True. Default is 1478 False. 1479 is_gamma_center: bool, default False 1480 Uniform mesh grids are generated centring at Gamma point but not 1481 the Monkhorst-Pack scheme. When type(mesh) is float, this parameter 1482 setting is ignored and it is forced to set is_gamma_center=True. 1483 use_iter_mesh: bool 1484 Use IterMesh instead of Mesh class not to store phonon properties 1485 in its instance to save memory consumption. This is used with 1486 ThermalDisplacements and ThermalDisplacementMatrices. 1487 Default is False. 1488 1489 """ 1490 if self._dynamical_matrix is None: 1491 msg = "Dynamical matrix has not yet built." 1492 raise RuntimeError(msg) 1493 1494 _mesh = np.array(mesh) 1495 mesh_nums = None 1496 if _mesh.shape: 1497 if _mesh.shape == (3,): 1498 mesh_nums = mesh 1499 _is_gamma_center = is_gamma_center 1500 else: 1501 if self._primitive_symmetry is not None: 1502 rots = self._primitive_symmetry.get_pointgroup_operations() 1503 mesh_nums = length2mesh(mesh, 1504 self._primitive.cell, 1505 rotations=rots) 1506 else: 1507 mesh_nums = length2mesh(mesh, self._primitive.cell) 1508 _is_gamma_center = True 1509 if mesh_nums is None: 1510 msg = "mesh has inappropriate type." 1511 raise TypeError(msg) 1512 1513 if with_group_velocities: 1514 if self._group_velocity is None: 1515 self._set_group_velocity() 1516 group_velocity = self._group_velocity 1517 else: 1518 group_velocity = None 1519 1520 if use_iter_mesh: 1521 self._mesh = IterMesh( 1522 self._dynamical_matrix, 1523 mesh_nums, 1524 shift=shift, 1525 is_time_reversal=is_time_reversal, 1526 is_mesh_symmetry=is_mesh_symmetry, 1527 with_eigenvectors=with_eigenvectors, 1528 is_gamma_center=is_gamma_center, 1529 rotations=self._primitive_symmetry.get_pointgroup_operations(), 1530 factor=self._factor) 1531 else: 1532 self._mesh = Mesh( 1533 self._dynamical_matrix, 1534 mesh_nums, 1535 shift=shift, 1536 is_time_reversal=is_time_reversal, 1537 is_mesh_symmetry=is_mesh_symmetry, 1538 with_eigenvectors=with_eigenvectors, 1539 is_gamma_center=_is_gamma_center, 1540 group_velocity=group_velocity, 1541 rotations=self._primitive_symmetry.get_pointgroup_operations(), 1542 factor=self._factor, 1543 use_lapack_solver=self._use_lapack_solver) 1544 1545 def run_mesh(self, 1546 mesh=100.0, 1547 shift=None, 1548 is_time_reversal=True, 1549 is_mesh_symmetry=True, 1550 with_eigenvectors=False, 1551 with_group_velocities=False, 1552 is_gamma_center=False): 1553 """Run mesh sampling phonon calculation. 1554 1555 See the parameter details in Phonopy.init_mesh. 1556 1557 """ 1558 self.init_mesh(mesh=mesh, 1559 shift=shift, 1560 is_time_reversal=is_time_reversal, 1561 is_mesh_symmetry=is_mesh_symmetry, 1562 with_eigenvectors=with_eigenvectors, 1563 with_group_velocities=with_group_velocities, 1564 is_gamma_center=is_gamma_center) 1565 self._mesh.run() 1566 1567 def set_mesh(self, 1568 mesh, 1569 shift=None, 1570 is_time_reversal=True, 1571 is_mesh_symmetry=True, 1572 is_eigenvectors=False, 1573 is_gamma_center=False, 1574 run_immediately=True): 1575 """Run or initialize phonon calculations on sampling mesh grids. 1576 1577 Parameters 1578 ---------- 1579 mesh: array_like 1580 Mesh numbers along a, b, c axes. 1581 dtype='intc' 1582 shape=(3,) 1583 shift: array_like, optional, default None (no shift) 1584 Mesh shifts along a*, b*, c* axes with respect to neighboring grid 1585 points from the original mesh (Monkhorst-Pack or Gamma center). 1586 0.5 gives half grid shift. Normally 0 or 0.5 is given. 1587 Otherwise q-points symmetry search is not performed. 1588 dtype='double' 1589 shape=(3, ) 1590 is_time_reversal: bool, optional, default True 1591 Time reversal symmetry is considered in symmetry search. By this, 1592 inversion symmetry is always included. 1593 is_mesh_symmetry: bool, optional, default True 1594 Wheather symmetry search is done or not. 1595 is_eigenvectors: bool, optional, default False 1596 Eigenvectors are stored by setting True. 1597 is_gamma_center: bool, default False 1598 Uniform mesh grids are generated centring at Gamma point but not 1599 the Monkhorst-Pack scheme. 1600 run_immediately: bool, default True 1601 With True, phonon calculations are performed immediately, which is 1602 usual usage. 1603 1604 """ 1605 warnings.warn("Phonopy.set_mesh is deprecated. " 1606 "Use Phonopy.run_mesh.", DeprecationWarning) 1607 1608 if self._group_velocity is None: 1609 with_group_velocities = False 1610 else: 1611 with_group_velocities = True 1612 if run_immediately: 1613 self.run_mesh(mesh, 1614 shift=shift, 1615 is_time_reversal=is_time_reversal, 1616 is_mesh_symmetry=is_mesh_symmetry, 1617 with_eigenvectors=is_eigenvectors, 1618 with_group_velocities=with_group_velocities, 1619 is_gamma_center=is_gamma_center) 1620 else: 1621 self.init_mesh(mesh, 1622 shift=shift, 1623 is_time_reversal=is_time_reversal, 1624 is_mesh_symmetry=is_mesh_symmetry, 1625 with_eigenvectors=is_eigenvectors, 1626 with_group_velocities=with_group_velocities, 1627 is_gamma_center=is_gamma_center) 1628 1629 def get_mesh_dict(self): 1630 """Return phonon properties calculated by mesh sampling. 1631 1632 Returns 1633 ------- 1634 dict 1635 keys: qpoints, weights, frequencies, eigenvectors, and 1636 group_velocities 1637 1638 Each value for the corresponding key is explained as below. 1639 1640 qpoints: ndarray 1641 q-points in reduced coordinates of reciprocal lattice 1642 dtype='double' 1643 shape=(ir-grid points, 3) 1644 weights: ndarray 1645 Geometric q-point weights. Its sum is the number of grid 1646 points. 1647 dtype='intc' 1648 shape=(ir-grid points,) 1649 frequencies: ndarray 1650 Phonon frequencies at ir-grid points. Imaginary frequenies are 1651 represented by negative real numbers. 1652 dtype='double' 1653 shape=(ir-grid points, bands) 1654 eigenvectors: ndarray 1655 Phonon eigenvectors at ir-grid points. See the data structure 1656 at np.linalg.eigh. 1657 shape=(ir-grid points, bands, bands) 1658 dtype=complex of "c%d" % (np.dtype('double').itemsize * 2) 1659 order='C' 1660 group_velocities: ndarray 1661 Phonon group velocities at ir-grid points. 1662 dtype='double' 1663 shape=(ir-grid points, bands, 3) 1664 1665 """ 1666 if self._mesh is None: 1667 msg = ("run_mesh has to be done.") 1668 raise RuntimeError(msg) 1669 1670 retdict = {'qpoints': self._mesh.qpoints, 1671 'weights': self._mesh.weights, 1672 'frequencies': self._mesh.frequencies, 1673 'eigenvectors': self._mesh.eigenvectors, 1674 'group_velocities': self._mesh.group_velocities} 1675 1676 return retdict 1677 1678 def get_mesh(self): 1679 """Return phonon properties calculated by mesh sampling.""" 1680 warnings.warn("Phonopy.get_mesh() is deprecated. " 1681 "Use Phonopy.get_mesh_dict().", 1682 DeprecationWarning) 1683 1684 if self._mesh is None: 1685 msg = ("run_mesh has to be done.") 1686 raise RuntimeError(msg) 1687 1688 mesh_dict = self.get_mesh_dict() 1689 1690 return (mesh_dict['qpoints'], 1691 mesh_dict['weights'], 1692 mesh_dict['frequencies'], 1693 mesh_dict['eigenvectors']) 1694 1695 def get_mesh_grid_info(self): 1696 """Return grid point information of mesh sampling.""" 1697 warnings.warn("Phonopy.get_mesh_grid_info() is deprecated. " 1698 "Use attributes of phonon.mesh instance.", 1699 DeprecationWarning) 1700 if self._mesh is None: 1701 msg = ("run_mesh has to be done.") 1702 raise RuntimeError(msg) 1703 1704 return (self._mesh.grid_address, 1705 self._mesh.ir_grid_points, 1706 self._mesh.grid_mapping_table) 1707 1708 def write_hdf5_mesh(self): 1709 """Write mesh calculation results in hdf5 format.""" 1710 self._mesh.write_hdf5() 1711 1712 def write_yaml_mesh(self): 1713 """Write mesh calculation results in yaml format.""" 1714 self._mesh.write_yaml() 1715 1716 # Sampling mesh: 1717 # Solving dynamical matrices at q-points one-by-one as an iterator 1718 def set_iter_mesh(self, 1719 mesh, 1720 shift=None, 1721 is_time_reversal=True, 1722 is_mesh_symmetry=True, 1723 is_eigenvectors=False, 1724 is_gamma_center=False): 1725 """Create an IterMesh instance. 1726 1727 See set_mesh method. 1728 1729 """ 1730 warnings.warn("Phonopy.set_iter_mesh() is deprecated. " 1731 "Use Phonopy.run_mesh() with use_iter_mesh=True.", 1732 DeprecationWarning) 1733 1734 self.run_mesh(mesh=mesh, 1735 shift=shift, 1736 is_time_reversal=is_time_reversal, 1737 is_mesh_symmetry=is_mesh_symmetry, 1738 with_eigenvectors=is_eigenvectors, 1739 is_gamma_center=is_gamma_center, 1740 use_iter_mesh=True) 1741 1742 # Plot band structure and DOS (PDOS) together 1743 def plot_band_structure_and_dos(self, pdos_indices=None): 1744 """Plot band structure and DOS.""" 1745 import matplotlib.pyplot as plt 1746 if self._band_structure.labels: 1747 from matplotlib import rc 1748 rc('text', usetex=True) 1749 1750 if self._band_structure.is_legacy_plot: 1751 import matplotlib.gridspec as gridspec 1752 # plt.figure(figsize=(10, 6)) 1753 gs = gridspec.GridSpec(1, 2, width_ratios=[3, 1]) 1754 ax2 = plt.subplot(gs[0, 1]) 1755 if pdos_indices is None: 1756 self._total_dos.plot(ax2, 1757 ylabel="", 1758 draw_grid=False, 1759 flip_xy=True) 1760 else: 1761 self._pdos.plot(ax2, 1762 indices=pdos_indices, 1763 ylabel="", 1764 draw_grid=False, 1765 flip_xy=True) 1766 ax2.set_xlim((0, None)) 1767 plt.setp(ax2.get_yticklabels(), visible=False) 1768 1769 ax1 = plt.subplot(gs[0, 0], sharey=ax2) 1770 self._band_structure.plot(ax1) 1771 1772 plt.subplots_adjust(wspace=0.03) 1773 plt.tight_layout() 1774 else: 1775 from mpl_toolkits.axes_grid1 import ImageGrid 1776 n = len([x for x in self._band_structure.path_connections 1777 if not x]) + 1 1778 fig = plt.figure() 1779 axs = ImageGrid(fig, 111, # similar to subplot(111) 1780 nrows_ncols=(1, n), 1781 axes_pad=0.11, 1782 label_mode="L") 1783 self._band_structure.plot(axs[:-1]) 1784 1785 if pdos_indices is None: 1786 self._total_dos.plot(axs[-1], 1787 xlabel="", 1788 ylabel="", 1789 draw_grid=False, 1790 flip_xy=True) 1791 else: 1792 self._pdos.plot(axs[-1], 1793 indices=pdos_indices, 1794 xlabel="", 1795 ylabel="", 1796 draw_grid=False, 1797 flip_xy=True) 1798 xlim = axs[-1].get_xlim() 1799 ylim = axs[-1].get_ylim() 1800 aspect = (xlim[1] - xlim[0]) / (ylim[1] - ylim[0]) * 3 1801 axs[-1].set_aspect(aspect) 1802 axs[-1].axhline(y=0, linestyle=':', linewidth=0.5, color='b') 1803 axs[-1].set_xlim((0, None)) 1804 1805 return plt 1806 1807 # Sampling at q-points 1808 def run_qpoints(self, 1809 q_points, 1810 with_eigenvectors=False, 1811 with_group_velocities=False, 1812 with_dynamical_matrices=False, 1813 nac_q_direction=None): 1814 """Run phonon calculation at specified q-points. 1815 1816 Parameters 1817 ---------- 1818 q_points: array_like or float, optional 1819 q-points in reduced coordinates. 1820 dtype='double', shape=(q-points, 3) 1821 with_eigenvectors: bool, optional 1822 Eigenvectors are stored by setting True. Default False. 1823 with_group_velocities : bool, optional 1824 Group velocities are calculated by setting True. Default is False. 1825 with_dynamical_matrices : bool, optional 1826 Calculated dynamical matrices are stored by setting True. 1827 Default is False. 1828 nac_q_direction : array_like 1829 q-point direction from Gamma-point in fractional coordinates of 1830 reciprocal basis vectors. Only the direction is used, i.e., 1831 (q_direction / |q_direction|) is computed and used. This parameter 1832 is activated only at q=(0, 0, 0). 1833 shape=(3,), dtype='double' 1834 1835 """ 1836 if self._dynamical_matrix is None: 1837 msg = ("Dynamical matrix has not yet built.") 1838 raise RuntimeError(msg) 1839 1840 if with_group_velocities: 1841 if self._group_velocity is None: 1842 self._set_group_velocity() 1843 group_velocity = self._group_velocity 1844 else: 1845 group_velocity = None 1846 1847 self._qpoints = QpointsPhonon( 1848 np.reshape(q_points, (-1, 3)), 1849 self._dynamical_matrix, 1850 nac_q_direction=nac_q_direction, 1851 with_eigenvectors=with_eigenvectors, 1852 group_velocity=group_velocity, 1853 with_dynamical_matrices=with_dynamical_matrices, 1854 factor=self._factor) 1855 1856 def set_qpoints_phonon(self, 1857 q_points, 1858 nac_q_direction=None, 1859 is_eigenvectors=False, 1860 write_dynamical_matrices=False): 1861 """Run phonon calculation at specified q-points.""" 1862 warnings.warn("Phonopy.set_qpoints_phonon() is deprecated. " 1863 "Use Phonopy.run_qpoints().", DeprecationWarning) 1864 if self._group_velocity is None: 1865 with_group_velocities = False 1866 else: 1867 with_group_velocities = True 1868 self.run_qpoints( 1869 q_points, 1870 with_eigenvectors=is_eigenvectors, 1871 with_group_velocities=with_group_velocities, 1872 with_dynamical_matrices=write_dynamical_matrices, 1873 nac_q_direction=nac_q_direction) 1874 1875 def get_qpoints_dict(self): 1876 """Return calculated phonon properties at q-points. 1877 1878 Returns 1879 ------- 1880 dict 1881 keys: frequencies, eigenvectors, and dynamical_matrices 1882 1883 frequencies : ndarray 1884 Phonon frequencies. Imaginary frequenies are represented by 1885 negative real numbers. 1886 shape=(qpoints, bands), dtype='double' 1887 eigenvectors : ndarray 1888 Phonon eigenvectors. None if eigenvectors are not stored. 1889 shape=(qpoints, bands, bands) 1890 dtype=complex of "c%d" % (np.dtype('double').itemsize * 2) 1891 order='C' 1892 group_velocities : ndarray 1893 Phonon group velocities. None if group velocities are not 1894 calculated. 1895 shape=(qpoints, bands, 3), dtype='double' 1896 dynamical_matrices : ndarray 1897 Dynamical matrices at q-points. 1898 shape=(qpoints, bands, bands), dtype='double' 1899 1900 """ 1901 if self._qpoints is None: 1902 msg = ("Phonopy.run_qpoints() has to be done.") 1903 raise RuntimeError(msg) 1904 1905 return {'frequencies': self._qpoints.frequencies, 1906 'eigenvectors': self._qpoints.eigenvectors, 1907 'group_velocities': self._qpoints.group_velocities, 1908 'dynamical_matrices': self._qpoints.dynamical_matrices} 1909 1910 def get_qpoints_phonon(self): 1911 """Return phonon properties calculated at q-points.""" 1912 warnings.warn("Phonopy.get_qpoints_phonon() is deprecated. " 1913 "Use Phonopy.run_get_qpoints_dict().", 1914 DeprecationWarning) 1915 qpt = self.get_qpoints_dict() 1916 return (qpt['frequencies'], qpt['eigenvectors']) 1917 1918 def write_hdf5_qpoints_phonon(self): 1919 """Write phonon properties calculated at q-points in hdf5 format.""" 1920 self._qpoints.write_hdf5() 1921 1922 def write_yaml_qpoints_phonon(self): 1923 """Write phonon properties calculated at q-points in yaml format.""" 1924 self._qpoints.write_yaml() 1925 1926 # DOS 1927 def run_total_dos(self, 1928 sigma=None, 1929 freq_min=None, 1930 freq_max=None, 1931 freq_pitch=None, 1932 use_tetrahedron_method=True): 1933 """Run total DOS calculation. 1934 1935 Parameters 1936 ---------- 1937 sigma : float, optional 1938 Smearing width for smearing method. Default is None 1939 freq_min, freq_max, freq_pitch : float, optional 1940 Minimum and maximum frequencies in which range DOS is computed 1941 with the specified interval (freq_pitch). 1942 Defaults are None and they are automatically determined. 1943 use_tetrahedron_method : float, optional 1944 Use tetrahedron method when this is True. When sigma is set, 1945 smearing method is used. 1946 1947 """ 1948 if self._mesh is None: 1949 msg = "run_mesh has to be done before DOS calculation." 1950 raise RuntimeError(msg) 1951 1952 total_dos = TotalDos(self._mesh, 1953 sigma=sigma, 1954 use_tetrahedron_method=use_tetrahedron_method) 1955 total_dos.set_draw_area(freq_min, freq_max, freq_pitch) 1956 total_dos.run() 1957 self._total_dos = total_dos 1958 1959 def set_total_DOS(self, 1960 sigma=None, 1961 freq_min=None, 1962 freq_max=None, 1963 freq_pitch=None, 1964 tetrahedron_method=False): 1965 """Run total DOS calculation.""" 1966 warnings.warn("Phonopy.set_total_DOS() is deprecated. " 1967 "Use Phonopy.run_total_DOS()", DeprecationWarning) 1968 1969 self.run_total_dos(sigma=sigma, 1970 freq_min=freq_min, 1971 freq_max=freq_max, 1972 freq_pitch=freq_pitch, 1973 use_tetrahedron_method=tetrahedron_method) 1974 1975 def auto_total_dos(self, 1976 mesh=100.0, 1977 is_time_reversal=True, 1978 is_mesh_symmetry=True, 1979 is_gamma_center=False, 1980 plot=False, 1981 write_dat=False, 1982 filename="total_dos.dat"): 1983 """Conveniently calculate and draw total DOS.""" 1984 self.run_mesh(mesh=mesh, 1985 is_time_reversal=is_time_reversal, 1986 is_mesh_symmetry=is_mesh_symmetry, 1987 is_gamma_center=is_gamma_center) 1988 self.run_total_dos() 1989 if write_dat: 1990 self.write_total_dos(filename=filename) 1991 if plot: 1992 return self.plot_total_dos() 1993 1994 def get_total_dos_dict(self): 1995 """Return total DOS. 1996 1997 Returns 1998 ------- 1999 A dictionary with keys of 'frequency_points' and 'total_dos'. 2000 Each value of corresponding key is as follows: 2001 2002 frequency_points: ndarray 2003 shape=(frequency_sampling_points, ), dtype='double' 2004 total_dos: 2005 shape=(frequency_sampling_points, ), dtype='double' 2006 2007 """ 2008 return {'frequency_points': self._total_dos.frequency_points, 2009 'total_dos': self._total_dos.dos} 2010 2011 def get_total_DOS(self): 2012 """Return total DOS. 2013 2014 Returns 2015 ------- 2016 A tuple with (frequency_points, total_dos). 2017 2018 frequency_points: ndarray 2019 shape=(frequency_sampling_points, ), dtype='double' 2020 total_dos: 2021 shape=(frequency_sampling_points, ), dtype='double' 2022 2023 """ 2024 warnings.warn("Phonopy.get_total_DOS() is deprecated. " 2025 "Use Phonopy.get_total_dos_dict().", DeprecationWarning) 2026 2027 dos = self.get_total_dos_dict() 2028 2029 return dos['frequency_points'], dos['total_dos'] 2030 2031 def set_Debye_frequency(self, freq_max_fit=None): 2032 """Calculate Debye frequency on top of total DOS.""" 2033 self._total_dos.set_Debye_frequency( 2034 len(self._primitive), freq_max_fit=freq_max_fit) 2035 2036 def get_Debye_frequency(self): 2037 """Return Debye frequency.""" 2038 return self._total_dos.get_Debye_frequency() 2039 2040 def plot_total_DOS(self): 2041 """Plot total DOS.""" 2042 warnings.warn("Phonopy.plot_total_DOS() is deprecated. " 2043 "Use Phonopy.plot_total_dos() (lowercase on DOS).", 2044 DeprecationWarning) 2045 return self.plot_total_dos() 2046 2047 def plot_total_dos(self): 2048 """Plot total DOS.""" 2049 if self._total_dos is None: 2050 msg = ("run_total_dos has to be done before plotting " 2051 "total DOS.") 2052 raise RuntimeError(msg) 2053 2054 import matplotlib.pyplot as plt 2055 2056 fig, ax = plt.subplots() 2057 self._total_dos.plot(ax, draw_grid=False) 2058 ax.set_ylim((0, None)) 2059 2060 return plt 2061 2062 def write_total_DOS(self, filename="total_dos.dat"): 2063 """Write total DOS to text file.""" 2064 warnings.warn("Phonopy.write_total_DOS() is deprecated. " 2065 "Use Phonopy.write_total_dos() (lowercase on DOS).", 2066 DeprecationWarning) 2067 self.write_total_dos(filename=filename) 2068 2069 def write_total_dos(self, filename="total_dos.dat"): 2070 """Write total DOS to text file.""" 2071 self._total_dos.write(filename=filename) 2072 2073 # PDOS 2074 def run_projected_dos(self, 2075 sigma=None, 2076 freq_min=None, 2077 freq_max=None, 2078 freq_pitch=None, 2079 use_tetrahedron_method=True, 2080 direction=None, 2081 xyz_projection=False): 2082 """Run projected DOS calculation. 2083 2084 Parameters 2085 ---------- 2086 sigma : float, optional 2087 Smearing width for smearing method. Default is None 2088 freq_min, freq_max, freq_pitch : float, optional 2089 Minimum and maximum frequencies in which range DOS is computed 2090 with the specified interval (freq_pitch). 2091 Defaults are None and they are automatically determined. 2092 use_tetrahedron_method : float, optional 2093 Use tetrahedron method when this is True. When sigma is set, 2094 smearing method is used. 2095 direction : array_like, optional 2096 Specific projection direction. This is specified three values 2097 along basis vectors or the primitive cell. Default is None, 2098 i.e., no projection. 2099 xyz_projection : bool, optional 2100 This determines whether projected along Cartesian directions or 2101 not. Default is False, i.e., no projection. 2102 2103 """ 2104 self._pdos = None 2105 2106 if self._mesh is None: 2107 msg = "run_mesh has to be done before PDOS calculation." 2108 raise RuntimeError(msg) 2109 2110 if not self._mesh.with_eigenvectors: 2111 msg = "run_mesh has to be called with with_eigenvectors=True." 2112 raise RuntimeError(msg) 2113 2114 if np.prod(self._mesh.mesh_numbers) != len(self._mesh.ir_grid_points): 2115 msg = "run_mesh has to be done with is_mesh_symmetry=False." 2116 raise RuntimeError(msg) 2117 2118 if direction is not None: 2119 direction_cart = np.dot(direction, self._primitive.cell) 2120 else: 2121 direction_cart = None 2122 self._pdos = PartialDos(self._mesh, 2123 sigma=sigma, 2124 use_tetrahedron_method=use_tetrahedron_method, 2125 direction=direction_cart, 2126 xyz_projection=xyz_projection) 2127 self._pdos.set_draw_area(freq_min, freq_max, freq_pitch) 2128 self._pdos.run() 2129 2130 def set_partial_DOS(self, 2131 sigma=None, 2132 freq_min=None, 2133 freq_max=None, 2134 freq_pitch=None, 2135 tetrahedron_method=False, 2136 direction=None, 2137 xyz_projection=False): 2138 """Run projected DOS calculation.""" 2139 warnings.warn("Phonopy.set_partial_DOS() is deprecated. " 2140 "Use Phonopy.run_projected_dos()", DeprecationWarning) 2141 2142 self.run_projected_dos(sigma=sigma, 2143 freq_min=freq_min, 2144 freq_max=freq_max, 2145 freq_pitch=freq_pitch, 2146 use_tetrahedron_method=tetrahedron_method, 2147 direction=direction, 2148 xyz_projection=xyz_projection) 2149 2150 def auto_projected_dos(self, 2151 mesh=100.0, 2152 is_time_reversal=True, 2153 is_gamma_center=False, 2154 plot=False, 2155 pdos_indices=None, 2156 legend=None, 2157 write_dat=False, 2158 filename="projected_dos.dat"): 2159 """Conveniently calculate and draw projected DOS. 2160 2161 Parameters 2162 ---------- 2163 See docstring of ``Phonopy.init_mesh`` for the parameters of ``mesh`` 2164 (default is 100.0), ``is_time_reversal`` (default is True), 2165 and ``is_gamma_center`` (default is False). 2166 See docstring of ``Phonopy.plot_projected_dos`` for the parameters 2167 ``pdos_indices`` and ``legend``. 2168 2169 plot : Bool, optional 2170 With setting True, PDOS is plotted using matplotlib and 2171 the matplotlib module (plt) is returned. To watch the result, 2172 usually ``show()`` has to be called. Default is False. 2173 write_dat : Bool 2174 With setting True, ``projected_dos.dat`` like file is written out. 2175 The file name can be specified with the ``filename`` parameter. 2176 Default is False. 2177 filename : str, optional 2178 File name used to write ``projected_dos.dat`` like file. Default 2179 is ``projected_dos.dat``. 2180 2181 """ 2182 self.run_mesh(mesh=mesh, 2183 is_time_reversal=is_time_reversal, 2184 is_mesh_symmetry=False, 2185 with_eigenvectors=True, 2186 is_gamma_center=is_gamma_center) 2187 self.run_projected_dos() 2188 if write_dat: 2189 self.write_projected_dos(filename=filename) 2190 if plot: 2191 return self.plot_projected_dos(pdos_indices=pdos_indices, 2192 legend=legend) 2193 2194 def get_projected_dos_dict(self): 2195 """Return projected DOS. 2196 2197 Projection is done to atoms and may be also done along directions 2198 depending on the parameters at run_projected_dos. 2199 2200 Returns 2201 ------- 2202 A dictionary with keys of 'frequency_points' and 'projected_dos'. 2203 Each value of corresponding key is as follows: 2204 2205 frequency_points: ndarray 2206 shape=(frequency_sampling_points, ), dtype='double' 2207 partial_dos: 2208 shape=(frequency_sampling_points, projections), dtype='double' 2209 2210 """ 2211 return {'frequency_points': self._pdos.frequency_points, 2212 'projected_dos': self._pdos.partial_dos} 2213 2214 def get_partial_DOS(self): 2215 """Return projected DOS. 2216 2217 Projection is done to atoms and may be also done along directions 2218 depending on the parameters at run_partial_dos. 2219 2220 Returns 2221 ------- 2222 A tuple with (frequency_points, partial_dos). 2223 2224 frequency_points: ndarray 2225 shape=(frequency_sampling_points, ), dtype='double' 2226 partial_dos: 2227 shape=(frequency_sampling_points, projections), dtype='double' 2228 2229 """ 2230 warnings.warn("Phonopy.get_partial_DOS() is deprecated. " 2231 "Use Phonopy.get_projected_dos_dict().", 2232 DeprecationWarning) 2233 2234 pdos = self.get_projected_dos_dict() 2235 2236 return pdos['frequency_points'], pdos['projected_dos'] 2237 2238 def plot_partial_DOS(self, pdos_indices=None, legend=None): 2239 """Plot projected DOS.""" 2240 warnings.warn("Phonopy.plot_partial_DOS() is deprecated. " 2241 "Use Phonopy.plot_projected_dos() (lowercase on DOS).", 2242 DeprecationWarning) 2243 2244 return self.plot_projected_dos(pdos_indices=pdos_indices, 2245 legend=legend) 2246 2247 def plot_projected_dos(self, pdos_indices=None, legend=None): 2248 """Plot projected DOS. 2249 2250 Parameters 2251 ---------- 2252 pdos_indices : list of list, optional 2253 Sets of indices of atoms whose projected DOS are summed over. 2254 The indices start with 0. An example is as follwos: 2255 pdos_indices=[[0, 1], [2, 3, 4, 5]] 2256 Default is None, which means 2257 pdos_indices=[[i] for i in range(natom)] 2258 legend : list of instances such as str or int, optional 2259 The str(instance) are shown in legend. 2260 It has to be len(pdos_indices)==len(legend). Default is None. 2261 When None, legend is not shown. 2262 2263 """ 2264 import matplotlib.pyplot as plt 2265 2266 fig, ax = plt.subplots() 2267 ax.xaxis.set_ticks_position('both') 2268 ax.yaxis.set_ticks_position('both') 2269 ax.xaxis.set_tick_params(which='both', direction='in') 2270 ax.yaxis.set_tick_params(which='both', direction='in') 2271 2272 self._pdos.plot(ax, 2273 indices=pdos_indices, 2274 legend=legend, 2275 draw_grid=False) 2276 2277 ax.set_ylim((0, None)) 2278 2279 return plt 2280 2281 def write_partial_DOS(self, filename="partial_dos.dat"): 2282 """Write projected DOS to text file.""" 2283 warnings.warn("Phonopy.write_partial_DOS() is deprecated. " 2284 "Use Phonopy.write_projected_dos() (lowercase on DOS).", 2285 DeprecationWarning) 2286 self.write_projected_dos(filename=filename) 2287 2288 def write_projected_dos(self, filename="projected_dos.dat"): 2289 """Write projected DOS to text file.""" 2290 self._pdos.write(filename=filename) 2291 2292 # Thermal property 2293 def run_thermal_properties(self, 2294 t_min=0, 2295 t_max=1000, 2296 t_step=10, 2297 temperatures=None, 2298 is_projection=False, 2299 band_indices=None, 2300 cutoff_frequency=None, 2301 pretend_real=False): 2302 """Run calculation of thermal properties at constant volume. 2303 2304 Parameters 2305 ---------- 2306 t_min, t_max, t_step : float, optional 2307 Minimum and maximum temperatures and the interval in this 2308 temperature range. Default values are 0, 1000, and 10. 2309 temperatures : array_like, optional 2310 Temperature points where thermal properties are calculated. 2311 When this is set, t_min, t_max, and t_step are ignored. 2312 2313 """ 2314 if self._mesh is None: 2315 msg = ("run_mesh has to be done before" 2316 "run_thermal_properties.") 2317 raise RuntimeError(msg) 2318 2319 tp = ThermalProperties(self._mesh, 2320 is_projection=is_projection, 2321 band_indices=band_indices, 2322 cutoff_frequency=cutoff_frequency, 2323 pretend_real=pretend_real) 2324 if temperatures is None: 2325 tp.set_temperature_range(t_step=t_step, 2326 t_max=t_max, 2327 t_min=t_min) 2328 else: 2329 tp.temperatures = temperatures 2330 tp.run() 2331 self._thermal_properties = tp 2332 2333 def set_thermal_properties(self, 2334 t_step=10, 2335 t_max=1000, 2336 t_min=0, 2337 temperatures=None, 2338 is_projection=False, 2339 band_indices=None, 2340 cutoff_frequency=None, 2341 pretend_real=False): 2342 """Run calculation of thermal properties at constant volume.""" 2343 warnings.warn("Phonopy.set_thermal_properties() is deprecated. " 2344 "Use Phonopy.run_thermal_properties()", 2345 DeprecationWarning) 2346 self.run_thermal_properties(t_step=t_step, 2347 t_max=t_max, 2348 t_min=t_min, 2349 temperatures=temperatures, 2350 is_projection=is_projection, 2351 band_indices=band_indices, 2352 cutoff_frequency=cutoff_frequency, 2353 pretend_real=pretend_real) 2354 2355 def get_thermal_properties_dict(self): 2356 """Return thermal properties. 2357 2358 Returns 2359 ------- 2360 A dictionary of thermal properties with keys of 'temperatures', 2361 'free_energy', 'entropy', and 'heat_capacity'. 2362 Each value of corresponding key is as follows: 2363 2364 temperatures: ndarray 2365 shape=(temperatures, ), dtype='double' 2366 free_energy : ndarray 2367 shape=(temperatures, ), dtype='double' 2368 entropy : ndarray 2369 shape=(temperatures, ), dtype='double' 2370 heat_capacity : ndarray 2371 shape=(temperatures, ), dtype='double' 2372 2373 """ 2374 keys = ('temperatures', 'free_energy', 'entropy', 'heat_capacity') 2375 return dict(zip(keys, self._thermal_properties.thermal_properties)) 2376 2377 def get_thermal_properties(self): 2378 """Return thermal properties. 2379 2380 Returns 2381 ------- 2382 (temperatures, free energy, entropy, heat capacity) 2383 2384 """ 2385 warnings.warn("Phonopy.get_thermal_properties() is deprecated. " 2386 "Use Phonopy.get_thermal_properties_dict().", 2387 DeprecationWarning) 2388 2389 tp = self.get_thermal_properties_dict() 2390 return (tp['temperatures'], 2391 tp['free_energy'], 2392 tp['entropy'], 2393 tp['heat_capacity']) 2394 2395 def plot_thermal_properties(self): 2396 """Plot thermal properties.""" 2397 import matplotlib.pyplot as plt 2398 plt.rcParams['pdf.fonttype'] = 42 2399 plt.rcParams['font.family'] = 'serif' 2400 2401 fig, ax = plt.subplots() 2402 ax.xaxis.set_ticks_position('both') 2403 ax.yaxis.set_ticks_position('both') 2404 ax.xaxis.set_tick_params(which='both', direction='in') 2405 ax.yaxis.set_tick_params(which='both', direction='in') 2406 2407 self._thermal_properties.plot(plt) 2408 2409 temps = self._thermal_properties.temperatures 2410 ax.set_xlim((0, temps[-1])) 2411 2412 return plt 2413 2414 def write_yaml_thermal_properties(self, 2415 filename='thermal_properties.yaml'): 2416 """Write thermal properties in yaml format.""" 2417 self._thermal_properties.write_yaml(filename=filename) 2418 2419 # Thermal displacement 2420 def run_thermal_displacements(self, 2421 t_min=0, 2422 t_max=1000, 2423 t_step=10, 2424 temperatures=None, 2425 direction=None, 2426 freq_min=None, 2427 freq_max=None): 2428 """Run thermal displacements calculation. 2429 2430 Parameters 2431 ---------- 2432 t_min, t_max, t_step : float, optional 2433 Minimum and maximum temperatures and the interval in this 2434 temperature range. Default valuues are 0, 1000, and 10. 2435 temperatures : array_like, optional 2436 Temperature points where thermal properties are calculated. 2437 When this is set, t_min, t_max, and t_step are ignored. 2438 direction : array_like, optional 2439 Projection direction in reduced coordinates. Default is None, 2440 i.e., no projection. 2441 dtype=float, shape=(3,) 2442 freq_min, freq_max : float, optional 2443 Phonon frequencies larger than freq_min and smaller than 2444 freq_max are included. Default is None, i.e., all phonons. 2445 2446 """ 2447 if self._dynamical_matrix is None: 2448 msg = ("Dynamical matrix has not yet built.") 2449 raise RuntimeError(msg) 2450 if self._mesh is None: 2451 msg = ("run_mesh has to be done.") 2452 raise RuntimeError(msg) 2453 mesh_nums = self._mesh.mesh_numbers 2454 ir_grid_points = self._mesh.ir_grid_points 2455 if not self._mesh.with_eigenvectors: 2456 msg = ("run_mesh has to be done with with_eigenvectors=True.") 2457 raise RuntimeError(msg) 2458 if np.prod(mesh_nums) != len(ir_grid_points): 2459 msg = ("run_mesh has to be done with is_mesh_symmetry=False.") 2460 raise RuntimeError(msg) 2461 2462 if direction is not None: 2463 projection_direction = np.dot(direction, self._primitive.cell) 2464 td = ThermalDisplacements( 2465 self._mesh, 2466 projection_direction=projection_direction, 2467 freq_min=freq_min, 2468 freq_max=freq_max) 2469 else: 2470 td = ThermalDisplacements(self._mesh, 2471 freq_min=freq_min, 2472 freq_max=freq_max) 2473 2474 if temperatures is None: 2475 td.set_temperature_range(t_min, t_max, t_step) 2476 else: 2477 td.set_temperatures(temperatures) 2478 td.run() 2479 2480 self._thermal_displacements = td 2481 2482 def set_thermal_displacements(self, 2483 t_step=10, 2484 t_max=1000, 2485 t_min=0, 2486 temperatures=None, 2487 direction=None, 2488 freq_min=None, 2489 freq_max=None): 2490 """Run thermal displacements calculation.""" 2491 warnings.warn("Phonopy.set_thermal_displacements() is deprecated. " 2492 "Use Phonopy.run_thermal_displacements()", 2493 DeprecationWarning) 2494 self.run_thermal_displacements(t_min=t_min, 2495 t_max=t_max, 2496 t_step=t_step, 2497 temperatures=temperatures, 2498 direction=direction, 2499 freq_min=freq_min, 2500 freq_max=freq_max) 2501 2502 def get_thermal_displacements_dict(self): 2503 """Return thermal displacements.""" 2504 if self._thermal_displacements is None: 2505 msg = ("run_thermal_displacements has to be done.") 2506 raise RuntimeError(msg) 2507 2508 td = self._thermal_displacements 2509 return {'temperatures': td.temperatures, 2510 'thermal_displacements': td.thermal_displacements} 2511 2512 def get_thermal_displacements(self): 2513 """Return thermal displacements.""" 2514 warnings.warn("Phonopy.get_thermal_displacements() is deprecated. " 2515 "Use Phonopy.get_thermal_displacements_dict()", 2516 DeprecationWarning) 2517 td = self.get_thermal_displacements_dict() 2518 return (td['temperatures'], td['thermal_displacements']) 2519 2520 def plot_thermal_displacements(self, is_legend=False): 2521 """Plot thermal displacements.""" 2522 import matplotlib.pyplot as plt 2523 2524 fig, ax = plt.subplots() 2525 ax.xaxis.set_ticks_position('both') 2526 ax.yaxis.set_ticks_position('both') 2527 ax.xaxis.set_tick_params(which='both', direction='in') 2528 ax.yaxis.set_tick_params(which='both', direction='in') 2529 2530 self._thermal_displacements.plot(plt, is_legend=is_legend) 2531 2532 temps, _ = self._thermal_displacements.get_thermal_displacements() 2533 ax.set_xlim((0, temps[-1])) 2534 2535 return plt 2536 2537 def write_yaml_thermal_displacements(self): 2538 """Write thermal displacements in yaml format.""" 2539 self._thermal_displacements.write_yaml() 2540 2541 # Thermal displacement matrix 2542 def run_thermal_displacement_matrices(self, 2543 t_min=0, 2544 t_max=1000, 2545 t_step=10, 2546 temperatures=None, 2547 freq_min=None, 2548 freq_max=None): 2549 """Run thermal displacement matrices calculation. 2550 2551 Parameters 2552 ---------- 2553 t_min, t_max, t_step : float, optional 2554 Minimum and maximum temperatures and the interval in this 2555 temperature range. Default valuues are 0, 1000, and 10. 2556 freq_min, freq_max : float, optional 2557 Phonon frequencies larger than freq_min and smaller than 2558 freq_max are included. Default is None, i.e., all phonons. 2559 temperatures : array_like, optional 2560 Temperature points where thermal properties are calculated. 2561 When this is set, t_min, t_max, and t_step are ignored. 2562 Default is None. 2563 2564 """ 2565 if self._dynamical_matrix is None: 2566 msg = ("Dynamical matrix has not yet built.") 2567 raise RuntimeError(msg) 2568 if self._mesh is None: 2569 msg = ("run_mesh has to be done.") 2570 raise RuntimeError(msg) 2571 mesh_nums = self._mesh.mesh_numbers 2572 ir_grid_points = self._mesh.ir_grid_points 2573 if not self._mesh.with_eigenvectors: 2574 msg = ("run_mesh has to be done with with_eigenvectors=True.") 2575 raise RuntimeError(msg) 2576 if np.prod(mesh_nums) != len(ir_grid_points): 2577 msg = ("run_mesh has to be done with is_mesh_symmetry=False.") 2578 raise RuntimeError(msg) 2579 2580 tdm = ThermalDisplacementMatrices( 2581 self._mesh, 2582 freq_min=freq_min, 2583 freq_max=freq_max, 2584 lattice=self._primitive.cell.T) 2585 2586 if temperatures is None: 2587 tdm.set_temperature_range(t_min, t_max, t_step) 2588 else: 2589 tdm.set_temperatures(temperatures) 2590 tdm.run() 2591 2592 self._thermal_displacement_matrices = tdm 2593 2594 def set_thermal_displacement_matrices(self, 2595 t_step=10, 2596 t_max=1000, 2597 t_min=0, 2598 freq_min=None, 2599 freq_max=None, 2600 t_cif=None): 2601 """Run thermal displacement matrices calculation.""" 2602 warnings.warn("Phonopy.set_thermal_displacement_matrices() is " 2603 "deprecated. Use Phonopy.run_thermal_displacements()", 2604 DeprecationWarning) 2605 if t_cif is None: 2606 temperatures = None 2607 else: 2608 temperatures = [t_cif, ] 2609 self.run_thermal_displacement_matrices(t_min=t_min, 2610 t_max=t_max, 2611 t_step=t_step, 2612 temperatures=temperatures, 2613 freq_min=freq_min, 2614 freq_max=freq_max) 2615 2616 def get_thermal_displacement_matrices_dict(self): 2617 """Return thermal displacement matrices.""" 2618 if self._thermal_displacement_matrices is None: 2619 msg = ("run_thermal_displacement_matrices has to be done.") 2620 raise RuntimeError(msg) 2621 2622 tdm = self._thermal_displacement_matrices 2623 return {'temperatures': tdm.temperatures, 2624 'thermal_displacement_matrices': 2625 tdm.thermal_displacement_matrices, 2626 'thermal_displacement_matrices_cif': 2627 tdm.thermal_displacement_matrices_cif} 2628 2629 def get_thermal_displacement_matrices(self): 2630 """Return thermal displacement matrices.""" 2631 warnings.warn("Phonopy.get_thermal_displacement_matrices() is " 2632 "deprecated. Use " 2633 "Phonopy.get_thermal_displacement_matrices_dict()", 2634 DeprecationWarning) 2635 tdm = self.get_thermal_displacement_matrices_dict() 2636 return (tdm['temperatures'], 2637 tdm['thermal_displacement_matrices']) 2638 2639 def write_yaml_thermal_displacement_matrices(self): 2640 """Write thermal displacement matrices in yaml format.""" 2641 self._thermal_displacement_matrices.write_yaml() 2642 2643 def write_thermal_displacement_matrix_to_cif(self, temperature_index): 2644 """Write thermal displacement matrices at a termperature in cif.""" 2645 self._thermal_displacement_matrices.write_cif(self._primitive, 2646 temperature_index) 2647 2648 def write_animation(self, 2649 q_point=None, 2650 anime_type='v_sim', 2651 band_index=None, 2652 amplitude=None, 2653 num_div=None, 2654 shift=None, 2655 filename=None): 2656 """Write atomic modulations in animation format. 2657 2658 Returns 2659 ------- 2660 str 2661 Output filename. 2662 2663 """ 2664 if self._dynamical_matrix is None: 2665 msg = ("Dynamical matrix has not yet built.") 2666 raise RuntimeError(msg) 2667 2668 if anime_type in ('arc', 'xyz', 'jmol', 'poscar'): 2669 if band_index is None or amplitude is None or num_div is None: 2670 msg = ("Parameters are not correctly set for animation.") 2671 raise RuntimeError(msg) 2672 2673 return write_animation(self._dynamical_matrix, 2674 q_point=q_point, 2675 anime_type=anime_type, 2676 band_index=band_index, 2677 amplitude=amplitude, 2678 num_div=num_div, 2679 shift=shift, 2680 factor=self._factor, 2681 filename=filename) 2682 2683 # Atomic modulation of normal mode 2684 def set_modulations(self, 2685 dimension, 2686 phonon_modes, 2687 delta_q=None, 2688 derivative_order=None, 2689 nac_q_direction=None): 2690 """Generate atomic displacements of phonon modes. 2691 2692 The design of this feature is not very satisfactory, and thus API. 2693 Therefore it should be reconsidered someday in the fugure. 2694 2695 Parameters 2696 ---------- 2697 dimension : array_like 2698 Supercell dimension with respect to the primitive cell. 2699 dtype='intc', shape=(3, ), (3, 3), (9, ) 2700 phonon_modes : list of phonon mode settings 2701 Each element of the outer list gives one phonon mode information: 2702 2703 [q-point, band index (int), amplitude (float), phase (float)] 2704 2705 In each list of the phonon mode information, the first element is 2706 a list that represents q-point in reduced coordinates. The second, 2707 third, and fourth elements show the band index starting with 0, 2708 amplitude, and phase factor, respectively. 2709 nac_q_direction : array_like 2710 q-point direction from Gamma-point in fractional coordinates of 2711 reciprocal basis vectors. Only the direction is used, i.e., 2712 (q_direction / |q_direction|) is computed and used. This parameter 2713 is activated only at q=(0, 0, 0). 2714 shape=(3,), dtype='double' 2715 2716 """ 2717 if self._dynamical_matrix is None: 2718 msg = ("Dynamical matrix has not yet built.") 2719 raise RuntimeError(msg) 2720 2721 self._modulation = Modulation(self._dynamical_matrix, 2722 dimension, 2723 phonon_modes, 2724 delta_q=delta_q, 2725 derivative_order=derivative_order, 2726 nac_q_direction=nac_q_direction, 2727 factor=self._factor) 2728 self._modulation.run() 2729 2730 def get_modulated_supercells(self): 2731 """Return cells with atom modulations. 2732 2733 list of PhonopyAtoms 2734 Modulated structures. 2735 2736 """ 2737 return self._modulation.get_modulated_supercells() 2738 2739 def get_modulations_and_supercell(self): 2740 """Return atomic modulations and perfect supercell. 2741 2742 (modulations, supercell) 2743 2744 modulations: Atomic modulations of supercell in Cartesian coordinates 2745 supercell: Supercell as an PhonopyAtoms instance. 2746 2747 """ 2748 return self._modulation.get_modulations_and_supercell() 2749 2750 def write_modulations(self): 2751 """Write modulated structures to MPOSCAR's.""" 2752 self._modulation.write() 2753 2754 def write_yaml_modulations(self): 2755 """Write atomic modulations in yaml format.""" 2756 self._modulation.write_yaml() 2757 2758 # Irreducible representation 2759 def set_irreps(self, 2760 q, 2761 is_little_cogroup=False, 2762 nac_q_direction=None, 2763 degeneracy_tolerance=1e-4): 2764 """Identify ir-reps of phonon modes. 2765 2766 The design of this API is not very satisfactory and is expceted 2767 to be redesined in the next major versions once the use case 2768 of the API for ir-reps feature becomes clearer. 2769 2770 nac_q_direction : array_like 2771 q-point direction from Gamma-point in fractional coordinates of 2772 reciprocal basis vectors. Only the direction is used, i.e., 2773 (q_direction / |q_direction|) is computed and used. This parameter 2774 is activated only at q=(0, 0, 0). 2775 shape=(3,), dtype='double' 2776 2777 """ 2778 if self._dynamical_matrix is None: 2779 msg = ("Dynamical matrix has not yet built.") 2780 raise RuntimeError(msg) 2781 2782 self._irreps = IrReps( 2783 self._dynamical_matrix, 2784 q, 2785 is_little_cogroup=is_little_cogroup, 2786 nac_q_direction=nac_q_direction, 2787 factor=self._factor, 2788 symprec=self._symprec, 2789 degeneracy_tolerance=degeneracy_tolerance, 2790 log_level=self._log_level) 2791 2792 return self._irreps.run() 2793 2794 def get_irreps(self): 2795 """Return Ir-reps.""" 2796 return self._irreps 2797 2798 def show_irreps(self, show_irreps=False): 2799 """Show Ir-reps.""" 2800 self._irreps.show(show_irreps=show_irreps) 2801 2802 def write_yaml_irreps(self, show_irreps=False): 2803 """Write Ir-reps in yaml format.""" 2804 self._irreps.write_yaml(show_irreps=show_irreps) 2805 2806 # Group velocity 2807 def set_group_velocity(self, q_length=None): 2808 """Prepare group velocity calculation.""" 2809 warnings.warn("Phonopy.set_group_velocity() is deprecated. " 2810 "No need to call this. gv_delta_q " 2811 "(q_length) is set at Phonopy.__init__().", 2812 DeprecationWarning) 2813 self._gv_delta_q = q_length 2814 self._set_group_velocity() 2815 2816 def get_group_velocity(self): 2817 """Return group velocities.""" 2818 warnings.warn("Phonopy.get_group_velocities_on_bands is deprecated. " 2819 "Use Phonopy.[mode].group_velocities attribute or " 2820 "Phonopy.get_[mode]_dict()[group_velocities], where " 2821 "[mode] is band_structure, mesh, or qpoints.", 2822 DeprecationWarning) 2823 return self._group_velocity.get_group_velocity() 2824 2825 def get_group_velocity_at_q(self, q_point): 2826 """Return group velocity at a q-point.""" 2827 if self._group_velocity is None: 2828 self._set_group_velocity() 2829 self._group_velocity.run([q_point]) 2830 return self._group_velocity.group_velocities[0] 2831 2832 def get_group_velocities_on_bands(self): 2833 """Return group velocities calculated on band structure.""" 2834 warnings.warn( 2835 "Phonopy.get_group_velocities_on_bands is deprecated. " 2836 "Use Phonopy.get_band_structure_dict()['group_velocities'].", 2837 DeprecationWarning) 2838 return self._band_structure.group_velocities 2839 2840 # Moment 2841 def run_moment(self, 2842 order=1, 2843 is_projection=False, 2844 freq_min=None, 2845 freq_max=None): 2846 """Run moment calculation.""" 2847 if self._mesh is None: 2848 msg = ("run_mesh has to be done before run_moment.") 2849 raise RuntimeError(msg) 2850 else: 2851 if is_projection: 2852 if self._mesh.eigenvectors is None: 2853 return RuntimeError( 2854 "run_mesh has to be done with with_eigenvectors=True.") 2855 self._moment = PhononMoment( 2856 self._mesh.frequencies, 2857 weights=self._mesh.weights, 2858 eigenvectors=self._mesh.eigenvectors) 2859 else: 2860 self._moment = PhononMoment( 2861 self._mesh.get_frequencies(), 2862 weights=self._mesh.get_weights()) 2863 if freq_min is not None or freq_max is not None: 2864 self._moment.set_frequency_range(freq_min=freq_min, 2865 freq_max=freq_max) 2866 self._moment.run(order=order) 2867 2868 def set_moment(self, 2869 order=1, 2870 is_projection=False, 2871 freq_min=None, 2872 freq_max=None): 2873 """Run moment calculation.""" 2874 warnings.warn("Phonopy.set_moment() is deprecated. " 2875 "Use Phonopy.run_moment().", DeprecationWarning) 2876 self.run_moment(order=order, 2877 is_projection=is_projection, 2878 freq_min=freq_min, 2879 freq_max=freq_max) 2880 2881 def get_moment(self): 2882 """Return moment.""" 2883 return self._moment.moment 2884 2885 def init_dynamic_structure_factor(self, 2886 Qpoints, 2887 T, 2888 atomic_form_factor_func=None, 2889 scattering_lengths=None, 2890 freq_min=None, 2891 freq_max=None): 2892 """Initialize dynamic structure factor calculation. 2893 2894 ******************************************************************* 2895 This is still an experimental feature. API can be changed without 2896 notification. 2897 ******************************************************************* 2898 2899 Need to call DynamicStructureFactor.run() to start calculation. 2900 2901 Parameters 2902 ---------- 2903 Qpoints: array_like 2904 Q-points in any Brillouin zone. 2905 dtype='double' 2906 shape=(qpoints, 3) 2907 T: float 2908 Temperature in K. 2909 atomic_form_factor_func: Function object 2910 Function that returns atomic form factor (``func`` below): 2911 2912 f_params = {'Na': [3.148690, 2.594987, 4.073989, 6.046925, 2913 0.767888, 0.070139, 0.995612, 14.1226457, 2914 0.968249, 0.217037, 0.045300], 2915 'Cl': [1.061802, 0.144727, 7.139886, 1.171795, 2916 6.524271, 19.467656, 2.355626, 60.320301, 2917 35.829404, 0.000436, -34.916604],b| 2918 2919 def get_func_AFF(f_params): 2920 def func(symbol, Q): 2921 return atomic_form_factor_WK1995(Q, f_params[symbol]) 2922 return func 2923 2924 scattering_lengths: dictionary 2925 Coherent scattering lengths averaged over isotopes and spins. 2926 Supposed for INS. For example, {'Na': 3.63, 'Cl': 9.5770}. 2927 freq_min, freq_min: float 2928 Minimum and maximum phonon frequencies to determine whether 2929 phonons are included in the calculation. 2930 2931 """ 2932 if self._mesh is None: 2933 msg = ("run_mesh has to be done before initializing dynamic" 2934 "structure factor.") 2935 raise RuntimeError(msg) 2936 2937 if not self._mesh.with_eigenvectors: 2938 msg = "run_mesh has to be called with with_eigenvectors=True." 2939 raise RuntimeError(msg) 2940 2941 if np.prod(self._mesh.mesh_numbers) != len(self._mesh.ir_grid_points): 2942 msg = "run_mesh has to be done with is_mesh_symmetry=False." 2943 raise RuntimeError(msg) 2944 2945 self._dynamic_structure_factor = DynamicStructureFactor( 2946 self._mesh, 2947 Qpoints, 2948 T, 2949 atomic_form_factor_func=atomic_form_factor_func, 2950 scattering_lengths=scattering_lengths, 2951 freq_min=freq_min, 2952 freq_max=freq_max) 2953 2954 def run_dynamic_structure_factor(self, 2955 Qpoints, 2956 T, 2957 atomic_form_factor_func=None, 2958 scattering_lengths=None, 2959 freq_min=None, 2960 freq_max=None): 2961 """Run dynamic structure factor calculation. 2962 2963 See the detail of parameters at 2964 Phonopy.init_dynamic_structure_factor(). 2965 2966 """ 2967 self.init_dynamic_structure_factor( 2968 Qpoints, 2969 T, 2970 atomic_form_factor_func=atomic_form_factor_func, 2971 scattering_lengths=scattering_lengths, 2972 freq_min=freq_min, 2973 freq_max=freq_max) 2974 self._dynamic_structure_factor.run() 2975 2976 def set_dynamic_structure_factor(self, 2977 Qpoints, 2978 T, 2979 atomic_form_factor_func=None, 2980 scattering_lengths=None, 2981 freq_min=None, 2982 freq_max=None, 2983 run_immediately=True): 2984 """Run dynamic structure factor calculation.""" 2985 warnings.warn("Phonopy.set_dynamic_structure_factor() is deprecated. " 2986 "Use Phonopy.run_dynamic_structure_factor()", 2987 DeprecationWarning) 2988 if run_immediately: 2989 self.run_dynamic_structure_factor( 2990 Qpoints, 2991 T, 2992 atomic_form_factor_func=atomic_form_factor_func, 2993 scattering_lengths=scattering_lengths, 2994 freq_min=freq_min, 2995 freq_max=freq_max) 2996 else: 2997 self.init_dynamic_structure_factor( 2998 Qpoints, 2999 T, 3000 atomic_form_factor_func=atomic_form_factor_func, 3001 scattering_lengths=scattering_lengths, 3002 freq_min=freq_min, 3003 freq_max=freq_max) 3004 3005 def get_dynamic_structure_factor(self): 3006 """Return dynamic structure factors.""" 3007 return (self._dynamic_structure_factor.qpoints, 3008 self._dynamic_structure_factor.dynamic_structure_factors) 3009 3010 def run_random_displacements(self, 3011 temperature, 3012 number_of_snapshots=1, 3013 random_seed=None, 3014 dist_func=None, 3015 cutoff_frequency=None): 3016 """Generate random displacements from phonon structure. 3017 3018 Some more details are written at generate_displacements. 3019 3020 temperature : float 3021 Temperature. 3022 number_of_snapshots : int 3023 Number of snapshots with random displacements created. 3024 random_seed : 32bit unsigned int 3025 Random seed. 3026 dist_func : str, None 3027 Harmonic oscillator distribution function either by 'quantum' 3028 or 'classical'. Default is None, corresponding to 'quantum'. 3029 cutoff_frequency : float 3030 Phonon frequency in THz below that phonons are ignored 3031 to generate random displacements. 3032 3033 """ 3034 self._random_displacements = RandomDisplacements( 3035 self._supercell, 3036 self._primitive, 3037 self._force_constants, 3038 dist_func=dist_func, 3039 cutoff_frequency=cutoff_frequency, 3040 factor=self._factor) 3041 self._random_displacements.run( 3042 temperature, 3043 number_of_snapshots=number_of_snapshots, 3044 random_seed=random_seed) 3045 3046 def save(self, 3047 filename="phonopy_params.yaml", 3048 settings=None, 3049 hdf5_settings=None): 3050 """Save phonopy parameters into file. 3051 3052 Parameters 3053 ---------- 3054 filename: str, optional 3055 File name. Default is "phonopy_params.yaml" 3056 settings: dict, optional 3057 It is described which parameters are written out. Only the 3058 settings expected to be updated from the following default 3059 settings are needed to be set in the dictionary. The 3060 possible parameters and their default settings are: 3061 {'force_sets': True, 3062 'displacements': True, 3063 'force_constants': False, 3064 'born_effective_charge': True, 3065 'dielectric_constant': True} 3066 This default settings are updated by {'force_constants': True} 3067 when dataset is None and force_constants is not None unless 3068 {'force_constants': False} is specified. 3069 hdf5_settings: dict, optional (To be implemented) 3070 Force constants and force_sets are stored in hdf5 file when 3071 they are activated in the dict. The dict has the following keys. 3072 The default filename is the filename of yaml file where '.yaml' is 3073 replaced by '.hdf5'. 3074 'filename' : str 3075 'force_constants': bool (default=False) 3076 'force_sets': bool (default=False) 3077 3078 """ 3079 if hdf5_settings is not None: 3080 msg = "hdf5_settings parameter has not yet been implemented." 3081 raise NotImplementedError(msg) 3082 3083 if settings is None: 3084 _settings = {} 3085 else: 3086 _settings = settings.copy() 3087 if _settings.get('force_constants') is False: 3088 pass 3089 elif (not forces_in_dataset(self.dataset) and 3090 self.force_constants is not None): 3091 _settings.update({'force_constants': True}) 3092 phpy_yaml = PhonopyYaml(settings=_settings) 3093 phpy_yaml.set_phonon_info(self) 3094 with open(filename, 'w') as w: 3095 w.write(str(phpy_yaml)) 3096 3097 ################### 3098 # private methods # 3099 ################### 3100 def _run_force_constants_from_forces(self, 3101 distributed_atom_list=None, 3102 fc_calculator=None, 3103 fc_calculator_options=None, 3104 decimals=None): 3105 if self._displacement_dataset is not None: 3106 if fc_calculator is not None: 3107 disps, forces = get_displacements_and_forces( 3108 self._displacement_dataset) 3109 self._force_constants = get_fc2( 3110 self._supercell, 3111 self._primitive, 3112 disps, 3113 forces, 3114 fc_calculator=fc_calculator, 3115 fc_calculator_options=fc_calculator_options, 3116 atom_list=distributed_atom_list, 3117 log_level=self._log_level, 3118 symprec=self._symprec) 3119 else: 3120 if 'displacements' in self._displacement_dataset: 3121 lines = [ 3122 "Type-II dataset for displacements and forces was " 3123 "given. fc_calculator", 3124 "(external force constants calculator) is required " 3125 "to produce force constants."] 3126 raise RuntimeError("\n".join(lines)) 3127 self._force_constants = get_phonopy_fc2( 3128 self._supercell, 3129 self._symmetry, 3130 self._displacement_dataset, 3131 atom_list=distributed_atom_list, 3132 decimals=decimals) 3133 3134 def _set_dynamical_matrix(self): 3135 self._dynamical_matrix = None 3136 3137 if self._is_symmetry and self._nac_params is not None: 3138 borns, epsilon = symmetrize_borns_and_epsilon( 3139 self._nac_params['born'], 3140 self._nac_params['dielectric'], 3141 self._primitive, 3142 symprec=self._symprec) 3143 nac_params = self._nac_params.copy() 3144 nac_params.update({'born': borns, 'dielectric': epsilon}) 3145 else: 3146 nac_params = self._nac_params 3147 3148 if self._supercell is None or self._primitive is None: 3149 raise RuntimeError("Supercell or primitive is not created.") 3150 if self._force_constants is None: 3151 raise RuntimeError("Force constants are not prepared.") 3152 if self._primitive.get_masses() is None: 3153 raise RuntimeError("Atomic masses are not correctly set.") 3154 self._dynamical_matrix = get_dynamical_matrix( 3155 self._force_constants, 3156 self._supercell, 3157 self._primitive, 3158 nac_params, 3159 self._frequency_scale_factor, 3160 self._dynamical_matrix_decimals, 3161 symprec=self._symprec, 3162 log_level=self._log_level) 3163 # DynamialMatrix instance transforms force constants in correct 3164 # type of numpy array. 3165 self._force_constants = self._dynamical_matrix.force_constants 3166 3167 if self._group_velocity is not None: 3168 self._set_group_velocity() 3169 3170 def _set_group_velocity(self): 3171 if self._dynamical_matrix is None: 3172 raise RuntimeError("Dynamical matrix has not yet built.") 3173 3174 if (self._dynamical_matrix.is_nac() and 3175 self._dynamical_matrix.nac_method == 'gonze' and 3176 self._gv_delta_q is None): # noqa: E129 3177 self._gv_delta_q = 1e-5 3178 if self._log_level: 3179 msg = "Group velocity calculation:\n" 3180 text = ("Analytical derivative of dynamical matrix is not " 3181 "implemented for NAC by Gonze et al. Instead " 3182 "numerical derivative of it is used with dq=1e-5 " 3183 "for group velocity calculation.") 3184 msg += textwrap.fill(text, 3185 initial_indent=" ", 3186 subsequent_indent=" ", 3187 width=70) 3188 print(msg) 3189 3190 self._group_velocity = GroupVelocity( 3191 self._dynamical_matrix, 3192 q_length=self._gv_delta_q, 3193 symmetry=self._primitive_symmetry, 3194 frequency_factor_to_THz=self._factor) 3195 3196 def _search_symmetry(self): 3197 self._symmetry = Symmetry(self._supercell, 3198 self._symprec, 3199 self._is_symmetry) 3200 3201 def _search_primitive_symmetry(self): 3202 self._primitive_symmetry = Symmetry(self._primitive, 3203 self._symprec, 3204 self._is_symmetry) 3205 3206 if (len(self._symmetry.get_pointgroup_operations()) != 3207 len(self._primitive_symmetry.get_pointgroup_operations())): # noqa: E129 E501 3208 print("Warning: Point group symmetries of supercell and primitive" 3209 "cell are different.") 3210 3211 def _build_supercell(self): 3212 self._supercell = get_supercell(self._unitcell, 3213 self._supercell_matrix, 3214 self._symprec) 3215 3216 def _build_supercells_with_displacements(self): 3217 all_positions = [] 3218 if 'first_atoms' in self._displacement_dataset: # type-1 3219 for disp in self._displacement_dataset['first_atoms']: 3220 positions = self._supercell.get_positions() 3221 positions[disp['number']] += disp['displacement'] 3222 all_positions.append(positions) 3223 elif 'displacements' in self._displacement_dataset: 3224 for disp in self._displacement_dataset['displacements']: 3225 all_positions.append(self._supercell.positions + disp) 3226 else: 3227 raise RuntimeError("displacement_dataset is not set.") 3228 3229 supercells = [] 3230 for positions in all_positions: 3231 supercells.append(PhonopyAtoms( 3232 numbers=self._supercell.numbers, 3233 masses=self._supercell.masses, 3234 magmoms=self._supercell.magnetic_moments, 3235 positions=positions, 3236 cell=self._supercell.cell, 3237 pbc=True)) 3238 self._supercells_with_displacements = supercells 3239 3240 def _build_primitive_cell(self): 3241 """Create primitive cell. 3242 3243 primitive_matrix: 3244 Relative axes of primitive cell to the input unit cell. 3245 Relative axes to the supercell is calculated by: 3246 supercell_matrix^-1 * primitive_matrix 3247 Therefore primitive cell lattice is finally calculated by: 3248 (supercell_lattice * (supercell_matrix)^-1 * primitive_matrix)^T 3249 3250 """ 3251 inv_supercell_matrix = np.linalg.inv(self._supercell_matrix) 3252 if self._primitive_matrix is None: 3253 trans_mat = inv_supercell_matrix 3254 else: 3255 trans_mat = np.dot(inv_supercell_matrix, self._primitive_matrix) 3256 3257 try: 3258 self._primitive = get_primitive( 3259 self._supercell, trans_mat, self._symprec, 3260 store_dense_svecs=self._store_dense_svecs) 3261 except ValueError: 3262 msg = ("Creating primitive cell is failed. " 3263 "PRIMITIVE_AXIS may be incorrectly specified.") 3264 raise RuntimeError(msg) 3265 3266 def _set_primitive_matrix(self, primitive_matrix): 3267 pmat = get_primitive_matrix(primitive_matrix, symprec=self._symprec) 3268 if isinstance(pmat, str) and pmat == 'auto': 3269 return guess_primitive_matrix(self._unitcell, 3270 symprec=self._symprec) 3271 else: 3272 return pmat 3273 3274 def _shape_supercell_matrix(self, smat): 3275 return shape_supercell_matrix(smat) 3276