1# Copyright (C) 2020 Atsushi Togo 2# All rights reserved. 3# 4# This file is part of phono3py. 5# 6# Redistribution and use in source and binary forms, with or without 7# modification, are permitted provided that the following conditions 8# are met: 9# 10# * Redistributions of source code must retain the above copyright 11# notice, this list of conditions and the following disclaimer. 12# 13# * Redistributions in binary form must reproduce the above copyright 14# notice, this list of conditions and the following disclaimer in 15# the documentation and/or other materials provided with the 16# distribution. 17# 18# * Neither the name of the phonopy project nor the names of its 19# contributors may be used to endorse or promote products derived 20# from this software without specific prior written permission. 21# 22# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 23# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 24# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS 25# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE 26# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, 27# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, 28# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 29# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 30# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 31# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN 32# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 33# POSSIBILITY OF SUCH DAMAGE. 34 35import os 36import sys 37import numpy as np 38from phonopy.harmonic.force_constants import ( 39 show_drift_force_constants, 40 symmetrize_force_constants, 41 symmetrize_compact_force_constants) 42from phonopy.file_IO import get_dataset_type2, parse_FORCE_SETS 43from phonopy.cui.phonopy_script import print_error, file_exists 44from phonopy.interface.calculator import get_default_physical_units 45from phono3py.phonon3.fc3 import show_drift_fc3 46from phono3py.file_IO import ( 47 parse_disp_fc3_yaml, parse_disp_fc2_yaml, parse_FORCES_FC2, 48 parse_FORCES_FC3, read_fc3_from_hdf5, read_fc2_from_hdf5, 49 write_fc3_to_hdf5, write_fc2_to_hdf5, get_length_of_first_line) 50from phono3py.cui.show_log import show_phono3py_force_constants_settings 51from phono3py.phonon3.fc3 import ( 52 set_permutation_symmetry_fc3, set_translational_invariance_fc3) 53from phono3py.interface.phono3py_yaml import Phono3pyYaml 54 55 56def create_phono3py_force_constants(phono3py, 57 settings, 58 ph3py_yaml=None, 59 input_filename=None, 60 output_filename=None, 61 phono3py_yaml_filename=None, 62 log_level=1): 63 if settings.fc_calculator is None: 64 symmetrize_fc3r = (settings.is_symmetrize_fc3_r or 65 settings.fc_symmetry) 66 symmetrize_fc2 = (settings.is_symmetrize_fc2 or 67 settings.fc_symmetry) 68 else: # Rely on fc calculator the symmetrization of fc. 69 symmetrize_fc2 = False 70 symmetrize_fc3r = False 71 72 if log_level: 73 show_phono3py_force_constants_settings(settings) 74 75 ####### 76 # fc3 # 77 ####### 78 if (settings.is_joint_dos or 79 (settings.is_isotope and 80 not (settings.is_bterta or settings.is_lbte)) or 81 settings.read_gamma or 82 settings.read_pp or 83 (not settings.is_bterta and settings.write_phonon) or 84 settings.constant_averaged_pp_interaction is not None): 85 pass 86 else: 87 if settings.read_fc3: 88 _read_phono3py_fc3(phono3py, 89 symmetrize_fc3r, 90 input_filename, 91 log_level) 92 else: # fc3 from FORCES_FC3 or ph3py_yaml 93 _create_phono3py_fc3(phono3py, 94 ph3py_yaml, 95 symmetrize_fc3r, 96 symmetrize_fc2, 97 input_filename, 98 output_filename, 99 settings.is_compact_fc, 100 settings.cutoff_pair_distance, 101 settings.fc_calculator, 102 settings.fc_calculator_options, 103 log_level) 104 if output_filename is None: 105 filename = 'fc3.hdf5' 106 else: 107 filename = 'fc3.' + output_filename + '.hdf5' 108 if log_level: 109 print("Writing fc3 to \"%s\"." % filename) 110 write_fc3_to_hdf5(phono3py.fc3, 111 filename=filename, 112 p2s_map=phono3py.primitive.p2s_map, 113 compression=settings.hdf5_compression) 114 115 cutoff_distance = settings.cutoff_fc3_distance 116 if cutoff_distance is not None and cutoff_distance > 0: 117 if log_level: 118 print("Cutting-off fc3 by zero (cut-off distance: %f)" % 119 cutoff_distance) 120 phono3py.cutoff_fc3_by_zero(cutoff_distance) 121 122 if log_level: 123 show_drift_fc3(phono3py.fc3, primitive=phono3py.primitive) 124 125 ####### 126 # fc2 # 127 ####### 128 phonon_primitive = phono3py.phonon_primitive 129 p2s_map = phonon_primitive.p2s_map 130 if settings.read_fc2: 131 _read_phono3py_fc2(phono3py, 132 symmetrize_fc2, 133 input_filename, 134 log_level) 135 else: 136 if phono3py.phonon_supercell_matrix is None: 137 if (settings.fc_calculator == 'alm' and phono3py.fc2 is not None): 138 if log_level: 139 print("fc2 that was fit simultaneously with fc3 " 140 "by ALM is used.") 141 else: 142 _create_phono3py_fc2(phono3py, 143 ph3py_yaml, 144 symmetrize_fc2, 145 input_filename, 146 settings.is_compact_fc, 147 settings.fc_calculator, 148 settings.fc_calculator_options, 149 log_level) 150 else: 151 _create_phono3py_phonon_fc2(phono3py, 152 ph3py_yaml, 153 symmetrize_fc2, 154 input_filename, 155 settings.is_compact_fc, 156 settings.fc_calculator, 157 settings.fc_calculator_options, 158 log_level) 159 if output_filename is None: 160 filename = 'fc2.hdf5' 161 else: 162 filename = 'fc2.' + output_filename + '.hdf5' 163 if log_level: 164 print("Writing fc2 to \"%s\"." % filename) 165 write_fc2_to_hdf5(phono3py.fc2, 166 filename=filename, 167 p2s_map=p2s_map, 168 physical_unit='eV/angstrom^2', 169 compression=settings.hdf5_compression) 170 171 if log_level: 172 show_drift_force_constants(phono3py.fc2, 173 primitive=phonon_primitive, 174 name='fc2') 175 176 177def parse_forces(phono3py, 178 ph3py_yaml=None, 179 cutoff_pair_distance=None, 180 force_filename="FORCES_FC3", 181 disp_filename=None, 182 fc_type=None, 183 log_level=0): 184 filename_read_from = None 185 186 if fc_type == 'phonon_fc2': 187 natom = len(phono3py.phonon_supercell) 188 else: 189 natom = len(phono3py.supercell) 190 191 # Get dataset from ph3py_yaml. dataset can be None. 192 dataset = _extract_datast_from_ph3py_yaml(ph3py_yaml, fc_type) 193 if dataset: 194 filename_read_from = ph3py_yaml.yaml_filename 195 196 # Try to read FORCES_FC* if type-2 and return dataset. 197 # None is returned unless type-2. 198 # can emit FileNotFoundError. 199 if (dataset is None or 200 dataset is not None and not forces_in_dataset(dataset)): 201 _dataset = _get_type2_dataset(natom, 202 phono3py.calculator, 203 filename=force_filename, 204 log_level=log_level) 205 # Do not overwrite dataset when _dataset is None. 206 if _dataset: 207 filename_read_from = force_filename 208 dataset = _dataset 209 210 if dataset is None: 211 # Displacement dataset is obtained from disp_filename. 212 # can emit FileNotFoundError. 213 dataset = _read_disp_fc_yaml(disp_filename, fc_type) 214 filename_read_from = disp_filename 215 216 if 'natom' in dataset and dataset['natom'] != natom: 217 msg = ("Number of atoms in supercell is not consistent with " 218 "\"%s\"." % filename_read_from) 219 raise RuntimeError(msg) 220 221 if log_level and filename_read_from is not None: 222 print("Displacement dataset for %s was read from \"%s\"." 223 % (fc_type, filename_read_from)) 224 225 if cutoff_pair_distance: 226 if ('cutoff_distance' not in dataset or 227 ('cutoff_distance' in dataset and 228 cutoff_pair_distance < dataset['cutoff_distance'])): 229 dataset['cutoff_distance'] = cutoff_pair_distance 230 if log_level: 231 print("Cutoff-pair-distance: %f" % cutoff_pair_distance) 232 233 # Type-1 FORCES_FC*. 234 # dataset comes either from disp_fc*.yaml or phono3py*.yaml. 235 if not forces_in_dataset(dataset): 236 if fc_type == 'phonon_fc2': 237 parse_FORCES_FC2(dataset, filename=force_filename) 238 else: 239 parse_FORCES_FC3(dataset, filename=force_filename) 240 241 if log_level: 242 print("Sets of supercell forces were read from \"%s\"." 243 % force_filename) 244 sys.stdout.flush() 245 246 _convert_unit_in_dataset(dataset, phono3py.calculator) 247 248 return dataset 249 250 251def forces_in_dataset(dataset): 252 return ('forces' in dataset or 253 ('first_atoms' in dataset and 254 'forces' in dataset['first_atoms'][0])) 255 256 257def _read_disp_fc_yaml(disp_filename, fc_type): 258 if fc_type == 'phonon_fc2': 259 dataset = parse_disp_fc2_yaml(filename=disp_filename) 260 else: 261 dataset = parse_disp_fc3_yaml(filename=disp_filename) 262 263 return dataset 264 265 266def _read_phono3py_fc3(phono3py, 267 symmetrize_fc3r, 268 input_filename, 269 log_level): 270 if input_filename is None: 271 filename = 'fc3.hdf5' 272 else: 273 filename = 'fc3.' + input_filename + '.hdf5' 274 file_exists(filename, log_level) 275 if log_level: 276 print("Reading fc3 from \"%s\"." % filename) 277 278 p2s_map = phono3py.primitive.p2s_map 279 try: 280 fc3 = read_fc3_from_hdf5(filename=filename, p2s_map=p2s_map) 281 except RuntimeError: 282 import traceback 283 traceback.print_exc() 284 if log_level: 285 print_error() 286 sys.exit(1) 287 num_atom = phono3py.supercell.get_number_of_atoms() 288 if fc3.shape[1] != num_atom: 289 print("Matrix shape of fc3 doesn't agree with supercell size.") 290 if log_level: 291 print_error() 292 sys.exit(1) 293 294 if symmetrize_fc3r: 295 set_translational_invariance_fc3(fc3) 296 set_permutation_symmetry_fc3(fc3) 297 298 phono3py.fc3 = fc3 299 300 301def _read_phono3py_fc2(phono3py, 302 symmetrize_fc2, 303 input_filename, 304 log_level): 305 if input_filename is None: 306 filename = 'fc2.hdf5' 307 else: 308 filename = 'fc2.' + input_filename + '.hdf5' 309 file_exists(filename, log_level) 310 if log_level: 311 print("Reading fc2 from \"%s\"." % filename) 312 313 num_atom = phono3py.phonon_supercell.get_number_of_atoms() 314 p2s_map = phono3py.phonon_primitive.p2s_map 315 try: 316 phonon_fc2 = read_fc2_from_hdf5(filename=filename, p2s_map=p2s_map) 317 except RuntimeError: 318 import traceback 319 traceback.print_exc() 320 if log_level: 321 print_error() 322 sys.exit(1) 323 324 if phonon_fc2.shape[1] != num_atom: 325 print("Matrix shape of fc2 doesn't agree with supercell size.") 326 if log_level: 327 print_error() 328 sys.exit(1) 329 330 if symmetrize_fc2: 331 if phonon_fc2.shape[0] == phonon_fc2.shape[1]: 332 symmetrize_force_constants(phonon_fc2) 333 else: 334 symmetrize_compact_force_constants(phonon_fc2, 335 phono3py.phonon_primitive) 336 337 phono3py.fc2 = phonon_fc2 338 339 340def _get_type2_dataset(natom, calculator, filename="FORCES_FC3", log_level=0): 341 if not os.path.isfile(filename): 342 return None 343 344 with open(filename, 'r') as f: 345 len_first_line = get_length_of_first_line(f) 346 if len_first_line == 6: 347 dataset = get_dataset_type2(f, natom) 348 if log_level: 349 print("%d snapshots were found in %s." 350 % (len(dataset['displacements']), "FORCES_FC3")) 351 else: 352 dataset = None 353 return dataset 354 355 356def _create_phono3py_fc3(phono3py, 357 ph3py_yaml, 358 symmetrize_fc3r, 359 symmetrize_fc2, 360 input_filename, 361 output_filename, 362 is_compact_fc, 363 cutoff_pair_distance, 364 fc_calculator, 365 fc_calculator_options, 366 log_level): 367 """ 368 369 Note 370 ---- 371 cutoff_pair_distance is the parameter to determine each displaced 372 supercell is included to the computation of fc3. It is assumed that 373 cutoff_pair_distance is stored in the step to create sets of 374 displacements and the value is stored n the displacement dataset and 375 also as the parameter 'included': True or False for each displacement. 376 The parameter cutoff_pair_distance here can be used in the step to 377 create fc3 by overwriting original cutoff_pair_distance value only 378 when the former value is smaller than the later. 379 380 """ 381 if input_filename is None: 382 disp_filename = 'disp_fc3.yaml' 383 else: 384 disp_filename = 'disp_fc3.' + input_filename + '.yaml' 385 386 _ph3py_yaml = _get_ph3py_yaml(disp_filename, ph3py_yaml) 387 388 try: 389 dataset = parse_forces(phono3py, 390 ph3py_yaml=_ph3py_yaml, 391 cutoff_pair_distance=cutoff_pair_distance, 392 force_filename="FORCES_FC3", 393 disp_filename=disp_filename, 394 fc_type='fc3', 395 log_level=log_level) 396 except RuntimeError as e: 397 # from _parse_forces_type1 398 if log_level: 399 print(str(e)) 400 print_error() 401 sys.exit(1) 402 except FileNotFoundError as e: 403 # from _get_type2_dataset 404 file_exists(e.filename, log_level) 405 406 phono3py.dataset = dataset 407 phono3py.produce_fc3(symmetrize_fc3r=symmetrize_fc3r, 408 is_compact_fc=is_compact_fc, 409 fc_calculator=fc_calculator, 410 fc_calculator_options=fc_calculator_options) 411 412 413def _create_phono3py_fc2(phono3py, 414 ph3py_yaml, 415 symmetrize_fc2, 416 input_filename, 417 is_compact_fc, 418 fc_calculator, 419 fc_calculator_options, 420 log_level): 421 if input_filename is None: 422 disp_filename = 'disp_fc3.yaml' 423 else: 424 disp_filename = 'disp_fc3.' + input_filename + '.yaml' 425 426 _ph3py_yaml = _get_ph3py_yaml(disp_filename, ph3py_yaml) 427 428 try: 429 dataset = parse_forces(phono3py, 430 ph3py_yaml=_ph3py_yaml, 431 force_filename="FORCES_FC3", 432 disp_filename=disp_filename, 433 fc_type='fc2', 434 log_level=log_level) 435 except RuntimeError as e: 436 if log_level: 437 print(str(e)) 438 print_error() 439 sys.exit(1) 440 except FileNotFoundError as e: 441 file_exists(e.filename, log_level) 442 443 phono3py.phonon_dataset = dataset 444 phono3py.produce_fc2( 445 symmetrize_fc2=symmetrize_fc2, 446 is_compact_fc=is_compact_fc, 447 fc_calculator=fc_calculator, 448 fc_calculator_options=fc_calculator_options) 449 450 451def _get_ph3py_yaml(disp_filename, ph3py_yaml): 452 _ph3py_yaml = ph3py_yaml 453 # Try to use phono3py.phonon_dataset when the disp file not found 454 if not os.path.isfile(disp_filename): 455 disp_filename = None 456 if _ph3py_yaml is None and os.path.isfile('phono3py_disp.yaml'): 457 _ph3py_yaml = Phono3pyYaml() 458 _ph3py_yaml.read('phono3py_disp.yaml') 459 return _ph3py_yaml 460 461 462def _create_phono3py_phonon_fc2(phono3py, 463 ph3py_yaml, 464 symmetrize_fc2, 465 input_filename, 466 is_compact_fc, 467 fc_calculator, 468 fc_calculator_options, 469 log_level): 470 if input_filename is None: 471 disp_filename = 'disp_fc2.yaml' 472 else: 473 disp_filename = 'disp_fc2.' + input_filename + '.yaml' 474 475 _ph3py_yaml = _get_ph3py_yaml(disp_filename, ph3py_yaml) 476 477 try: 478 dataset = parse_forces(phono3py, 479 ph3py_yaml=_ph3py_yaml, 480 force_filename="FORCES_FC2", 481 disp_filename=disp_filename, 482 fc_type='phonon_fc2', 483 log_level=log_level) 484 except RuntimeError as e: 485 if log_level: 486 print(str(e)) 487 print_error() 488 sys.exit(1) 489 except FileNotFoundError as e: 490 file_exists(e.filename, log_level) 491 492 phono3py.phonon_dataset = dataset 493 phono3py.produce_fc2( 494 symmetrize_fc2=symmetrize_fc2, 495 is_compact_fc=is_compact_fc, 496 fc_calculator=fc_calculator, 497 fc_calculator_options=fc_calculator_options) 498 499 500def _convert_unit_in_dataset(dataset, calculator): 501 physical_units = get_default_physical_units(calculator) 502 force_to_eVperA = physical_units['force_to_eVperA'] 503 distance_to_A = physical_units['distance_to_A'] 504 505 if 'first_atoms' in dataset: 506 for d1 in dataset['first_atoms']: 507 if distance_to_A is not None: 508 disp = _to_ndarray(d1['displacement']) 509 d1['displacement'] = disp * distance_to_A 510 if force_to_eVperA is not None and 'forces' in d1: 511 forces = _to_ndarray(d1['forces']) 512 d1['forces'] = forces * force_to_eVperA 513 if 'second_atoms' in d1: 514 for d2 in d1['second_atoms']: 515 if distance_to_A is not None: 516 disp = _to_ndarray(d2['displacement']) 517 d2['displacement'] = disp * distance_to_A 518 if force_to_eVperA is not None and 'forces' in d2: 519 forces = _to_ndarray(d2['forces']) 520 d2['forces'] = forces * force_to_eVperA 521 else: 522 if distance_to_A is not None and 'displacements' in dataset: 523 disp = _to_ndarray(dataset['displacements']) 524 dataset['displacements'] = disp * distance_to_A 525 if force_to_eVperA is not None and 'forces' in dataset: 526 forces = _to_ndarray(dataset['forces']) 527 dataset['forces'] = forces * force_to_eVperA 528 529 530def _to_ndarray(array, dtype='double'): 531 if type(array) is not np.ndarray: 532 return np.array(array, dtype=dtype, order='C') 533 else: 534 return array 535 536 537def _extract_datast_from_ph3py_yaml(ph3py_yaml, fc_type): 538 dataset = None 539 if fc_type == 'phonon_fc2': 540 if ph3py_yaml and ph3py_yaml.phonon_dataset is not None: 541 # copy dataset 542 # otherwise unit conversion can be applied multiple times. 543 _ph3py_yaml = Phono3pyYaml() 544 _ph3py_yaml.yaml_data = ph3py_yaml.yaml_data 545 _ph3py_yaml.parse() 546 dataset = _ph3py_yaml.phonon_dataset 547 else: 548 if ph3py_yaml and ph3py_yaml.dataset is not None: 549 # copy dataset 550 # otherwise unit conversion can be applied multiple times. 551 _ph3py_yaml = Phono3pyYaml() 552 _ph3py_yaml.yaml_data = ph3py_yaml.yaml_data 553 _ph3py_yaml.parse() 554 dataset = _ph3py_yaml.dataset 555 return dataset 556