1"""Phonopy command user interface.""" 2# Copyright (C) 2020 Atsushi Togo 3# All rights reserved. 4# 5# This file is part of phonopy. 6# 7# Redistribution and use in source and binary forms, with or without 8# modification, are permitted provided that the following conditions 9# are met: 10# 11# * Redistributions of source code must retain the above copyright 12# notice, this list of conditions and the following disclaimer. 13# 14# * Redistributions in binary form must reproduce the above copyright 15# notice, this list of conditions and the following disclaimer in 16# the documentation and/or other materials provided with the 17# distribution. 18# 19# * Neither the name of the phonopy project nor the names of its 20# contributors may be used to endorse or promote products derived 21# from this software without specific prior written permission. 22# 23# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 24# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 25# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS 26# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE 27# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, 28# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, 29# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 30# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 31# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 32# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN 33# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 34# POSSIBILITY OF SUCH DAMAGE. 35 36import sys 37import os 38import numpy as np 39from phonopy import Phonopy, __version__ 40from phonopy.harmonic.force_constants import compact_fc_to_full_fc 41from phonopy.structure.cells import print_cell 42from phonopy.structure.cells import isclose as cells_isclose 43from phonopy.structure.atoms import atom_data, symbol_map 44from phonopy.structure.dataset import forces_in_dataset 45from phonopy.interface.phonopy_yaml import PhonopyYaml 46from phonopy.interface.fc_calculator import fc_calculator_names 47from phonopy.interface.calculator import ( 48 write_supercells_with_displacements, 49 get_default_physical_units, get_default_displacement_distance) 50from phonopy.interface.vasp import create_FORCE_CONSTANTS 51from phonopy.phonon.band_structure import ( 52 get_band_qpoints, get_band_qpoints_by_seekpath) 53from phonopy.phonon.dos import get_pdos_indices 54from phonopy.file_IO import ( 55 parse_FORCE_CONSTANTS, parse_FORCE_SETS, 56 write_FORCE_CONSTANTS, write_force_constants_to_hdf5, 57 get_born_parameters, parse_QPOINTS, is_file_phonopy_yaml) 58from phonopy.cui.create_force_sets import create_FORCE_SETS 59from phonopy.cui.load_helper import ( 60 read_force_constants_from_hdf5, get_nac_params, 61 set_dataset_and_force_constants) 62from phonopy.cui.show_symmetry import check_symmetry 63from phonopy.cui.settings import PhonopyConfParser 64from phonopy.cui.collect_cell_info import collect_cell_info 65from phonopy.cui.phonopy_argparse import ( 66 get_parser, show_deprecated_option_warnings) 67 68 69# AA is created at http://www.network-science.de/ascii/. 70def print_phonopy(): 71 """Show phonopy logo.""" 72 print(r""" _ 73 _ __ | |__ ___ _ __ ___ _ __ _ _ 74 | '_ \| '_ \ / _ \| '_ \ / _ \ | '_ \| | | | 75 | |_) | | | | (_) | | | | (_) || |_) | |_| | 76 | .__/|_| |_|\___/|_| |_|\___(_) .__/ \__, | 77 |_| |_| |___/""") 78 79 80def print_version(version, package_name="phonopy"): 81 """Show phonopy version number.""" 82 try: 83 version_text = ('%s' % version).rjust(44) 84 import pkg_resources 85 dist = pkg_resources.get_distribution(package_name) 86 if dist.has_version(): 87 ver = dist.version.split('.') 88 if len(ver) > 3: 89 rev = ver[3] 90 version_text = ('%s-r%s' % (version, rev)).rjust(44) 91 except ImportError: 92 pass 93 except Exception as err: 94 if (err.__module__ == 'pkg_resources' and 95 err.__class__.__name__ == 'DistributionNotFound'): # noqa E129 96 pass 97 else: 98 raise 99 finally: 100 print(version_text) 101 print('') 102 103 104def print_end(): 105 """Show end banner.""" 106 print(r""" _ 107 ___ _ __ __| | 108 / _ \ '_ \ / _` | 109 | __/ | | | (_| | 110 \___|_| |_|\__,_| 111""") 112 113 114def print_error(): 115 """Show error banner.""" 116 print(r""" ___ _ __ _ __ ___ _ __ 117 / _ \ '__| '__/ _ \| '__| 118| __/ | | | | (_) | | 119 \___|_| |_| \___/|_| 120""") 121 122 123def print_attention(attention_text): 124 """Show attentinal information.""" 125 print("*" * 67) 126 print(attention_text) 127 print("*" * 67) 128 print('') 129 130 131def print_error_message(message): 132 """Show error message.""" 133 print('') 134 print(message) 135 136 137def file_exists(filename, log_level, is_any=False): 138 """Check existence of file.""" 139 if os.path.exists(filename): 140 return True 141 else: 142 if is_any: 143 return False 144 else: 145 error_text = "\"%s\" was not found." % filename 146 print_error_message(error_text) 147 if log_level > 0: 148 print_error() 149 sys.exit(1) 150 151 152def files_exist(filename_list, log_level, is_any=False): 153 """Check existence of files.""" 154 filenames = [] 155 for filename in filename_list: 156 if file_exists(filename, log_level, is_any=is_any): 157 filenames.append(filename) 158 159 if filenames: 160 return filenames 161 else: 162 if len(filenames) == 2: 163 all_filenames = "\"%s\" or \"%s\"" % tuple(filenames) 164 else: 165 all_filenames = ", ".join(["\"%s\"" % 166 fn for fn in filename_list[:-1]]) 167 all_filenames += " or \"%s\"" % filename_list[-1] 168 error_text = "Any of %s was not found." % all_filenames 169 print_error_message(error_text) 170 if log_level > 0: 171 print_error() 172 sys.exit(1) 173 174 175def finalize_phonopy(log_level, 176 settings, 177 confs, 178 phonon, 179 filename="phonopy.yaml"): 180 """Finalize phonopy.""" 181 units = get_default_physical_units(phonon.calculator) 182 183 if settings.save_params: 184 exists_fc_only = (not forces_in_dataset(phonon.dataset) and 185 phonon.force_constants is not None) 186 yaml_settings = { 187 'displacements': not exists_fc_only, 188 'force_sets': not exists_fc_only, 189 'force_constants': exists_fc_only, 190 'born_effective_charge': True, 191 'dielectric_constant': True 192 } 193 _filename = "phonopy_params.yaml" 194 else: 195 yaml_settings = { 196 'force_sets': settings.include_force_sets, 197 'force_constants': settings.include_force_constants, 198 'born_effective_charge': settings.include_nac_params, 199 'dielectric_constant': settings.include_nac_params, 200 'displacements': settings.include_displacements 201 } 202 _filename = filename 203 204 phpy_yaml = PhonopyYaml(configuration=confs, 205 physical_units=units, 206 settings=yaml_settings) 207 phpy_yaml.set_phonon_info(phonon) 208 with open(_filename, 'w') as w: 209 w.write(str(phpy_yaml)) 210 211 if log_level > 0: 212 print("") 213 if settings.save_params: 214 print("Summary of calculation and parameters were written " 215 "in \"%s\"." % _filename) 216 else: 217 print("Summary of calculation was written in \"%s\"." % _filename) 218 print_end() 219 sys.exit(0) 220 221 222def print_cells(phonon, unitcell_filename): 223 """Show cells.""" 224 supercell = phonon.get_supercell() 225 unitcell = phonon.get_unitcell() 226 primitive = phonon.get_primitive() 227 p2p_map = primitive.get_primitive_to_primitive_map() 228 mapping = np.array( 229 [p2p_map[x] for x in primitive.get_supercell_to_primitive_map()], 230 dtype='intc') 231 s_indep_atoms = phonon.get_symmetry().get_independent_atoms() 232 p_indep_atoms = mapping[s_indep_atoms] 233 u2s_map = supercell.get_unitcell_to_supercell_map() 234 print("-" * 30 + " primitive cell " + "-" * 30) 235 print_cell(primitive, stars=p_indep_atoms) 236 print("-" * 32 + " unit cell " + "-" * 33) # 32 + 11 + 33 = 76 237 u2u_map = supercell.get_unitcell_to_unitcell_map() 238 u_indep_atoms = [u2u_map[x] for x in s_indep_atoms] 239 print_cell(unitcell, mapping=mapping[u2s_map], stars=u_indep_atoms) 240 print("-" * 32 + " super cell " + "-" * 32) 241 print_cell(supercell, mapping=mapping, stars=s_indep_atoms) 242 print("-" * 76) 243 244 245def print_settings(settings, 246 phonon, 247 is_primitive_axes_auto, 248 unitcell_filename, 249 load_phonopy_yaml): 250 """Show setting info.""" 251 primitive_matrix = phonon.primitive_matrix 252 supercell_matrix = phonon.supercell_matrix 253 interface_mode = phonon.calculator 254 run_mode = settings.run_mode 255 if interface_mode: 256 print("Calculator interface: %s" % interface_mode) 257 if (settings.cell_filename is not None and 258 settings.cell_filename != unitcell_filename): # noqa E129 259 print("\"%s\" was not able to be used." 260 % settings.cell_filename) 261 print("Crystal structure was read from \"%s\"." % unitcell_filename) 262 physical_units = get_default_physical_units(interface_mode) 263 print("Unit of length: %s" % physical_units['length_unit']) 264 if is_band_auto(settings) and not is_primitive_axes_auto: 265 print("Automatic band structure mode forced automatic choice " 266 "of primitive axes.") 267 if run_mode == 'band': 268 if is_band_auto(settings): 269 print("Band structure mode (Auto)") 270 else: 271 print("Band structure mode") 272 if run_mode == 'mesh': 273 print("Mesh sampling mode") 274 if run_mode == 'band_mesh': 275 print("Band structure and mesh sampling mode") 276 if run_mode == 'anime': 277 print("Animation mode") 278 if run_mode == 'modulation': 279 print("Modulation mode") 280 if run_mode == 'irreps': 281 print("Ir-representation mode") 282 if run_mode == 'qpoints': 283 if settings.write_dynamical_matrices: 284 print("QPOINTS mode (dynamical matrices written out)") 285 else: 286 print("QPOINTS mode") 287 if ((run_mode == 'band' or run_mode == 'mesh' or run_mode == 'qpoints') and 288 settings.is_group_velocity): # noqa 129 289 gv_delta_q = settings.group_velocity_delta_q 290 if gv_delta_q is not None: 291 print(" With group velocity calculation (dq=%3.1e)" % gv_delta_q) 292 else: 293 print('') 294 if settings.create_displacements: 295 print("Displacements creation mode") 296 if not settings.is_plusminus_displacement == 'auto': 297 if settings.is_plusminus_displacement: 298 print(" Plus Minus displacement: full plus minus directions") 299 else: 300 print(" Plus Minus displacement: only one direction") 301 if not settings.is_diagonal_displacement: 302 print(" Diagonal displacement: off") 303 if settings.random_displacements: 304 print(" Number of supercells with random displacements: %d" 305 % settings.random_displacements) 306 if settings.random_seed is not None: 307 print(" Random seed: %d" % settings.random_seed) 308 309 print("Settings:") 310 if load_phonopy_yaml: 311 if not settings.is_nac: 312 print(" Non-analytical term correction (NAC): off") 313 else: 314 if settings.is_nac: 315 print(" Non-analytical term correction (NAC): on") 316 if settings.nac_q_direction is not None: 317 print(" NAC q-direction: %s" % settings.nac_q_direction) 318 if settings.fc_spg_symmetry: 319 print(" Enforce space group symmetry to force constants: on") 320 if load_phonopy_yaml: 321 if not settings.fc_symmetry: 322 print(" Force constants symmetrization: off") 323 else: 324 if settings.fc_symmetry: 325 print(" Force constants symmetrization: on") 326 if settings.lapack_solver: 327 print(" Use Lapack solver via Lapacke: on") 328 if settings.symmetry_tolerance is not None: 329 print(" Symmetry tolerance: %5.2e" 330 % settings.symmetry_tolerance) 331 if run_mode == 'mesh' or run_mode == 'band_mesh': 332 mesh = settings.mesh_numbers 333 if type(mesh) is float: 334 print(" Length for sampling mesh: %.1f" % mesh) 335 elif mesh is not None: 336 print(" Sampling mesh: %s" % np.array(mesh)) 337 if settings.is_thermal_properties: 338 cutoff_freq = settings.cutoff_frequency 339 if cutoff_freq is None: 340 pass 341 elif cutoff_freq < 0: 342 print(" - Thermal properties are calculatd with " 343 "absolute phonon frequnecy.") 344 else: 345 print(" - Phonon frequencies > %5.3f are used to calculate " 346 "thermal properties." % cutoff_freq) 347 elif (settings.is_thermal_displacements or 348 settings.is_thermal_displacement_matrices): 349 fmin = settings.min_frequency 350 fmax = settings.max_frequency 351 text = None 352 if (fmin is not None) and (fmax is not None): 353 text = " - Phonon frequency from %5.3f to %5.3f" % (fmin, 354 fmax) 355 text += " are used to calculate\n" 356 text += " thermal displacements." 357 elif (fmin is None) and (fmax is not None): 358 text = "Phonon frequency < %5.3f" % fmax 359 text = (" - %s are used to calculate thermal displacements." % 360 text) 361 elif (fmin is not None) and (fmax is None): 362 text = "Phonon frequency > %5.3f" % fmin 363 text = (" - %s are used to calculate thermal displacements." % 364 text) 365 if text: 366 print(text) 367 if (np.diag(np.diag(supercell_matrix)) - supercell_matrix).any(): 368 print(" Supercell matrix:") 369 for v in supercell_matrix: 370 print(" %s" % v) 371 else: 372 print(" Supercell: %s" % np.diag(supercell_matrix)) 373 if is_primitive_axes_auto or is_band_auto(settings): 374 print(" Primitive matrix (Auto):") 375 for v in primitive_matrix: 376 print(" %s" % v) 377 elif primitive_matrix is not None: 378 print(" Primitive matrix:") 379 for v in primitive_matrix: 380 print(" %s" % v) 381 382 383def write_displacements_files_then_exit(phonon, 384 settings, 385 confs, 386 optional_structure_info, 387 log_level): 388 """Write supercells with displacements and displacement dataset. 389 390 Note 391 ---- 392 From phonopy v1.15.0, displacement dataset is written into 393 phonopy_disp.yaml. 394 395 """ 396 cells_with_disps = phonon.supercells_with_displacements 397 additional_info = {'supercell_matrix': phonon.supercell_matrix} 398 write_supercells_with_displacements(phonon.calculator, 399 phonon.supercell, 400 cells_with_disps, 401 optional_structure_info, 402 additional_info=additional_info) 403 404 if log_level > 0: 405 print("\"phonopy_disp.yaml\" and supercells have been created.") 406 407 settings.set_include_displacements(True) 408 409 finalize_phonopy(log_level, 410 settings, 411 confs, 412 phonon, 413 filename="phonopy_disp.yaml") 414 415 416def create_FORCE_SETS_from_settings(settings, symprec, log_level): 417 """Create FORCE_SETS.""" 418 if settings.create_force_sets: 419 filenames = settings.create_force_sets 420 force_sets_zero_mode = False 421 elif settings.create_force_sets_zero: 422 filenames = settings.create_force_sets_zero 423 force_sets_zero_mode = True 424 else: 425 print_error_message("Something wrong for parsing arguments.") 426 sys.exit(0) 427 428 disp_filenames = files_exist( 429 ['phonopy_disp.yaml', 'disp.yaml'], log_level, is_any=True) 430 interface_mode = settings.calculator 431 432 if disp_filenames[0] == 'phonopy_disp.yaml': 433 try: 434 phpy_yaml = PhonopyYaml() 435 phpy_yaml.read('phonopy_disp.yaml') 436 if phpy_yaml.calculator is not None: 437 interface_mode = phpy_yaml.calculator # overwrite 438 disp_filename = 'phonopy_disp.yaml' 439 except KeyError: 440 file_exists('disp.yaml', log_level) 441 if log_level > 0: 442 print("\"phonopy_disp.yaml\" was found but wasn't used " 443 "because of the old-style format.") 444 disp_filename = 'disp.yaml' 445 else: 446 disp_filename = disp_filenames[0] 447 448 files_exist(filenames, log_level) 449 450 create_FORCE_SETS( 451 interface_mode, 452 filenames, 453 symmetry_tolerance=symprec, 454 force_sets_zero_mode=force_sets_zero_mode, 455 disp_filename=disp_filename, 456 log_level=log_level) 457 458 459def produce_force_constants(phonon, 460 settings, 461 phpy_yaml, 462 unitcell_filename, 463 log_level): 464 """Run force constants calculation.""" 465 num_satom = len(phonon.supercell) 466 p2s_map = phonon.primitive.p2s_map 467 is_full_fc = settings.fc_spg_symmetry or settings.is_full_fc 468 469 if settings.read_force_constants: 470 if settings.is_hdf5 or settings.readfc_format == 'hdf5': 471 try: 472 import h5py 473 except ImportError: 474 print_error_message("You need to install python-h5py.") 475 if log_level: 476 print_error() 477 sys.exit(1) 478 479 file_exists("force_constants.hdf5", log_level) 480 fc = read_force_constants_from_hdf5( 481 filename="force_constants.hdf5", 482 p2s_map=p2s_map, 483 calculator=phonon.calculator) 484 fc_filename = "force_constants.hdf5" 485 else: 486 file_exists("FORCE_CONSTANTS", log_level) 487 fc = parse_FORCE_CONSTANTS(filename="FORCE_CONSTANTS", 488 p2s_map=p2s_map) 489 fc_filename = "FORCE_CONSTANTS" 490 491 if log_level: 492 print("Force constants are read from \"%s\"." % fc_filename) 493 494 if fc.shape[1] != num_satom: 495 error_text = ("Number of atoms in supercell is not consistent " 496 "with the matrix shape of\nforce constants read " 497 "from ") 498 if settings.is_hdf5 or settings.readfc_format == 'hdf5': 499 error_text += "force_constants.hdf5.\n" 500 else: 501 error_text += "FORCE_CONSTANTS.\n" 502 error_text += ("Please carefully check DIM, FORCE_CONSTANTS, " 503 "and %s.") % unitcell_filename 504 print_error_message(error_text) 505 if log_level: 506 print_error() 507 sys.exit(1) 508 509 # Compact fc is expanded to full fc when full fc is required. 510 if is_full_fc and fc.shape[0] != fc.shape[1]: 511 fc = compact_fc_to_full_fc(phonon, fc, log_level=log_level) 512 513 phonon.force_constants = fc 514 else: 515 def read_force_sets_from_phonopy_yaml(phpy_yaml): 516 if (phpy_yaml.dataset is not None and 517 ('forces' in phpy_yaml.dataset or 518 ('first_atoms' in phpy_yaml.dataset and 519 'forces' in phpy_yaml.dataset['first_atoms'][0]))): 520 return phpy_yaml.dataset 521 else: 522 return None 523 524 force_sets = None 525 526 if phpy_yaml is not None: 527 force_sets = read_force_sets_from_phonopy_yaml(phpy_yaml) 528 if log_level: 529 if force_sets is None: 530 print("Force sets were not found in \"%s\"." 531 % unitcell_filename) 532 else: 533 print("Forces and displacements were read from \"%s\"." 534 % unitcell_filename) 535 536 if force_sets is None: 537 file_exists("FORCE_SETS", log_level) 538 force_sets = parse_FORCE_SETS(natom=num_satom) 539 if log_level: 540 print("Forces and displacements were read from \"%s\"." 541 % "FORCE_SETS") 542 543 if (log_level and 544 force_sets is not None and 545 'displacements' in force_sets): 546 print("%d snapshots were found." 547 % len(force_sets['displacements'])) 548 549 if 'natom' in force_sets: 550 natom = force_sets['natom'] 551 else: 552 natom = force_sets['forces'].shape[1] 553 if natom != num_satom: 554 error_text = "Number of atoms in supercell is not consistent with " 555 error_text += "the data in FORCE_SETS.\n" 556 error_text += ("Please carefully check DIM, FORCE_SETS," 557 " and %s") % unitcell_filename 558 print_error_message(error_text) 559 if log_level: 560 print_error() 561 sys.exit(1) 562 563 (fc_calculator, 564 fc_calculator_options) = get_fc_calculator_params(settings) 565 566 phonon.dataset = force_sets 567 if log_level: 568 if fc_calculator is not None: 569 print("Force constants calculation by %s starts." 570 % fc_calculator_names[fc_calculator]) 571 else: 572 print("Computing force constants...") 573 574 if is_full_fc: 575 # Need to calculate full force constant tensors 576 phonon.produce_force_constants( 577 fc_calculator=fc_calculator, 578 fc_calculator_options=fc_calculator_options) 579 else: 580 # Only force constants between atoms in primitive cell and 581 # supercell 582 phonon.produce_force_constants( 583 calculate_full_force_constants=False, 584 fc_calculator=fc_calculator, 585 fc_calculator_options=fc_calculator_options) 586 587 588def store_force_constants(phonon, 589 settings, 590 phpy_yaml, 591 unitcell_filename, 592 load_phonopy_yaml, 593 log_level): 594 """Calculate or read force constants.""" 595 physical_units = get_default_physical_units(phonon.calculator) 596 p2s_map = phonon.primitive.p2s_map 597 598 if load_phonopy_yaml: 599 is_full_fc = settings.fc_spg_symmetry or settings.is_full_fc 600 if phpy_yaml.force_constants is not None: 601 if log_level: 602 print("Force constants were read from \"%s\"." 603 % unitcell_filename) 604 605 fc = phpy_yaml.force_constants 606 607 if fc.shape[1] != len(phonon.supercell): 608 error_text = ("Number of atoms in supercell is not consistent " 609 "with the matrix shape of\nforce constants read " 610 "from %s." % unitcell_filename) 611 print_error_message(error_text) 612 if log_level: 613 print_error() 614 sys.exit(1) 615 616 # Compact fc is expanded to full fc when full fc is required. 617 if is_full_fc and fc.shape[0] != fc.shape[1]: 618 fc = compact_fc_to_full_fc(phonon, fc, log_level=log_level) 619 620 phonon.set_force_constants(fc, show_drift=(log_level > 0)) 621 622 if forces_in_dataset(phpy_yaml.dataset): 623 if log_level: 624 text = "Force sets were read from \"%s\"" % unitcell_filename 625 if phpy_yaml.force_constants is not None: 626 text += " but not to be used." 627 else: 628 text += "." 629 print(text) 630 631 if phpy_yaml.force_constants is None: 632 (fc_calculator, 633 fc_calculator_options) = get_fc_calculator_params(settings) 634 try: 635 set_dataset_and_force_constants( 636 phonon, 637 phpy_yaml.dataset, 638 None, 639 fc_calculator=fc_calculator, 640 fc_calculator_options=fc_calculator_options, 641 produce_fc=True, 642 symmetrize_fc=False, 643 is_compact_fc=(not is_full_fc), 644 log_level=log_level) 645 except RuntimeError as e: 646 print_error_message(str(e)) 647 if log_level: 648 print_error() 649 sys.exit(1) 650 else: 651 produce_force_constants(phonon, 652 settings, 653 phpy_yaml, 654 unitcell_filename, 655 log_level) 656 657 # Impose cutoff radius on force constants 658 cutoff_radius = settings.cutoff_radius 659 if cutoff_radius: 660 phonon.set_force_constants_zero_with_radius(cutoff_radius) 661 662 # Enforce space group symmetry to force constants 663 if settings.fc_spg_symmetry: 664 if log_level: 665 print("Force constants are symmetrized by space group operations.") 666 print("This may take some time...") 667 phonon.symmetrize_force_constants_by_space_group() 668 if not load_phonopy_yaml: 669 write_FORCE_CONSTANTS(phonon.get_force_constants(), 670 filename='FORCE_CONSTANTS_SPG') 671 if log_level: 672 print("Symmetrized force constants are written into " 673 "\"FORCE_CONSTANTS_SPG\".") 674 675 # Imporse translational invariance and index permulation symmetry to 676 # force constants 677 if settings.fc_symmetry: 678 phonon.symmetrize_force_constants() 679 680 # Write FORCE_CONSTANTS 681 if settings.write_force_constants: 682 if settings.is_hdf5 or settings.writefc_format == 'hdf5': 683 fc_unit = physical_units['force_constants_unit'] 684 write_force_constants_to_hdf5( 685 phonon.get_force_constants(), 686 p2s_map=p2s_map, 687 physical_unit=fc_unit, 688 compression=settings.hdf5_compression) 689 if log_level: 690 print("Force constants are written into " 691 "\"force_constants.hdf5\".") 692 else: 693 fc = phonon.force_constants 694 write_FORCE_CONSTANTS(fc, p2s_map=p2s_map) 695 if log_level: 696 print("Force constants are written into \"FORCE_CONSTANTS\".") 697 print(" Array shape of force constants is %s." 698 % str(fc.shape)) 699 if fc.shape[0] != fc.shape[1]: 700 print(" Use --full-fc option for full array of force " 701 "constants.") 702 703 if log_level: 704 print("") 705 706 707def store_nac_params(phonon, 708 settings, 709 phpy_yaml, 710 unitcell_filename, 711 log_level, 712 nac_factor=None, 713 load_phonopy_yaml=False): 714 """Calculate or read NAC params.""" 715 if nac_factor is None: 716 physical_units = get_default_physical_units(phonon.calculator) 717 _nac_factor = physical_units['nac_factor'] 718 else: 719 _nac_factor = nac_factor 720 721 if settings.is_nac: 722 def read_BORN(phonon): 723 with open("BORN") as f: 724 return get_born_parameters( 725 f, phonon.primitive, phonon.primitive_symmetry) 726 727 nac_params = None 728 729 if load_phonopy_yaml: 730 nac_params = get_nac_params(primitive=phonon.primitive, 731 nac_params=phpy_yaml.nac_params, 732 log_level=log_level) 733 if phpy_yaml.nac_params is not None and log_level: 734 print("NAC parameters were read from \"%s\"." 735 % unitcell_filename) 736 else: 737 if phpy_yaml: 738 nac_params = phpy_yaml.nac_params 739 if log_level: 740 if nac_params is None: 741 print("NAC parameters were not found in \"%s\"." 742 % unitcell_filename) 743 else: 744 print("NAC parameters were read from \"%s\"." 745 % unitcell_filename) 746 747 if nac_params is None and file_exists("BORN", log_level): 748 nac_params = read_BORN(phonon) 749 if nac_params is not None and log_level: 750 print("NAC parameters were read from \"%s\"." % "BORN") 751 752 if not nac_params: 753 error_text = "BORN file could not be read correctly." 754 print_error_message(error_text) 755 if log_level: 756 print_error() 757 sys.exit(1) 758 759 if nac_params is not None: 760 if nac_params['factor'] is None: 761 nac_params['factor'] = _nac_factor 762 if settings.nac_method is not None: 763 nac_params['method'] = settings.nac_method 764 phonon.nac_params = nac_params 765 if log_level: 766 dm = phonon.dynamical_matrix 767 if dm is not None: 768 if dm.is_nac(): 769 dm.show_nac_message() 770 print("") 771 772 if log_level > 1: 773 print("-" * 27 + " Dielectric constant " + "-" * 28) 774 for v in nac_params['dielectric']: 775 print(" %12.7f %12.7f %12.7f" % tuple(v)) 776 print("-" * 26 + " Born effective charges " + "-" * 26) 777 symbols = phonon.primitive.symbols 778 for i, (z, s) in enumerate(zip(nac_params['born'], symbols)): 779 for j, v in enumerate(z): 780 if j == 0: 781 text = "%5d %-2s" % (i + 1, s) 782 else: 783 text = " " 784 print("%s %12.7f %12.7f %12.7f" % ((text,) + tuple(v))) 785 print("-" * 76) 786 787 788def run(phonon, settings, plot_conf, log_level): 789 """Run phonon calculations.""" 790 interface_mode = phonon.calculator 791 physical_units = get_default_physical_units(interface_mode) 792 run_mode = settings.run_mode 793 794 # 795 # QPOINTS mode 796 # 797 if run_mode == 'qpoints': 798 if settings.read_qpoints: 799 q_points = parse_QPOINTS() 800 if log_level: 801 print("Frequencies at q-points given by QPOINTS:") 802 elif settings.qpoints: 803 q_points = settings.qpoints 804 if log_level: 805 print("Q-points that will be calculated at:") 806 for q in q_points: 807 print(" %s" % q) 808 else: 809 print_error_message("Q-points are not properly specified.") 810 if log_level: 811 print_error() 812 sys.exit(1) 813 phonon.run_qpoints( 814 q_points, 815 with_eigenvectors=settings.is_eigenvectors, 816 with_group_velocities=settings.is_group_velocity, 817 with_dynamical_matrices=settings.write_dynamical_matrices, 818 nac_q_direction=settings.nac_q_direction) 819 820 if settings.is_hdf5 or settings.qpoints_format == "hdf5": 821 phonon.write_hdf5_qpoints_phonon() 822 else: 823 phonon.write_yaml_qpoints_phonon() 824 825 # 826 # Band structure 827 # 828 if run_mode == 'band' or run_mode == 'band_mesh': 829 if settings.band_points is None: 830 npoints = 51 831 else: 832 npoints = settings.band_points 833 band_paths = settings.band_paths 834 835 if is_band_auto(settings): 836 print("SeeK-path is used to generate band paths.") 837 print(" About SeeK-path https://seekpath.readthedocs.io/ " 838 "(citation there-in)") 839 is_legacy_plot = False 840 bands, labels, path_connections = get_band_qpoints_by_seekpath( 841 phonon.primitive, npoints, 842 is_const_interval=settings.is_band_const_interval) 843 else: 844 is_legacy_plot = settings.is_legacy_plot 845 if settings.is_band_const_interval: 846 reclat = np.linalg.inv(phonon.primitive.cell) 847 bands = get_band_qpoints(band_paths, npoints=npoints, 848 rec_lattice=reclat) 849 else: 850 bands = get_band_qpoints(band_paths, npoints=npoints) 851 path_connections = [] 852 for paths in band_paths: 853 path_connections += [True, ] * (len(paths) - 2) 854 path_connections.append(False) 855 labels = settings.band_labels 856 857 if log_level: 858 print("Reciprocal space paths in reduced coordinates:") 859 for band in bands: 860 print("[%6.3f %6.3f %6.3f] --> [%6.3f %6.3f %6.3f]" % 861 (tuple(band[0]) + tuple(band[-1]))) 862 863 phonon.run_band_structure( 864 bands, 865 with_eigenvectors=settings.is_eigenvectors, 866 with_group_velocities=settings.is_group_velocity, 867 is_band_connection=settings.is_band_connection, 868 path_connections=path_connections, 869 labels=labels, 870 is_legacy_plot=is_legacy_plot) 871 if interface_mode is None: 872 comment = None 873 else: 874 comment = {'calculator': interface_mode, 875 'length_unit': physical_units['length_unit']} 876 877 if settings.is_hdf5 or settings.band_format == 'hdf5': 878 phonon.write_hdf5_band_structure(comment=comment) 879 else: 880 phonon.write_yaml_band_structure(comment=comment) 881 882 if plot_conf['plot_graph'] and run_mode != 'band_mesh': 883 plot = phonon.plot_band_structure() 884 if plot_conf['save_graph']: 885 plot.savefig('band.pdf') 886 else: 887 plot.show() 888 889 # 890 # mesh sampling 891 # 892 if run_mode == 'mesh' or run_mode == 'band_mesh': 893 mesh_numbers = settings.mesh_numbers 894 if mesh_numbers is None: 895 mesh_numbers = 50.0 896 mesh_shift = settings.mesh_shift 897 t_symmetry = settings.is_time_reversal_symmetry 898 q_symmetry = settings.is_mesh_symmetry 899 is_gamma_center = settings.is_gamma_center 900 901 if (settings.is_thermal_displacements or 902 settings.is_thermal_displacement_matrices): # noqa E129 903 if settings.cutoff_frequency is not None: 904 if log_level: 905 print_error_message( 906 "Use FMIN (--fmin) instead of CUTOFF_FREQUENCY " 907 "(--cutoff-freq).") 908 print_error() 909 sys.exit(1) 910 911 phonon.init_mesh(mesh=mesh_numbers, 912 shift=mesh_shift, 913 is_time_reversal=t_symmetry, 914 is_mesh_symmetry=q_symmetry, 915 with_eigenvectors=settings.is_eigenvectors, 916 is_gamma_center=is_gamma_center, 917 use_iter_mesh=True) 918 if log_level: 919 print("Mesh numbers: %s" % phonon.mesh_numbers) 920 else: 921 phonon.init_mesh( 922 mesh=mesh_numbers, 923 shift=mesh_shift, 924 is_time_reversal=t_symmetry, 925 is_mesh_symmetry=q_symmetry, 926 with_eigenvectors=settings.is_eigenvectors, 927 with_group_velocities=settings.is_group_velocity, 928 is_gamma_center=is_gamma_center) 929 if log_level: 930 print("Mesh numbers: %s" % phonon.mesh_numbers) 931 weights = phonon.mesh.weights 932 if q_symmetry: 933 print( 934 "Number of irreducible q-points on sampling mesh: " 935 "%d/%d" % (weights.shape[0], 936 np.prod(phonon.mesh_numbers))) 937 else: 938 print("Number of q-points on sampling mesh: %d" % 939 weights.shape[0]) 940 print("Calculating phonons on sampling mesh...") 941 942 phonon.mesh.run() 943 944 if settings.write_mesh: 945 if settings.is_hdf5 or settings.mesh_format == 'hdf5': 946 phonon.write_hdf5_mesh() 947 else: 948 phonon.write_yaml_mesh() 949 950 # 951 # Thermal property 952 # 953 if settings.is_thermal_properties: 954 955 if log_level: 956 if settings.is_projected_thermal_properties: 957 print("Calculating projected thermal properties...") 958 else: 959 print("Calculating thermal properties...") 960 t_step = settings.temperature_step 961 t_max = settings.max_temperature 962 t_min = settings.min_temperature 963 phonon.run_thermal_properties( 964 t_min=t_min, 965 t_max=t_max, 966 t_step=t_step, 967 is_projection=settings.is_projected_thermal_properties, 968 band_indices=settings.band_indices, 969 cutoff_frequency=settings.cutoff_frequency, 970 pretend_real=settings.pretend_real) 971 phonon.write_yaml_thermal_properties() 972 973 if log_level: 974 print("#%11s %15s%15s%15s%15s" % ('T [K]', 975 'F [kJ/mol]', 976 'S [J/K/mol]', 977 'C_v [J/K/mol]', 978 'E [kJ/mol]')) 979 tp = phonon.get_thermal_properties_dict() 980 temps = tp['temperatures'] 981 fe = tp['free_energy'] 982 entropy = tp['entropy'] 983 heat_capacity = tp['heat_capacity'] 984 for T, F, S, CV in zip(temps, fe, entropy, heat_capacity): 985 print(("%12.3f " + "%15.7f" * 4) % 986 (T, F, S, CV, F + T * S / 1000)) 987 988 if plot_conf['plot_graph']: 989 plot = phonon.plot_thermal_properties() 990 if plot_conf['save_graph']: 991 plot.savefig('thermal_properties.pdf') 992 else: 993 plot.show() 994 995 # 996 # Thermal displacements 997 # 998 elif (settings.is_thermal_displacements and 999 run_mode in ('mesh', 'band_mesh')): 1000 p_direction = settings.projection_direction 1001 if log_level and p_direction is not None: 1002 c_direction = np.dot(p_direction, phonon.primitive.cell) 1003 c_direction /= np.linalg.norm(c_direction) 1004 print("Projection direction: [%6.3f %6.3f %6.3f] " 1005 "(fractional)" % tuple(p_direction)) 1006 print(" [%6.3f %6.3f %6.3f] " 1007 "(Cartesian)" % tuple(c_direction)) 1008 if log_level: 1009 print("Calculating thermal displacements...") 1010 t_step = settings.temperature_step 1011 t_max = settings.max_temperature 1012 t_min = settings.min_temperature 1013 phonon.run_thermal_displacements( 1014 t_min=t_min, 1015 t_max=t_max, 1016 t_step=t_step, 1017 direction=p_direction, 1018 freq_min=settings.min_frequency, 1019 freq_max=settings.max_frequency) 1020 phonon.write_yaml_thermal_displacements() 1021 1022 if plot_conf['plot_graph']: 1023 plot = phonon.plot_thermal_displacements( 1024 plot_conf['with_legend']) 1025 if plot_conf['save_graph']: 1026 plot.savefig('thermal_displacement.pdf') 1027 else: 1028 plot.show() 1029 1030 # 1031 # Thermal displacement matrices 1032 # 1033 elif (settings.is_thermal_displacement_matrices and 1034 run_mode in ('mesh', 'band_mesh')): 1035 if log_level: 1036 print("Calculating thermal displacement matrices...") 1037 t_step = settings.temperature_step 1038 t_max = settings.max_temperature 1039 t_min = settings.min_temperature 1040 t_cif = settings.thermal_displacement_matrix_temperatue 1041 if t_cif is None: 1042 temperatures = None 1043 else: 1044 temperatures = [t_cif, ] 1045 phonon.run_thermal_displacement_matrices( 1046 t_step=t_step, 1047 t_max=t_max, 1048 t_min=t_min, 1049 temperatures=temperatures, 1050 freq_min=settings.min_frequency, 1051 freq_max=settings.max_frequency) 1052 phonon.write_yaml_thermal_displacement_matrices() 1053 if t_cif is not None: 1054 phonon.write_thermal_displacement_matrix_to_cif(0) 1055 1056 # 1057 # Projected DOS 1058 # 1059 elif (settings.pdos_indices is not None and 1060 run_mode in ('mesh', 'band_mesh')): 1061 p_direction = settings.projection_direction 1062 if (log_level and 1063 p_direction is not None and 1064 not settings.xyz_projection): # noqa E129 1065 c_direction = np.dot(p_direction, phonon.primitive.cell) 1066 c_direction /= np.linalg.norm(c_direction) 1067 print("Projection direction: [%6.3f %6.3f %6.3f] " 1068 "(fractional)" % tuple(p_direction)) 1069 print(" [%6.3f %6.3f %6.3f] " 1070 "(Cartesian)" % tuple(c_direction)) 1071 if log_level: 1072 print("Calculating projected DOS...") 1073 1074 phonon.run_projected_dos( 1075 sigma=settings.sigma, 1076 freq_min=settings.min_frequency, 1077 freq_max=settings.max_frequency, 1078 freq_pitch=settings.frequency_pitch, 1079 use_tetrahedron_method=settings.is_tetrahedron_method, 1080 direction=p_direction, 1081 xyz_projection=settings.xyz_projection) 1082 phonon.write_projected_dos() 1083 1084 if plot_conf['plot_graph']: 1085 pdos_indices = settings.pdos_indices 1086 if is_pdos_auto(settings): 1087 pdos_indices = get_pdos_indices( 1088 phonon.primitive_symmetry) 1089 legend = [phonon.primitive.symbols[x[0]] 1090 for x in pdos_indices] 1091 else: 1092 legend = [np.array(x) + 1 for x in pdos_indices] 1093 if run_mode != 'band_mesh': 1094 plot = phonon.plot_projected_dos( 1095 pdos_indices=pdos_indices, 1096 legend=legend) 1097 if plot_conf['save_graph']: 1098 plot.savefig('partial_dos.pdf') 1099 else: 1100 plot.show() 1101 1102 # 1103 # Total DOS 1104 # 1105 elif ((plot_conf['plot_graph'] or settings.is_dos_mode) and 1106 not is_pdos_auto(settings) and 1107 run_mode in ('mesh', 'band_mesh')): 1108 phonon.run_total_dos( 1109 sigma=settings.sigma, 1110 freq_min=settings.min_frequency, 1111 freq_max=settings.max_frequency, 1112 freq_pitch=settings.frequency_pitch, 1113 use_tetrahedron_method=settings.is_tetrahedron_method) 1114 1115 if log_level: 1116 print("Calculating DOS...") 1117 1118 if settings.fits_Debye_model: 1119 phonon.set_Debye_frequency() 1120 if log_level: 1121 debye_freq = phonon.get_Debye_frequency() 1122 print("Debye frequency: %10.5f" % debye_freq) 1123 phonon.write_total_dos() 1124 1125 if plot_conf['plot_graph'] and run_mode != 'band_mesh': 1126 plot = phonon.plot_total_dos() 1127 if plot_conf['save_graph']: 1128 plot.savefig('total_dos.pdf') 1129 else: 1130 plot.show() 1131 1132 # 1133 # Momemt 1134 # 1135 elif (settings.is_moment and run_mode in ('mesh', 'band_mesh')): 1136 freq_min = settings.min_frequency 1137 freq_max = settings.max_frequency 1138 if log_level: 1139 text = "Calculating moment of phonon states distribution" 1140 if freq_min is None and freq_max is None: 1141 text += "..." 1142 elif freq_min is None and freq_max is not None: 1143 text += "\nbelow frequency %.3f..." % freq_max 1144 elif freq_min is not None and freq_max is None: 1145 text += "\nabove frequency %.3f..." % freq_min 1146 elif freq_min is not None and freq_max is not None: 1147 text += ("\nbetween frequencies %.3f and %.3f..." % 1148 (freq_min, freq_max)) 1149 print(text) 1150 print('') 1151 print("Order| Total | Projected to atoms") 1152 if settings.moment_order is not None: 1153 phonon.run_moment(order=settings.moment_order, 1154 freq_min=freq_min, 1155 freq_max=freq_max, 1156 is_projection=False) 1157 total_moment = phonon.get_moment() 1158 phonon.run_moment(order=settings.moment_order, 1159 freq_min=freq_min, 1160 freq_max=freq_max, 1161 is_projection=True) 1162 text = " %3d |%10.5f | " % (settings.moment_order, 1163 total_moment) 1164 for m in phonon.get_moment(): 1165 text += "%10.5f " % m 1166 print(text) 1167 else: 1168 for i in range(3): 1169 phonon.run_moment(order=i, 1170 freq_min=freq_min, 1171 freq_max=freq_max, 1172 is_projection=False) 1173 total_moment = phonon.get_moment() 1174 phonon.run_moment(order=i, 1175 freq_min=freq_min, 1176 freq_max=freq_max, 1177 is_projection=True) 1178 text = " %3d |%10.5f | " % (i, total_moment) 1179 for m in phonon.get_moment(): 1180 text += "%10.5f " % m 1181 print(text) 1182 1183 # 1184 # Band structure and DOS are plotted simultaneously. 1185 # 1186 if (run_mode == 'band_mesh' and 1187 plot_conf['plot_graph'] and 1188 not settings.is_thermal_properties and 1189 not settings.is_thermal_displacements and 1190 not settings.is_thermal_displacement_matrices and 1191 not settings.is_thermal_distances): # noqa E129 1192 if settings.pdos_indices is not None: 1193 plot = phonon.plot_band_structure_and_dos( 1194 pdos_indices=pdos_indices) 1195 else: 1196 plot = phonon.plot_band_structure_and_dos() 1197 if plot_conf['save_graph']: 1198 plot.savefig('band_dos.pdf') 1199 else: 1200 plot.show() 1201 1202 # 1203 # Animation 1204 # 1205 elif run_mode == 'anime': 1206 anime_type = settings.anime_type 1207 if anime_type == "v_sim": 1208 q_point = settings.anime_qpoint 1209 amplitude = settings.anime_amplitude 1210 fname_out = phonon.write_animation(q_point=q_point, 1211 anime_type='v_sim', 1212 amplitude=amplitude) 1213 if log_level: 1214 print("Animation type: v_sim") 1215 print("q-point: [%6.3f %6.3f %6.3f]" % tuple(q_point)) 1216 else: 1217 amplitude = settings.anime_amplitude 1218 band_index = settings.anime_band_index 1219 division = settings.anime_division 1220 shift = settings.anime_shift 1221 fname_out = phonon.write_animation(anime_type=anime_type, 1222 band_index=band_index, 1223 amplitude=amplitude, 1224 num_div=division, 1225 shift=shift) 1226 if log_level: 1227 print("Animation type: %s" % anime_type) 1228 print("amplitude: %f" % amplitude) 1229 if anime_type != "jmol": 1230 print("band index: %d" % band_index) 1231 print("Number of images: %d" % division) 1232 if log_level: 1233 print("Animation was written in \"%s\". " % fname_out) 1234 1235 # 1236 # Modulation 1237 # 1238 elif run_mode == 'modulation': 1239 mod_setting = settings.modulation 1240 phonon_modes = mod_setting['modulations'] 1241 dimension = mod_setting['dimension'] 1242 if 'delta_q' in mod_setting: 1243 delta_q = mod_setting['delta_q'] 1244 else: 1245 delta_q = None 1246 derivative_order = mod_setting['order'] 1247 num_band = len(phonon.primitive) * 3 1248 1249 if log_level: 1250 if len(phonon_modes) == 1: 1251 print("Modulated structure with %s multiplicity was created." 1252 % dimension) 1253 else: 1254 print("Modulated structures with %s multiplicity were created." 1255 % dimension) 1256 1257 error_indices = [] 1258 for i, ph_mode in enumerate(phonon_modes): 1259 if ph_mode[1] < 0 or ph_mode[1] >= num_band: 1260 error_indices.append(i) 1261 if log_level: 1262 text = ("%d: q%s, band index=%d, amplitude=%f" 1263 % (i + 1, ph_mode[0], ph_mode[1] + 1, ph_mode[2])) 1264 if len(ph_mode) > 3: 1265 text += ", phase=%f" % ph_mode[3] 1266 print(text) 1267 1268 if error_indices: 1269 if log_level: 1270 lines = ["Band index of modulation %d is out of range." 1271 % (i + 1) for i in error_indices] 1272 print_error_message('\n'.join(lines)) 1273 print_error() 1274 sys.exit(1) 1275 1276 phonon.set_modulations(dimension, 1277 phonon_modes, 1278 delta_q=delta_q, 1279 derivative_order=derivative_order, 1280 nac_q_direction=settings.nac_q_direction) 1281 phonon.write_modulations() 1282 phonon.write_yaml_modulations() 1283 1284 # 1285 # Ir-representation 1286 # 1287 elif run_mode == 'irreps': 1288 if phonon.set_irreps( 1289 settings.irreps_q_point, 1290 is_little_cogroup=settings.is_little_cogroup, 1291 nac_q_direction=settings.nac_q_direction, 1292 degeneracy_tolerance=settings.irreps_tolerance): 1293 phonon.show_irreps(settings.show_irreps) 1294 phonon.write_yaml_irreps(settings.show_irreps) 1295 1296 1297def start_phonopy(**argparse_control): 1298 """Parse arguments and set some basic parameters.""" 1299 parser, deprecated = get_parser(**argparse_control) 1300 args = parser.parse_args() 1301 1302 # Set log level 1303 log_level = 1 1304 if args.verbose: 1305 log_level = 2 1306 if args.quiet or args.is_check_symmetry: 1307 log_level = 0 1308 if args.loglevel is not None: 1309 log_level = args.loglevel 1310 1311 if args.is_graph_save: 1312 import matplotlib 1313 matplotlib.use('Agg') 1314 1315 # Show phonopy logo 1316 if log_level: 1317 print_phonopy() 1318 print_version(__version__) 1319 if argparse_control.get('load_phonopy_yaml', False): 1320 print("Running in phonopy.load mode.") 1321 print("Python version %d.%d.%d" % sys.version_info[:3]) 1322 import spglib 1323 print("Spglib version %d.%d.%d" % spglib.get_version()) 1324 print("") 1325 1326 if deprecated: 1327 show_deprecated_option_warnings(deprecated) 1328 1329 return (args, log_level, 1330 {'plot_graph': args.is_graph_plot, 1331 'save_graph': args.is_graph_save, 1332 'with_legend': args.is_legend}) 1333 1334 1335def read_phonopy_settings(args, argparse_control, log_level): 1336 """Read phonopy settings.""" 1337 load_phonopy_yaml = argparse_control.get('load_phonopy_yaml', False) 1338 conf_filename = None 1339 1340 if load_phonopy_yaml: 1341 if args.conf_filename: 1342 conf_filename = args.conf_filename 1343 phonopy_conf_parser = PhonopyConfParser( 1344 filename=args.conf_filename, args=args, 1345 default_settings=argparse_control) 1346 else: 1347 phonopy_conf_parser = PhonopyConfParser( 1348 args=args, default_settings=argparse_control) 1349 if len(args.filename) > 0: 1350 file_exists(args.filename[0], log_level) 1351 cell_filename = args.filename[0] 1352 else: 1353 cell_filename = phonopy_conf_parser.settings.cell_filename 1354 else: 1355 if len(args.filename) > 0: 1356 if is_file_phonopy_yaml(args.filename[0]): 1357 phonopy_conf_parser = PhonopyConfParser(args=args) 1358 cell_filename = args.filename[0] 1359 else: 1360 conf_filename = args.filename[0] 1361 phonopy_conf_parser = PhonopyConfParser( 1362 filename=args.filename[0], args=args) 1363 cell_filename = phonopy_conf_parser.settings.cell_filename 1364 else: 1365 phonopy_conf_parser = PhonopyConfParser(args=args) 1366 cell_filename = phonopy_conf_parser.settings.cell_filename 1367 1368 confs = phonopy_conf_parser.confs.copy() 1369 settings = phonopy_conf_parser.settings 1370 1371 if log_level > 0 and conf_filename is not None: 1372 print("Phonopy configuration was read from \"%s\"." % 1373 conf_filename) 1374 1375 return settings, confs, cell_filename 1376 1377 1378def is_band_auto(settings): 1379 """Check whether automatic band paths setting or not.""" 1380 return (type(settings.band_paths) is str and settings.band_paths == 'auto') 1381 1382 1383def is_pdos_auto(settings): 1384 """Check whether automatic PDOS setting or not.""" 1385 return (settings.pdos_indices == 'auto') 1386 1387 1388def auto_primitive_axes(primitive_matrix): 1389 """Check whether automatic primitive matrix setting or not.""" 1390 return (type(primitive_matrix) is str and primitive_matrix == 'auto') 1391 1392 1393def get_fc_calculator_params(settings): 1394 """Return fc_calculator and fc_calculator_params from settings.""" 1395 fc_calculator = None 1396 if settings.fc_calculator is not None: 1397 if settings.fc_calculator.lower() in fc_calculator_names: 1398 fc_calculator = settings.fc_calculator.lower() 1399 1400 fc_calculator_options = None 1401 if settings.fc_calculator_options is not None: 1402 fc_calculator_options = settings.fc_calculator_options 1403 1404 return fc_calculator, fc_calculator_options 1405 1406 1407def get_cell_info(settings, cell_filename, symprec, log_level): 1408 """Return calculator interface and crystal structure information.""" 1409 cell_info = collect_cell_info( 1410 supercell_matrix=settings.supercell_matrix, 1411 primitive_matrix=settings.primitive_matrix, 1412 interface_mode=settings.calculator, 1413 cell_filename=cell_filename, 1414 chemical_symbols=settings.chemical_symbols, 1415 enforce_primitive_matrix_auto=is_band_auto(settings), 1416 symprec=symprec, 1417 return_dict=True) 1418 if type(cell_info) is str: 1419 print_error_message(cell_info) 1420 if log_level: 1421 print_error() 1422 sys.exit(1) 1423 1424 # (unitcell, supercell_matrix, primitive_matrix, 1425 # optional_structure_info, interface_mode, 1426 # phpy_yaml) = cell_info 1427 # unitcell_filename = optional_structure_info[0] 1428 1429 # Set magnetic moments 1430 magmoms = settings.magnetic_moments 1431 if magmoms is not None: 1432 unitcell = cell_info['unitcell'] 1433 if len(magmoms) == len(unitcell): 1434 unitcell.magnetic_moments = magmoms 1435 else: 1436 error_text = "Invalid MAGMOM setting" 1437 print_error_message(error_text) 1438 if log_level: 1439 print_error() 1440 sys.exit(1) 1441 1442 if auto_primitive_axes(cell_info['primitive_matrix']): 1443 error_text = ("'PRIMITIVE_AXES = auto', 'BAND = auto', or no DIM " 1444 "setting is not allowed with MAGMOM.") 1445 print_error_message(error_text) 1446 if log_level: 1447 print_error() 1448 sys.exit(1) 1449 1450 cell_info['magmoms'] = magmoms 1451 1452 return cell_info 1453 1454 1455def show_symmetry_info_then_exit(cell_info, symprec): 1456 """Show crystal structure information in yaml style.""" 1457 phonon = Phonopy(cell_info['unitcell'], 1458 np.eye(3, dtype=int), 1459 primitive_matrix=cell_info['primitive_matrix'], 1460 symprec=symprec, 1461 calculator=cell_info['interface_mode'], 1462 log_level=0) 1463 check_symmetry(phonon, cell_info['optional_structure_info']) 1464 sys.exit(0) 1465 1466 1467def check_supercell_in_yaml(cell_info, ph, log_level): 1468 """Check supercell size consistency.""" 1469 if (cell_info['phonopy_yaml'] is not None and 1470 cell_info['phonopy_yaml'].supercell is not None): 1471 if not cells_isclose(cell_info['phonopy_yaml'].supercell, 1472 ph.supercell): 1473 if log_level: 1474 print("Generated Supercell is inconsistent with that " 1475 "in \"%s\"." % 1476 cell_info['optional_structure_info'][0]) 1477 print_error() 1478 sys.exit(1) 1479 1480 1481def init_phonopy(settings, cell_info, symprec, log_level): 1482 """Prepare phonopy object.""" 1483 if settings.create_displacements and settings.temperatures is None: 1484 phonon = Phonopy(cell_info['unitcell'], 1485 cell_info['supercell_matrix'], 1486 primitive_matrix=cell_info['primitive_matrix'], 1487 symprec=symprec, 1488 is_symmetry=settings.is_symmetry, 1489 store_dense_svecs=settings.store_dense_svecs, 1490 calculator=cell_info['interface_mode'], 1491 log_level=log_level) 1492 else: # Read FORCE_SETS, FORCE_CONSTANTS, or force_constants.hdf5 1493 # Overwrite frequency unit conversion factor 1494 if settings.frequency_conversion_factor is not None: 1495 freq_factor = settings.frequency_conversion_factor 1496 else: 1497 physical_units = get_default_physical_units( 1498 cell_info['interface_mode']) 1499 freq_factor = physical_units['factor'] 1500 1501 phonon = Phonopy( 1502 cell_info['unitcell'], 1503 cell_info['supercell_matrix'], 1504 primitive_matrix=cell_info['primitive_matrix'], 1505 factor=freq_factor, 1506 frequency_scale_factor=settings.frequency_scale_factor, 1507 dynamical_matrix_decimals=settings.dm_decimals, 1508 force_constants_decimals=settings.fc_decimals, 1509 group_velocity_delta_q=settings.group_velocity_delta_q, 1510 symprec=symprec, 1511 is_symmetry=settings.is_symmetry, 1512 store_dense_svecs=settings.store_dense_svecs, 1513 calculator=cell_info['interface_mode'], 1514 use_lapack_solver=settings.lapack_solver, 1515 log_level=log_level) 1516 1517 check_supercell_in_yaml(cell_info, phonon, log_level) 1518 1519 # Set atomic masses of primitive cell 1520 if settings.masses is not None: 1521 phonon.masses = settings.masses 1522 1523 # Atomic species without mass case 1524 symbols_with_no_mass = [] 1525 if phonon.primitive.masses is None: 1526 for s in phonon.primitive.symbols: 1527 if (atom_data[symbol_map[s]][3] is None and 1528 s not in symbols_with_no_mass): 1529 symbols_with_no_mass.append(s) 1530 print_error_message( 1531 "Atomic mass of \'%s\' is not implemented in phonopy." % s) 1532 print_error_message( 1533 "MASS tag can be used to set atomic masses.") 1534 1535 if len(symbols_with_no_mass) > 0: 1536 if log_level: 1537 print_end() 1538 sys.exit(1) 1539 1540 return phonon 1541 1542 1543def main(**argparse_control): 1544 """Start phonopy.""" 1545 ############################################ 1546 # Parse phonopy conf and crystal structure # 1547 ############################################ 1548 load_phonopy_yaml = argparse_control.get('load_phonopy_yaml', False) 1549 1550 args, log_level, plot_conf = start_phonopy(**argparse_control) 1551 1552 settings, confs, cell_filename = read_phonopy_settings( 1553 args, argparse_control, log_level) 1554 1555 # phonopy --symmetry 1556 run_symmetry_info = args.is_check_symmetry 1557 1558 # ----------------------------------------------------------------------- 1559 # ----------------- 'args' should not be used below. -------------------- 1560 # ----------------------------------------------------------------------- 1561 1562 ########################################################### 1563 # Symmetry tolerance. Distance unit depends on interface. # 1564 ########################################################### 1565 if settings.symmetry_tolerance is None: 1566 symprec = 1e-5 1567 else: 1568 symprec = settings.symmetry_tolerance 1569 1570 ########################################## 1571 # Create FORCE_SETS (-f or --force_sets) # 1572 ########################################## 1573 if settings.create_force_sets or settings.create_force_sets_zero: 1574 create_FORCE_SETS_from_settings(settings, symprec, log_level) 1575 if log_level > 0: 1576 print_end() 1577 sys.exit(0) 1578 1579 #################################################################### 1580 # Create FORCE_CONSTANTS (--fc or --force_constants) only for VASP # 1581 #################################################################### 1582 if settings.create_force_constants: 1583 filename = settings.create_force_constants 1584 file_exists(filename, log_level) 1585 write_hdf5 = (settings.is_hdf5 or settings.writefc_format == 'hdf5') 1586 is_error = create_FORCE_CONSTANTS(filename, write_hdf5, log_level) 1587 if log_level: 1588 print_end() 1589 sys.exit(is_error) 1590 1591 ################################################################# 1592 # Parse crystal structure and optionally phonopy.yaml-like file # 1593 ################################################################# 1594 cell_info = get_cell_info(settings, cell_filename, symprec, log_level) 1595 unitcell_filename = cell_info['optional_structure_info'][0] 1596 1597 ########################################################### 1598 # Show crystal symmetry information and exit (--symmetry) # 1599 ########################################################### 1600 if run_symmetry_info: 1601 show_symmetry_info_then_exit(cell_info, symprec) 1602 1603 ###################### 1604 # Initialize phonopy # 1605 ###################### 1606 phonon = init_phonopy(settings, cell_info, symprec, log_level) 1607 1608 ################################################ 1609 # Show phonopy settings and crystal structures # 1610 ################################################ 1611 if log_level: 1612 print_settings(settings, 1613 phonon, 1614 auto_primitive_axes(cell_info['primitive_matrix']), 1615 unitcell_filename, 1616 load_phonopy_yaml) 1617 if cell_info['magmoms'] is None: 1618 print("Spacegroup: %s" % 1619 phonon.symmetry.get_international_table()) 1620 if log_level > 1: 1621 print_cells(phonon, unitcell_filename) 1622 else: 1623 print("Use -v option to watch primitive cell, unit cell, " 1624 "and supercell structures.") 1625 if log_level == 1: 1626 print("") 1627 1628 ######################################################### 1629 # Create constant amplitude displacements and then exit # 1630 ######################################################### 1631 if ((settings.create_displacements or settings.random_displacements) and 1632 settings.temperatures is None): 1633 if settings.displacement_distance is None: 1634 displacement_distance = get_default_displacement_distance( 1635 phonon.calculator) 1636 else: 1637 displacement_distance = settings.displacement_distance 1638 1639 phonon.generate_displacements( 1640 distance=displacement_distance, 1641 is_plusminus=settings.is_plusminus_displacement, 1642 is_diagonal=settings.is_diagonal_displacement, 1643 is_trigonal=settings.is_trigonal_displacement, 1644 number_of_snapshots=settings.random_displacements, 1645 random_seed=settings.random_seed) 1646 write_displacements_files_then_exit( 1647 phonon, 1648 settings, 1649 confs, 1650 cell_info['optional_structure_info'], 1651 log_level) 1652 1653 ################### 1654 # Force constants # 1655 ################### 1656 store_force_constants(phonon, 1657 settings, 1658 cell_info['phonopy_yaml'], 1659 unitcell_filename, 1660 load_phonopy_yaml, 1661 log_level) 1662 1663 ################################## 1664 # Non-analytical term correction # 1665 ################################## 1666 store_nac_params(phonon, 1667 settings, 1668 cell_info['phonopy_yaml'], 1669 unitcell_filename, 1670 log_level, 1671 load_phonopy_yaml=load_phonopy_yaml) 1672 1673 ################################################################### 1674 # Create random displacements at finite temperature and then exit # 1675 ################################################################### 1676 if settings.random_displacements and settings.temperatures is not None: 1677 phonon.generate_displacements( 1678 number_of_snapshots=settings.random_displacements, 1679 random_seed=settings.random_seed, 1680 temperature=settings.temperatures[0]) 1681 write_displacements_files_then_exit( 1682 phonon, 1683 settings, 1684 confs, 1685 cell_info['optional_structure_info'], 1686 log_level) 1687 1688 ####################### 1689 # Phonon calculations # 1690 ####################### 1691 if settings.run_mode not in ('band', 'mesh', 'band_mesh', 'anime', 1692 'modulation', 'irreps', 'qpoints'): 1693 print("-" * 76) 1694 print(" One of the following run modes may be specified for phonon " 1695 "calculations.") 1696 for mode in ['Mesh sampling (MESH, --mesh)', 1697 'Q-points (QPOINTS, --qpoints)', 1698 'Band structure (BAND, --band)', 1699 'Animation (ANIME, --anime)', 1700 'Modulation (MODULATION, --modulation)', 1701 'Characters of Irreps (IRREPS, --irreps)', 1702 'Create displacements (CREATE_DISPLACEMENTS, -d)']: 1703 print(" - %s" % mode) 1704 print("-" * 76) 1705 1706 run(phonon, settings, plot_conf, log_level) 1707 1708 ######################## 1709 # Phonopy finalization # 1710 ######################## 1711 finalize_phonopy(log_level, settings, confs, phonon) 1712