From 13c401dc761351c476eba841032f6e874b415d7e Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Sun, 28 Apr 2024 22:18:05 +0900 Subject: [PATCH] Write phono3py_params.yaml when parsing calculator result --- phono3py/cui/create_force_sets.py | 355 ++++++++++++++++------------ phono3py/cui/phono3py_script.py | 2 +- phono3py/interface/phono3py_yaml.py | 10 + 3 files changed, 217 insertions(+), 150 deletions(-) diff --git a/phono3py/cui/create_force_sets.py b/phono3py/cui/create_force_sets.py index 2956c78f..5173d4b2 100644 --- a/phono3py/cui/create_force_sets.py +++ b/phono3py/cui/create_force_sets.py @@ -38,9 +38,10 @@ import copy import sys -from typing import Optional, Union +from typing import Optional from phonopy.cui.create_force_sets import check_number_of_force_files +from phonopy.cui.load_helper import get_nac_params from phonopy.cui.phonopy_script import file_exists, files_exist, print_error from phonopy.file_IO import parse_FORCE_SETS, write_FORCE_SETS from phonopy.interface.calculator import get_calc_dataset @@ -60,172 +61,90 @@ def create_FORCES_FC3_and_FORCES_FC2( settings, cell_filename: Optional[str], - log_level: Union[bool, int], + log_level: int = 0, ): - """Create FORCES_FC3 and FORCES_FC2 from files.""" - interface_mode = settings.calculator - ph3py_yaml = None + """Create FORCES_FC3 and FORCES_FC2 or phono3py_params.yaml from files. - ##################### - # Create FORCES_FC3 # - ##################### - if settings.create_forces_fc3 or settings.create_forces_fc3_file: - disp_filename_candidates = [ - "phono3py_disp.yaml", - ] - if cell_filename is not None: - disp_filename_candidates.insert(0, cell_filename) - disp_filenames = files_exist(disp_filename_candidates, log_level, is_any=True) - disp_filename = disp_filenames[0] - ph3py_yaml = Phono3pyYaml() - ph3py_yaml.read(disp_filename) - if ph3py_yaml.calculator is not None: - interface_mode = ph3py_yaml.calculator # overwrite - disp_dataset = ph3py_yaml.dataset + With settings.save_params=True, phono3py_params.yaml is created instead of + FORCES_FC3 and FORCES_FC2. To create phono3py_params.yaml, at least fc3 + forces have to be collected. - if log_level: - print("") - print('Displacement dataset was read from "%s".' % disp_filename) - - num_atoms = disp_dataset["natom"] - num_disps = len(disp_dataset["first_atoms"]) - for d1 in disp_dataset["first_atoms"]: - for d2 in d1["second_atoms"]: - if "included" not in d2 or d2["included"]: - num_disps += 1 - - if settings.create_forces_fc3_file: - file_exists(settings.create_forces_fc3_file, log_level) - force_filenames = [x.strip() for x in open(settings.create_forces_fc3_file)] - else: - force_filenames = settings.create_forces_fc3 - - for filename in force_filenames: - file_exists(filename, log_level) - - if log_level > 0: - print("Number of displacements: %d" % num_disps) - print("Number of supercell files: %d" % len(force_filenames)) + """ + interface_mode = settings.calculator + disp_filename_candidates = [ + "phono3py_disp.yaml", + ] + if cell_filename is not None: + disp_filename_candidates.insert(0, cell_filename) + disp_filenames = files_exist(disp_filename_candidates, log_level, is_any=True) + disp_filename = disp_filenames[0] + ph3py_yaml = Phono3pyYaml(settings={"force_sets": True}) + ph3py_yaml.read(disp_filename) + if ph3py_yaml.calculator is not None: + interface_mode = ph3py_yaml.calculator # overwrite - if not check_number_of_force_files(num_disps, force_filenames, disp_filename): - force_sets = [] - else: - calc_dataset = get_calc_dataset( - interface_mode, - num_atoms, - force_filenames, - verbose=(log_level > 0), - ) - force_sets = calc_dataset["forces"] - - if settings.subtract_forces: - force_filename = settings.subtract_forces - file_exists(force_filename, log_level) - calc_dataset = get_calc_dataset( - interface_mode, - num_atoms, - [ - force_filename, - ], - verbose=(log_level > 0), - ) - force_set_zero = calc_dataset["forces"][0] - for fs in force_sets: - fs -= force_set_zero - - if log_level > 0: - print( - "Forces in '%s' were subtracted from supercell forces." - % force_filename - ) - - if force_sets: - write_FORCES_FC3(disp_dataset, forces_fc3=force_sets, filename="FORCES_FC3") - if log_level: - print("") - print("%s has been created." % "FORCES_FC3") - print("") - else: + if settings.create_forces_fc3 or settings.create_forces_fc3_file: + calc_dataset_fc3 = _get_force_sets_fc3( + settings, ph3py_yaml.dataset, disp_filename, interface_mode, log_level + ) + if not calc_dataset_fc3["forces"]: if log_level: - print("") print("%s could not be created." % "FORCES_FC3") print_error() sys.exit(1) - ##################### - # Create FORCES_FC2 # - ##################### if settings.create_forces_fc2: - disp_filename_candidates = [ - "phono3py_disp.yaml", - ] - if cell_filename is not None: - disp_filename_candidates.insert(0, cell_filename) - disp_filenames = files_exist(disp_filename_candidates, log_level, is_any=True) - disp_filename = disp_filenames[0] - - # ph3py_yaml is not None, phono3py_disp.yaml is already read. - if ph3py_yaml is None: - ph3py_yaml = Phono3pyYaml() - ph3py_yaml.read(disp_filename) - if ph3py_yaml.calculator is not None: - interface_mode = ph3py_yaml.calculator # overwrite - disp_dataset = ph3py_yaml.phonon_dataset - - if log_level: - print('Displacement dataset was read from "%s".' % disp_filename) - num_atoms = disp_dataset["natom"] - num_disps = len(disp_dataset["first_atoms"]) - force_filenames = settings.create_forces_fc2 - for filename in force_filenames: - file_exists(filename, log_level) - - if log_level > 0: - print("Number of displacements: %d" % num_disps) - print("Number of supercell files: %d" % len(force_filenames)) - - calc_dataset = get_calc_dataset( + calc_dataset_fc2 = _get_force_sets_fc2( + settings, + ph3py_yaml.phonon_dataset, + disp_filename, interface_mode, - num_atoms, - force_filenames, - verbose=(log_level > 0), + log_level, ) - force_sets = calc_dataset["forces"] - - if settings.subtract_forces: - force_filename = settings.subtract_forces - file_exists(force_filename, log_level) - calc_dataset_zero = get_calc_dataset( - interface_mode, - num_atoms, - [ - force_filename, - ], - verbose=(log_level > 0), - ) - force_set_zero = calc_dataset_zero["forces"][0] - for fs in force_sets: - fs -= force_set_zero - - if log_level > 0: - print( - "Forces in '%s' were subtracted from supercell forces." - % force_filename - ) - - if force_sets: - write_FORCES_FC2(disp_dataset, forces_fc2=force_sets, filename="FORCES_FC2") + if not calc_dataset_fc2["forces"]: if log_level: - print("") - print("%s has been created." % "FORCES_FC2") - print("") - else: + print("%s could not be created." % "FORCES_FC2") + print_error() + sys.exit(1) + + if settings.save_params: + fc3_yaml_filename = "phono3py_params.yaml" + if not (settings.create_forces_fc3 or settings.create_forces_fc3_file): if log_level: + print(f'When creating "{fc3_yaml_filename}" with force sets for fc2, ') + print("force sets for fc3 have to be collected simultaneously.") + print(f'"{fc3_yaml_filename}" could not be created.') print("") - print("%s could not be created." % "FORCES_FC2") print_error() sys.exit(1) + _set_forces_and_nac_params( + ph3py_yaml, settings, calc_dataset_fc3, calc_dataset_fc2 + ) + + with open(fc3_yaml_filename, "w") as w: + w.write(str(ph3py_yaml)) + if log_level: + print(f'"{fc3_yaml_filename}" has been created.') + else: + if settings.create_forces_fc3 or settings.create_forces_fc3_file: + write_FORCES_FC3( + ph3py_yaml.dataset, + forces_fc3=calc_dataset_fc3["forces"], + filename="FORCES_FC3", + ) + if log_level: + print("%s has been created." % "FORCES_FC3") + + if settings.create_forces_fc2: + write_FORCES_FC2( + ph3py_yaml.phonon_dataset, + forces_fc2=calc_dataset_fc2["forces"], + filename="FORCES_FC2", + ) + if log_level: + print("%s has been created." % "FORCES_FC2") + def create_FORCES_FC2_from_FORCE_SETS(log_level): """Convert FORCE_SETS to FORCES_FC2.""" @@ -300,3 +219,141 @@ def create_FORCE_SETS_from_FORCES_FCx( "The file format of %s is already readable by phonopy." % forces_filename ) + + +def _get_force_sets_fc2( + settings, disp_dataset, disp_filename, interface_mode, log_level +) -> dict: + interface_mode = settings.calculator + if log_level: + print(f'FC2 displacement dataset was read from "{disp_filename}".') + num_atoms = disp_dataset["natom"] + num_disps = len(disp_dataset["first_atoms"]) + force_filenames = settings.create_forces_fc2 + for filename in force_filenames: + file_exists(filename, log_level) + + if log_level > 0: + print(f" Number of displacements: {num_disps}") + print(f" Number of supercell files: {len(force_filenames)}") + + calc_dataset = get_calc_dataset( + interface_mode, + num_atoms, + force_filenames, + verbose=(log_level > 0), + ) + force_sets = calc_dataset["forces"] + + if settings.subtract_forces: + force_filename = settings.subtract_forces + file_exists(force_filename, log_level) + calc_dataset_zero = get_calc_dataset( + interface_mode, + num_atoms, + [ + force_filename, + ], + verbose=(log_level > 0), + ) + force_set_zero = calc_dataset_zero["forces"][0] + for fs in force_sets: + fs -= force_set_zero + + if log_level > 0: + print("Forces in '{force_filename}' were subtracted from supercell forces.") + + if log_level > 0: + print("") + + return calc_dataset + + +def _get_force_sets_fc3( + settings, disp_dataset, disp_filename, interface_mode, log_level +) -> dict: + if log_level: + print("") + print(f'FC3 Displacement dataset was read from "{disp_filename}".') + + num_atoms = disp_dataset["natom"] + num_disps = len(disp_dataset["first_atoms"]) + for d1 in disp_dataset["first_atoms"]: + for d2 in d1["second_atoms"]: + if "included" not in d2 or d2["included"]: + num_disps += 1 + + if settings.create_forces_fc3_file: + file_exists(settings.create_forces_fc3_file, log_level) + force_filenames = [x.strip() for x in open(settings.create_forces_fc3_file)] + else: + force_filenames = settings.create_forces_fc3 + + for filename in force_filenames: + file_exists(filename, log_level) + + if log_level > 0: + print(f" Number of displacements: {num_disps}") + print(f" Number of supercell files: {len(force_filenames)}") + + if not check_number_of_force_files(num_disps, force_filenames, disp_filename): + calc_dataset = {"forces": []} + else: + calc_dataset = get_calc_dataset( + interface_mode, + num_atoms, + force_filenames, + verbose=(log_level > 0), + ) + force_sets = calc_dataset["forces"] + + if settings.subtract_forces: + force_filename = settings.subtract_forces + file_exists(force_filename, log_level) + calc_dataset = get_calc_dataset( + interface_mode, + num_atoms, + [ + force_filename, + ], + verbose=(log_level > 0), + ) + force_set_zero = calc_dataset["forces"][0] + for fs in force_sets: + fs -= force_set_zero + + if log_level > 0: + print( + f"Forces in '{force_filename}' were subtracted from supercell forces." + ) + + if log_level > 0: + print("") + + return calc_dataset + + +def _set_forces_and_nac_params( + ph3py_yaml: Phono3pyYaml, settings, calc_dataset_fc3: dict, calc_dataset_fc2: dict +): + count = len(ph3py_yaml.dataset["first_atoms"]) + for i, d1 in enumerate(ph3py_yaml.dataset["first_atoms"]): + d1["forces"] = calc_dataset_fc3["forces"][i] + if "energies" in calc_dataset_fc3: + d1["energy"] = float(calc_dataset_fc3["energies"][i]) + for d2 in d1["second_atoms"]: + if "included" not in d2 or d2["included"]: + d2["forces"] = calc_dataset_fc3["forces"][count] + if "energies" in calc_dataset_fc3: + d2["energy"] = float(calc_dataset_fc3["energies"][count]) + count += 1 + + if settings.create_forces_fc2: + for i, d in enumerate(ph3py_yaml.phonon_dataset["first_atoms"]): + d["forces"] = calc_dataset_fc2["forces"][i] + if "energies" in calc_dataset_fc2: + d["energy"] = float(calc_dataset_fc2["energies"][i]) + + nac_params = get_nac_params(primitive=ph3py_yaml.primitive) + if nac_params: + ph3py_yaml.nac_params = nac_params diff --git a/phono3py/cui/phono3py_script.py b/phono3py/cui/phono3py_script.py index b4954eae..0850d0b2 100644 --- a/phono3py/cui/phono3py_script.py +++ b/phono3py/cui/phono3py_script.py @@ -879,7 +879,7 @@ def main(**argparse_control): or settings.create_forces_fc3_file or settings.create_forces_fc2 ): - create_FORCES_FC3_and_FORCES_FC2(settings, cell_filename, log_level) + create_FORCES_FC3_and_FORCES_FC2(settings, cell_filename, log_level=log_level) if log_level: print_end_phono3py() sys.exit(0) diff --git a/phono3py/interface/phono3py_yaml.py b/phono3py/interface/phono3py_yaml.py index 79c03256..81d972f0 100644 --- a/phono3py/interface/phono3py_yaml.py +++ b/phono3py/interface/phono3py_yaml.py @@ -192,6 +192,8 @@ def _parse_fc3_dataset_type1(self, natom): data1["id"] = d1_id if "forces" in d1: data1["forces"] = np.array(d1["forces"], dtype="double", order="C") + if "energy" in d1: + data1["energy"] = d1["energy"] d2_list = d1.get("paired_with") if d2_list is None: # backward compatibility d2_list = d1.get("second_atoms") @@ -226,6 +228,8 @@ def _parse_fc3_dataset_type1_with_forces(self, data1, d2, disp2_id): second_atom_dict["id"] = d2_id if "pair_distance" in d2: second_atom_dict["pair_distance"] = d2["pair_distance"] + if "energy" in d2: + second_atom_dict["energy"] = d2["energy"] disp2_id += 1 data1["second_atoms"].append(second_atom_dict) return disp2_id @@ -464,6 +468,8 @@ def displacements_yaml_lines_type1(dataset, with_forces=False, key="displacement lines.append(" forces:") for v in d1["forces"]: lines.append(" - [ %19.16f, %19.16f, %19.16f ]" % tuple(v)) + if "energy" in d1: + lines.append(" energy: {energy:.8f}".format(energy=d1["energy"])) if "second_atoms" in d1: ret_lines, id_offset = _second_displacements_yaml_lines( d1["second_atoms"], id_offset, with_forces=with_forces @@ -542,6 +548,10 @@ def _second_displacements_yaml_lines(dataset2, id_offset, with_forces=False): lines.append(" forces:") for v in dataset2[j]["forces"]: lines.append(" - [ %19.16f, %19.16f, %19.16f ]" % tuple(v)) + if "energy" in dataset2[j]: + lines.append( + " energy: {energy:.8f}".format(energy=dataset2[j]["energy"]) + ) else: lines.append(" - atom: %4d" % (i + 1)) lines.append(