diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index ab24f699..e8ee8feb 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -19,7 +19,7 @@ repos: - "--ignore=E203,W503" - repo: https://github.com/psf/black - rev: 24.4.0 + rev: 24.4.2 hooks: - id: black args: diff --git a/doc/changelog.md b/doc/changelog.md index 20942849..3bd42254 100644 --- a/doc/changelog.md +++ b/doc/changelog.md @@ -2,6 +2,10 @@ # Change Log +## May-4-2024: Version 3.0.3 + +- Release to follow the update of phonopy. + ## Apr-21-2024: Version 3.0.2 - New way of ph-ph interaction calculation (see {ref}`changelog_v290`) is used as diff --git a/doc/conf.py b/doc/conf.py index 2d38dbd6..e4c57601 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -60,7 +60,7 @@ # The short X.Y version. version = "3.0" # The full version, including alpha/beta/rc tags. -release = "3.0.2" +release = "3.0.3" # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/phono3py/api_phono3py.py b/phono3py/api_phono3py.py index 8de413d4..34b6be28 100644 --- a/phono3py/api_phono3py.py +++ b/phono3py/api_phono3py.py @@ -34,8 +34,9 @@ # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. +import copy from collections.abc import Sequence -from typing import Optional, Union +from typing import Literal, Optional, Union import numpy as np from phonopy.exception import ForceCalculatorRequiredError @@ -582,6 +583,7 @@ def dataset(self): {'number': atom index of second displaced atom, 'displacement': displacement in Cartesian coordinates}, 'forces': forces on atoms in supercell, + 'supercell_energy': energy of supercell, 'pair_distance': distance between paired atoms, 'included': with cutoff pair distance in displacement pair generation, this indicates if this @@ -591,9 +593,11 @@ def dataset(self): ... ] }, ... ] } Type 2. All atomic displacements in each supercell: {'displacements': ndarray, dtype='double', order='C', - shape=(supercells, atoms in supercell, 3) + shape=(supercells, atoms in supercell, 3), 'forces': ndarray, dtype='double',, order='C', - shape=(supercells, atoms in supercell, 3)} + shape=(supercells, atoms in supercell, 3), + 'supercell_energies': ndarray, dtype='double', + shape=(supercells,)} In type 2, displacements and forces can be given by numpy array with different shape but that can be reshaped to (supercells, natom, 3). @@ -607,7 +611,20 @@ def dataset(self): @dataset.setter def dataset(self, dataset): - self._dataset = dataset + if dataset is None: + self._dataset = None + elif "first_atoms" in dataset: + self._dataset = copy.deepcopy(dataset) + elif "displacements" in dataset: + self._dataset = {} + self.displacements = dataset["displacements"] + if "forces" in dataset: + self.forces = dataset["forces"] + if "supercell_energies" in dataset: + self.supercell_energies = dataset["supercell_energies"] + else: + raise RuntimeError("Data format of dataset is wrong.") + self._supercells_with_displacements = None self._phonon_supercells_with_displacements = None @@ -622,12 +639,15 @@ def phonon_dataset(self): 'first_atoms': [ {'number': atom index of first displaced atom, 'displacement': displacement in Cartesian coordinates, - 'forces': forces on atoms in supercell} ... ]} + 'forces': forces on atoms in supercell, + 'supercell_energy': energy of supercell}, ... ]} Type 2. All atomic displacements in each supercell: {'displacements': ndarray, dtype='double', order='C', - shape=(supercells, atoms in supercell, 3) + shape=(supercells, atoms in supercell, 3), 'forces': ndarray, dtype='double',, order='C', - shape=(supercells, atoms in supercell, 3)} + shape=(supercells, atoms in supercell, 3), + 'supercell_energies': ndarray, dtype='double', + shape=(supercells,)} In type 2, displacements and forces can be given by numpy array with different shape but that can be reshaped to (supercells, natom, 3). @@ -637,7 +657,21 @@ def phonon_dataset(self): @phonon_dataset.setter def phonon_dataset(self, dataset): - self._phonon_dataset = dataset + if dataset is None: + self._phonon_dataset = None + elif "first_atoms" in dataset: + self._phonon_dataset = copy.deepcopy(dataset) + elif "displacements" in dataset: + self._phonon_dataset = {} + self.phonon_displacements = dataset["displacements"] + if "forces" in dataset: + self.phonon_forces = dataset["forces"] + if "supercell_energies" in dataset: + self.phonon_supercell_energies = dataset["supercell_energies"] + else: + raise RuntimeError("Data format of dataset is wrong.") + + self._phonon_supercells_with_displacements = None @property def band_indices(self): @@ -816,45 +850,32 @@ def forces(self): be the same order of supercells_with_displacements. """ - dataset = self._dataset - if "forces" in dataset: - return dataset["forces"] - elif "first_atoms" in dataset: - num_scells = len(dataset["first_atoms"]) - for disp1 in dataset["first_atoms"]: - num_scells += len(disp1["second_atoms"]) - forces = np.zeros( - (num_scells, self._supercell.get_number_of_atoms(), 3), - dtype="double", - order="C", - ) - i = 0 - for disp1 in dataset["first_atoms"]: - forces[i] = disp1["forces"] - i += 1 - for disp1 in dataset["first_atoms"]: - for disp2 in disp1["second_atoms"]: - forces[i] = disp2["forces"] - i += 1 - return forces - else: - raise RuntimeError("displacement dataset has wrong format.") + return self._get_forces_energies(target="forces") @forces.setter - def forces(self, forces_fc3): - forces = np.array(forces_fc3, dtype="double", order="C") - dataset = self._dataset - if "first_atoms" in dataset: - i = 0 - for disp1 in dataset["first_atoms"]: - disp1["forces"] = forces[i] - i += 1 - for disp1 in dataset["first_atoms"]: - for disp2 in disp1["second_atoms"]: - disp2["forces"] = forces[i] - i += 1 - elif "displacements" in dataset or "forces" in dataset: - dataset["forces"] = forces + def forces(self, values): + self._set_forces_energies(values, target="forces") + + @property + def supercell_energies(self): + """Setter and getter of supercell energies in displacement dataset. + + A set of supercell energies of displaced supercells. The order of + displaced supercells has to match with that in displacement dataset. + shape=(displaced supercells,) + + getter : ndarray + + setter : array_like + The order of supercells used for calculating supercell energies has + to be the same order of supercells_with_displacements. + + """ + return self._get_forces_energies(target="supercell_energies") + + @supercell_energies.setter + def supercell_energies(self, values): + self._set_forces_energies(values, target="supercell_energies") @property def phonon_displacements(self): @@ -922,7 +943,7 @@ def phonon_displacements(self, displacements): @property def phonon_forces(self): - """Setter and getter of forces in displacement dataset for fc2. + """Setter and getter of forces in fc2 displacement dataset. A set of atomic forces in displaced supercells. The order of displaced supercells has to match with that in phonon displacement @@ -936,39 +957,30 @@ def phonon_forces(self): be the same order of phonon_supercells_with_displacements. """ - if self._phonon_dataset is None: - raise RuntimeError("phonon_displacement_dataset does not exist.") - - dataset = self._phonon_dataset - if "forces" in dataset: - return dataset["forces"] - elif "first_atoms" in dataset: - num_scells = len(dataset["first_atoms"]) - forces = np.zeros( - (num_scells, self._phonon_supercell.get_number_of_atoms(), 3), - dtype="double", - order="C", - ) - for i, disp1 in enumerate(dataset["first_atoms"]): - forces[i] = disp1["forces"] - return forces - else: - raise RuntimeError("displacement dataset has wrong format.") + return self._get_phonon_forces_energies(target="forces") @phonon_forces.setter - def phonon_forces(self, forces_fc2): - if self._phonon_dataset is None: - raise RuntimeError("phonon_displacement_dataset does not exist.") + def phonon_forces(self, values): + self._set_phonon_forces_energies(values, target="forces") - forces = np.array(forces_fc2, dtype="double", order="C") - dataset = self._phonon_dataset - if "first_atoms" in dataset: - i = 0 - for i, disp1 in enumerate(dataset["first_atoms"]): - disp1["forces"] = forces[i] - i += 1 - elif "displacements" in dataset or "forces" in dataset: - dataset["forces"] = forces + @property + def phonon_supercell_energies(self): + """Setter and getter of supercell energies in fc2 displacement dataset. + + shape=(displaced supercells,) + + getter : ndarray + + setter : array_like + The order of supercells used for calculating supercell energies has + to be the same order of phonon_supercells_with_displacements. + + """ + return self._get_phonon_forces_energies(target="supercell_energies") + + @phonon_supercell_energies.setter + def phonon_supercell_energies(self, values): + self._set_phonon_forces_energies(values, target="supercell_energies") @property def phph_interaction(self): @@ -2274,3 +2286,116 @@ def _extract_fc2_fc3_calculators(self, fc_calculator, fc_calculator_options, ord _fc_calculator_options = None return _fc_calculator, _fc_calculator_options + + def _get_forces_energies( + self, target: Literal["forces", "supercell_energies"] + ) -> Optional[np.ndarray]: + """Return fc3 forces and supercell energies. + + Return None if tagert data is not found rather than raising exception. + + """ + if target in self._dataset: # type-2 + return self._dataset[target] + elif "first_atoms" in self._dataset: # type-1 + num_scells = len(self._dataset["first_atoms"]) + for disp1 in self._dataset["first_atoms"]: + num_scells += len(disp1["second_atoms"]) + if target == "forces": + values = np.zeros( + (num_scells, len(self._supercell), 3), + dtype="double", + order="C", + ) + type1_target = "forces" + elif target == "supercell_energies": + values = np.zeros(num_scells, dtype="double") + type1_target = "supercell_energy" + count = 0 + for disp1 in self._dataset["first_atoms"]: + values[count] = disp1[type1_target] + count += 1 + for disp1 in self._dataset["first_atoms"]: + for disp2 in disp1["second_atoms"]: + values[count] = disp2[type1_target] + count += 1 + return values + return None + + def _set_forces_energies( + self, values, target: Literal["forces", "supercell_energies"] + ): + if "first_atoms" in self._dataset: # type-1 + count = 0 + for disp1 in self._dataset["first_atoms"]: + if target == "forces": + disp1[target] = np.array(values[count], dtype="double", order="C") + elif target == "supercell_energies": + disp1["supercell_energy"] = float(values[count]) + count += 1 + for disp1 in self._dataset["first_atoms"]: + for disp2 in disp1["second_atoms"]: + if target == "forces": + disp2[target] = np.array( + values[count], dtype="double", order="C" + ) + elif target == "supercell_energies": + disp2["supercell_energy"] = float(values[count]) + count += 1 + elif "displacements" in self._dataset or "forces" in self._dataset: # type-2 + self._dataset[target] = np.array(values, dtype="double", order="C") + else: + raise RuntimeError("Set of FC3 displacements is not available.") + + def _get_phonon_forces_energies( + self, target: Literal["forces", "supercell_energies"] + ) -> Optional[np.ndarray]: + """Return fc2 forces and supercell energies. + + Return None if tagert data is not found rather than raising exception. + + """ + if self._phonon_dataset is None: + raise RuntimeError("Dataset for fc2does not exist.") + + if target in self._phonon_dataset: # type-2 + return self._phonon_dataset[target] + elif "first_atoms" in self._phonon_dataset: # type-1 + values = [] + for disp in self._phonon_dataset["first_atoms"]: + if target == "forces": + if target in disp: + values.append(disp[target]) + elif target == "supercell_energies": + if "supercell_energy" in disp: + values.append(disp["supercell_energy"]) + if values: + return np.array(values, dtype="double", order="C") + return None + + def _set_phonon_forces_energies( + self, values, target: Literal["forces", "supercell_energies"] + ): + if self._phonon_dataset is None: + raise RuntimeError("Dataset for fc2 does not exist.") + + if "first_atoms" in self._phonon_dataset: + for disp, v in zip(self._phonon_dataset["first_atoms"], values): + if target == "forces": + disp[target] = np.array(v, dtype="double", order="C") + elif target == "supercell_energies": + disp["supercell_energy"] = float(v) + elif "displacements" in self._phonon_dataset: + _values = np.array(values, dtype="double", order="C") + natom = len(self._phonon_supercell) + ndisps = len(self._phonon_dataset["displacements"]) + if target == "forces" and ( + _values.ndim != 3 or _values.shape != (ndisps, natom, 3) + ): + raise RuntimeError(f"Array shape of input {target} is incorrect.") + elif target == "supercell_energies": + if _values.ndim != 1 or _values.shape != (ndisps,): + raise RuntimeError(f"Array shape of input {target} is incorrect.") + self._phonon_dataset[target] = _values + else: + raise RuntimeError("Set of FC2 displacements is not available.") diff --git a/phono3py/cui/create_force_sets.py b/phono3py/cui/create_force_sets.py new file mode 100644 index 00000000..8916aa17 --- /dev/null +++ b/phono3py/cui/create_force_sets.py @@ -0,0 +1,396 @@ +"""Command line user interface to create force sets files.""" + +# Copyright (C) 2024 Atsushi Togo +# All rights reserved. +# +# This file is part of phono3py. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions +# are met: +# +# * Redistributions of source code must retain the above copyright +# notice, this list of conditions and the following disclaimer. +# +# * Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in +# the documentation and/or other materials provided with the +# distribution. +# +# * Neither the name of the phonopy project nor the names of its +# contributors may be used to endorse or promote products derived +# from this software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS +# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE +# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. + +from __future__ import annotations + +import copy +import sys +from typing import Optional + +import numpy as np +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 + +from phono3py.file_IO import ( + get_length_of_first_line, + parse_FORCES_FC2, + write_FORCES_FC2, + write_FORCES_FC3, +) +from phono3py.interface.phono3py_yaml import ( + Phono3pyYaml, + displacements_yaml_lines_type1, +) + + +def create_FORCES_FC3_and_FORCES_FC2( + settings, + cell_filename: Optional[str], + log_level: int = 0, +): + """Create FORCES_FC3 and FORCES_FC2 or phono3py_params.yaml from files. + + 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. + + """ + 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 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("%s could not be created." % "FORCES_FC3") + print_error() + sys.exit(1) + + if settings.create_forces_fc2: + calc_dataset_fc2 = _get_force_sets_fc2( + settings, + ph3py_yaml.phonon_dataset, + disp_filename, + interface_mode, + log_level, + ) + if not calc_dataset_fc2["forces"]: + if log_level: + 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_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.""" + filename = "FORCE_SETS" + file_exists(filename, log_level) + disp_dataset = parse_FORCE_SETS(filename=filename) + write_FORCES_FC2(disp_dataset) + + if log_level: + print("") + print("FORCES_FC2 has been created from FORCE_SETS.") + print("The following yaml lines should replace respective part of") + print("phono3py_disp.yaml made with --dim-fc2=dim_of_FORCE_SETS.") + + print("") + print("\n".join(displacements_yaml_lines_type1(disp_dataset))) + + +def create_FORCE_SETS_from_FORCES_FCx( + phonon_smat, input_filename: Optional[str], cell_filename: Optional[str], log_level +): + """Convert FORCES_FC3 or FORCES_FC2 to FORCE_SETS.""" + if cell_filename is not None: + disp_filename = cell_filename + elif input_filename is None: + disp_filename = "phono3py_disp.yaml" + else: + disp_filename = f"phono3py_disp.{input_filename}.yaml" + if phonon_smat is not None: + forces_filename = "FORCES_FC2" + else: + forces_filename = "FORCES_FC3" + + if log_level: + print(f'Displacement dataset is read from "{disp_filename}".') + print(f'Forces are read from "{forces_filename}"') + + with open(forces_filename, "r") as f: + len_first_line = get_length_of_first_line(f) + + if len_first_line == 3: + file_exists(disp_filename, log_level) + file_exists(forces_filename, log_level) + ph3yml = Phono3pyYaml() + ph3yml.read(disp_filename) + if phonon_smat is None: + dataset = copy.deepcopy(ph3yml.dataset) + smat = ph3yml.supercell_matrix + else: + dataset = copy.deepcopy(ph3yml.phonon_dataset) + smat = ph3yml.phonon_supercell_matrix + + if smat is None or (phonon_smat is not None and (phonon_smat != smat).any()): + if log_level: + print("") + print("Supercell matrix is inconsistent.") + print(f'Supercell matrix read from "{disp_filename}":') + print(smat) + print("Supercell matrix given by --dim-fc2:") + print(phonon_smat) + print_error() + sys.exit(1) + + parse_FORCES_FC2(dataset, filename=forces_filename) + write_FORCE_SETS(dataset) + + if log_level: + print("FORCE_SETS has been created.") + else: + if log_level: + print( + "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}".') + + if "first_atoms" in disp_dataset: # type-1 + num_atoms = disp_dataset["natom"] + num_disps = len(disp_dataset["first_atoms"]) + elif "displacements" in disp_dataset: # type-2: + num_disps = len(disp_dataset["displacements"]) + num_atoms = len(disp_dataset["displacements"][0]) + else: + raise RuntimeError("FC2 displacement dataset is broken.") + + 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}".') + + if "first_atoms" in disp_dataset: # type-1 + 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 + elif "displacements" in disp_dataset: # type-2: + num_disps = len(disp_dataset["displacements"]) + num_atoms = len(disp_dataset["displacements"][0]) + else: + raise RuntimeError("FC3 displacement dataset is broken.") + + 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 +): + if "first_atoms" in ph3py_yaml.dataset: + 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 "supercell_energies" in calc_dataset_fc3: + d1["supercell_energy"] = float( + calc_dataset_fc3["supercell_energies"][i] + ) + for d2 in d1["second_atoms"]: + if "included" not in d2 or d2["included"]: + d2["forces"] = calc_dataset_fc3["forces"][count] + if "supercell_energies" in calc_dataset_fc3: + d2["supercell_energy"] = float( + calc_dataset_fc3["supercell_energies"][count] + ) + count += 1 + else: + ph3py_yaml.dataset["forces"] = np.array( + calc_dataset_fc3["forces"], dtype="double", order="C" + ) + ph3py_yaml.dataset["supercell_energies"] = np.array( + calc_dataset_fc3["supercell_energies"], dtype="double" + ) + + if settings.create_forces_fc2: + if "first_atoms" in ph3py_yaml.phonon_dataset: + for i, d in enumerate(ph3py_yaml.phonon_dataset["first_atoms"]): + d["forces"] = calc_dataset_fc2["forces"][i] + if "supercell_energies" in calc_dataset_fc2: + d["supercell_energy"] = float( + calc_dataset_fc2["supercell_energies"][i] + ) + else: + ph3py_yaml.phonon_dataset["forces"] = np.array( + calc_dataset_fc2["forces"], dtype="double", order="C" + ) + ph3py_yaml.phonon_dataset["supercell_energies"] = np.array( + calc_dataset_fc2["supercell_energies"], dtype="double" + ) + + nac_params = get_nac_params(primitive=ph3py_yaml.primitive) + if nac_params: + ph3py_yaml.nac_params = nac_params diff --git a/phono3py/cui/create_supercells.py b/phono3py/cui/create_supercells.py index 1b6cbbc4..3d2e3e84 100644 --- a/phono3py/cui/create_supercells.py +++ b/phono3py/cui/create_supercells.py @@ -1,4 +1,4 @@ -"""Utilities of main CUI script.""" +"""Command line user interface to create supercells.""" # Copyright (C) 2015 Atsushi Togo # All rights reserved. @@ -48,7 +48,6 @@ def create_phono3py_supercells( cell_info, settings, symprec, - output_filename=None, interface_mode="vasp", log_level=1, ): @@ -121,7 +120,7 @@ def create_phono3py_supercells( ) write_supercells_with_displacements( interface_mode, - phono3py.supercell, + phono3py.phonon_supercell, phono3py.phonon_supercells_with_displacements, optional_structure_info, zfill_width=5, diff --git a/phono3py/cui/load.py b/phono3py/cui/load.py index d9ae3c8a..820bb157 100644 --- a/phono3py/cui/load.py +++ b/phono3py/cui/load.py @@ -390,13 +390,14 @@ def set_dataset_and_force_constants( cutoff_pair_distance=cutoff_pair_distance, log_level=log_level, ) - read_fc["fc2"] = _set_dataset_phonon_dataset_or_fc2( - ph3py, - ph3py_yaml=ph3py_yaml, - fc2_filename=fc2_filename, - forces_fc2_filename=forces_fc2_filename, - log_level=log_level, - ) + if ph3py.phonon_supercell_matrix is not None: + read_fc["fc2"] = _set_dataset_phonon_dataset_or_fc2( + ph3py, + ph3py_yaml=ph3py_yaml, + fc2_filename=fc2_filename, + forces_fc2_filename=forces_fc2_filename, + log_level=log_level, + ) # Cases that dataset is in phono3py.yaml but not forces. if ph3py.dataset is None: diff --git a/phono3py/cui/phono3py_argparse.py b/phono3py/cui/phono3py_argparse.py index 6c5d2e74..960f8639 100644 --- a/phono3py/cui/phono3py_argparse.py +++ b/phono3py/cui/phono3py_argparse.py @@ -681,6 +681,14 @@ def get_parser(fc_symmetry=False, is_nac=False, load_phono3py_yaml=False): default=False, help="Calculate real part of self energy", ) + parser.add_argument( + "--sp", + "--save-params", + dest="save_params", + action="store_true", + default=None, + help="Save parameters that can run phono3py in phono3py_params.yaml.", + ) parser.add_argument( "--scattering-event-class", dest="scattering_event_class", diff --git a/phono3py/cui/phono3py_script.py b/phono3py/cui/phono3py_script.py index 138e6e5d..0850d0b2 100644 --- a/phono3py/cui/phono3py_script.py +++ b/phono3py/cui/phono3py_script.py @@ -36,18 +36,14 @@ from __future__ import annotations -import copy import datetime import sys -from typing import Optional, Union import numpy as np from phonopy.cui.collect_cell_info import collect_cell_info -from phonopy.cui.create_force_sets import check_number_of_force_files from phonopy.cui.phonopy_argparse import show_deprecated_option_warnings from phonopy.cui.phonopy_script import ( file_exists, - files_exist, print_end, print_error, print_error_message, @@ -56,9 +52,9 @@ store_nac_params, ) from phonopy.exception import ForceCalculatorRequiredError -from phonopy.file_IO import is_file_phonopy_yaml, parse_FORCE_SETS, write_FORCE_SETS +from phonopy.file_IO import is_file_phonopy_yaml from phonopy.harmonic.force_constants import show_drift_force_constants -from phonopy.interface.calculator import get_default_physical_units, get_force_sets +from phonopy.interface.calculator import get_default_physical_units from phonopy.phonon.band_structure import get_band_qpoints from phonopy.structure.cells import isclose as cells_isclose from phonopy.units import Bohr, Hartree, VaspToTHz @@ -68,6 +64,11 @@ create_phono3py_force_constants, get_fc_calculator_params, ) +from phono3py.cui.create_force_sets import ( + create_FORCE_SETS_from_FORCES_FCx, + create_FORCES_FC2_from_FORCE_SETS, + create_FORCES_FC3_and_FORCES_FC2, +) from phono3py.cui.create_supercells import create_phono3py_supercells from phono3py.cui.load import ( compute_force_constants_from_datasets, @@ -82,19 +83,12 @@ ) from phono3py.cui.triplets_info import show_num_triplets, write_grid_points from phono3py.file_IO import ( - get_length_of_first_line, - parse_FORCES_FC2, read_phonon_from_hdf5, write_fc2_to_hdf5, write_fc3_to_hdf5, - write_FORCES_FC2, - write_FORCES_FC3, write_phonon_to_hdf5, ) -from phono3py.interface.phono3py_yaml import ( - Phono3pyYaml, - displacements_yaml_lines_type1, -) +from phono3py.interface.phono3py_yaml import Phono3pyYaml from phono3py.phonon3.fc3 import show_drift_fc3 from phono3py.phonon3.gruneisen import run_gruneisen_parameters from phono3py.phonon.grid import get_grid_point_from_address, get_ir_grid_points @@ -299,255 +293,6 @@ def read_phono3py_settings(args, argparse_control, log_level): return settings, confs_dict, cell_filename -def create_FORCES_FC2_from_FORCE_SETS_then_exit(log_level): - """Convert FORCE_SETS to FORCES_FC2.""" - filename = "FORCE_SETS" - file_exists(filename, log_level) - disp_dataset = parse_FORCE_SETS(filename=filename) - write_FORCES_FC2(disp_dataset) - - if log_level: - print("") - print("FORCES_FC2 has been created from FORCE_SETS.") - print("The following yaml lines should replace respective part of") - print("phono3py_disp.yaml made with --dim-fc2=dim_of_FORCE_SETS.") - - print("") - print("\n".join(displacements_yaml_lines_type1(disp_dataset))) - - if log_level: - print_end_phono3py() - sys.exit(0) - - -def create_FORCE_SETS_from_FORCES_FCx_then_exit( - phonon_smat, input_filename: Optional[str], cell_filename: Optional[str], log_level -): - """Convert FORCES_FC3 or FORCES_FC2 to FORCE_SETS.""" - if cell_filename is not None: - disp_filename = cell_filename - elif input_filename is None: - disp_filename = "phono3py_disp.yaml" - else: - disp_filename = f"phono3py_disp.{input_filename}.yaml" - if phonon_smat is not None: - forces_filename = "FORCES_FC2" - else: - forces_filename = "FORCES_FC3" - - if log_level: - print(f'Displacement dataset is read from "{disp_filename}".') - print(f'Forces are read from "{forces_filename}"') - - with open(forces_filename, "r") as f: - len_first_line = get_length_of_first_line(f) - - if len_first_line == 3: - file_exists(disp_filename, log_level) - file_exists(forces_filename, log_level) - ph3yml = Phono3pyYaml() - ph3yml.read(disp_filename) - if phonon_smat is None: - dataset = copy.deepcopy(ph3yml.dataset) - smat = ph3yml.supercell_matrix - else: - dataset = copy.deepcopy(ph3yml.phonon_dataset) - smat = ph3yml.phonon_supercell_matrix - - if smat is None or (phonon_smat is not None and (phonon_smat != smat).any()): - if log_level: - print("") - print("Supercell matrix is inconsistent.") - print(f'Supercell matrix read from "{disp_filename}":') - print(smat) - print("Supercell matrix given by --dim-fc2:") - print(phonon_smat) - print_error() - sys.exit(1) - - parse_FORCES_FC2(dataset, filename=forces_filename) - write_FORCE_SETS(dataset) - - if log_level: - print("FORCE_SETS has been created.") - print_end_phono3py() - else: - if log_level: - print( - "The file format of %s is already readable by phonopy." - % forces_filename - ) - print_end_phono3py() - sys.exit(0) - - -def create_FORCES_FC3_and_FORCES_FC2_then_exit( - settings, - cell_filename: Optional[str], - log_level: Union[bool, int], -): - """Create FORCES_FC3 and FORCES_FC2 from files.""" - interface_mode = settings.calculator - ph3py_yaml = None - - ##################### - # 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 - - 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)) - - if not check_number_of_force_files(num_disps, force_filenames, disp_filename): - force_sets = [] - else: - force_sets = get_force_sets( - interface_mode, - num_atoms, - force_filenames, - verbose=(log_level > 0), - ) - - if settings.subtract_forces: - force_filename = settings.subtract_forces - file_exists(force_filename, log_level) - force_set_zero = get_force_sets( - interface_mode, - num_atoms, - [ - force_filename, - ], - verbose=(log_level > 0), - )[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_end_phono3py() - sys.exit(0) - else: - 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)) - force_sets = get_force_sets( - interface_mode, - num_atoms, - force_filenames, - verbose=(log_level > 0), - ) - - if settings.subtract_forces: - force_filename = settings.subtract_forces - file_exists(force_filename, log_level) - force_set_zero = get_force_sets( - interface_mode, - num_atoms, - [ - force_filename, - ], - verbose=(log_level > 0), - )[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 log_level: - print("") - print("%s has been created." % "FORCES_FC2") - print_end_phono3py() - sys.exit(0) - else: - if log_level: - print("") - print("%s could not be created." % "FORCES_FC2") - print_error() - sys.exit(1) - - def get_input_output_filenames_from_args(args): """Return strings inserted to input and output filenames.""" input_filename = args.input_filename @@ -1104,11 +849,17 @@ def main(**argparse_control): ) if args.force_sets_to_forces_fc2_mode: - create_FORCES_FC2_from_FORCE_SETS_then_exit(log_level) + create_FORCES_FC2_from_FORCE_SETS(log_level) + if log_level: + print_end_phono3py() + sys.exit(0) if args.force_sets_mode: - create_FORCE_SETS_from_FORCES_FCx_then_exit( + create_FORCE_SETS_from_FORCES_FCx( settings.phonon_supercell_matrix, input_filename, cell_filename, log_level ) + if log_level: + print_end_phono3py() + sys.exit(0) if args.write_grid_points: run_mode = "write_grid_info" elif args.show_num_triplets: @@ -1123,7 +874,15 @@ def main(**argparse_control): #################################### # Create FORCES_FC3 and FORCES_FC2 # #################################### - create_FORCES_FC3_and_FORCES_FC2_then_exit(settings, cell_filename, log_level) + if ( + settings.create_forces_fc3 + or settings.create_forces_fc3_file + or settings.create_forces_fc2 + ): + create_FORCES_FC3_and_FORCES_FC2(settings, cell_filename, log_level=log_level) + if log_level: + print_end_phono3py() + sys.exit(0) ########################################################### # Symmetry tolerance. Distance unit depends on interface. # @@ -1149,7 +908,6 @@ def main(**argparse_control): cell_info, settings, symprec, - output_filename=output_filename, interface_mode=interface_mode, log_level=log_level, ) diff --git a/phono3py/interface/phono3py_yaml.py b/phono3py/interface/phono3py_yaml.py index 79c03256..dc86c5d2 100644 --- a/phono3py/interface/phono3py_yaml.py +++ b/phono3py/interface/phono3py_yaml.py @@ -133,14 +133,11 @@ def _parse_dataset(self): This method override PhonopyYaml._parse_dataset. - key="phonon_displacements" at v2.2 or later. - key="displacements" at older version than v2.2. - """ self._data.phonon_dataset = self._get_dataset( - self._data.phonon_supercell, key="phonon_displacements" + self._data.phonon_supercell, key_prefix="phonon_" ) - if self._data.phonon_dataset is None: # key="displacements" + if self._data.phonon_dataset is None: self._data.phonon_dataset = self._get_dataset(self._data.phonon_supercell) def _parse_fc3_dataset(self): @@ -192,6 +189,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 "supercell_energy" in d1: + data1["supercell_energy"] = d1["supercell_energy"] d2_list = d1.get("paired_with") if d2_list is None: # backward compatibility d2_list = d1.get("second_atoms") @@ -226,6 +225,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 "supercell_energy" in d2: + second_atom_dict["supercell_energy"] = d2["supercell_energy"] disp2_id += 1 data1["second_atoms"].append(second_atom_dict) return disp2_id @@ -294,15 +295,6 @@ def _cell_info_yaml_lines(self): lines += self._phonon_supercell_yaml_lines() return lines - def _phonon_supercell_matrix_yaml_lines(self): - lines = [] - if self._data.phonon_supercell_matrix is not None: - lines.append("phonon_supercell_matrix:") - for v in self._data.supercell_matrix: - lines.append("- [ %3d, %3d, %3d ]" % tuple(v)) - lines.append("") - return lines - def _phonon_supercell_yaml_lines(self): lines = [] if self._data.phonon_supercell is not None: @@ -326,7 +318,7 @@ def _nac_yaml_lines(self): else: return self._nac_yaml_lines_given_symbols(self._data.primitive.symbols) - def _displacements_yaml_lines(self, with_forces=False, key="displacements"): + def _displacements_yaml_lines(self, with_forces: bool = False) -> list: """Get YAML lines for phonon_dataset and dataset. This method override PhonopyYaml._displacements_yaml_lines. @@ -337,16 +329,18 @@ def _displacements_yaml_lines(self, with_forces=False, key="displacements"): lines = [] if self._data.phonon_supercell_matrix is not None: lines += self._displacements_yaml_lines_2types( - self._data.phonon_dataset, with_forces=with_forces, key=f"phonon_{key}" + self._data.phonon_dataset, + with_forces=with_forces, + key_prefix="phonon_", ) lines += self._displacements_yaml_lines_2types( - self._data.dataset, with_forces=with_forces, key=key + self._data.dataset, with_forces=with_forces ) return lines def _displacements_yaml_lines_type1( - self, dataset, with_forces=False, key="displacements" - ): + self, dataset: dict, with_forces: bool = False, key_prefix: str = "" + ) -> list: """Get YAML lines for type1 phonon_dataset and dataset. This method override PhonopyYaml._displacements_yaml_lines_type1. @@ -354,7 +348,9 @@ def _displacements_yaml_lines_type1( Phono3pyYaml._displacements_yaml_lines_type1. """ - return displacements_yaml_lines_type1(dataset, with_forces=with_forces, key=key) + return displacements_yaml_lines_type1( + dataset, with_forces=with_forces, key_prefix=key_prefix + ) class Phono3pyYaml(PhonopyYaml): @@ -441,7 +437,9 @@ def set_phonon_info(self, phono3py: "Phono3py"): self._data.phonon_supercell = phono3py.phonon_supercell -def displacements_yaml_lines_type1(dataset, with_forces=False, key="displacements"): +def displacements_yaml_lines_type1( + dataset: dict, with_forces: bool = False, key_prefix: str = "" +) -> list: """Get YAML lines for type1 phonon_dataset and dataset. This is a function but not class method because used by other function. @@ -452,7 +450,7 @@ def displacements_yaml_lines_type1(dataset, with_forces=False, key="displacement if "second_atoms" in dataset["first_atoms"][0]: lines = ["displacement_pairs:"] else: - lines = [f"{key}:"] + lines = [f"{key_prefix}displacements:"] for disp1_id, d1 in enumerate(dataset["first_atoms"]): lines.append("- atom: %4d" % (d1["number"] + 1)) lines.append(" displacement:") @@ -464,6 +462,12 @@ 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 "supercell_energy" in d1: + lines.append( + " supercell_energy: {energy:.8f}".format( + energy=d1["supercell_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 +546,12 @@ 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 "supercell_energy" in dataset2[j]: + lines.append( + " supercell_energy: {energy:.8f}".format( + energy=dataset2[j]["supercell_energy"] + ) + ) else: lines.append(" - atom: %4d" % (i + 1)) lines.append( diff --git a/phono3py/version.py b/phono3py/version.py index c2305594..3426bb80 100644 --- a/phono3py/version.py +++ b/phono3py/version.py @@ -34,4 +34,4 @@ # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. -__version__ = "3.0.2" +__version__ = "3.0.3" diff --git a/requirements.txt b/requirements.txt index ed682251..a8ad2963 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,4 +2,4 @@ numpy >= 1.17.0 PyYAML >= 5.3 matplotlib >= 2.2.2 h5py >= 3.0 -phonopy >=2.22,<2.23 +phonopy >=2.23,<2.24 diff --git a/setup.py b/setup.py index 2f0cbba5..26c05a49 100644 --- a/setup.py +++ b/setup.py @@ -326,7 +326,7 @@ def main(build_dir): "matplotlib>=2.2.2", "h5py>=3.0", "spglib>=2.0", - "phonopy>=2.22,<2.23", + "phonopy>=2.23,<2.24", ], provides=["phono3py"], scripts=scripts_phono3py, diff --git a/test/Si-111-222-fd.tar.xz b/test/Si-111-222-fd.tar.xz new file mode 100644 index 00000000..9c395230 Binary files /dev/null and b/test/Si-111-222-fd.tar.xz differ diff --git a/test/Si-111-222-rd.tar.xz b/test/Si-111-222-rd.tar.xz new file mode 100644 index 00000000..871c3e46 Binary files /dev/null and b/test/Si-111-222-rd.tar.xz differ diff --git a/test/api/test_api_phono3py.py b/test/api/test_api_phono3py.py index 38687170..cbed29eb 100644 --- a/test/api/test_api_phono3py.py +++ b/test/api/test_api_phono3py.py @@ -1,7 +1,11 @@ """Tests of Phono3py API.""" +from __future__ import annotations + from pathlib import Path +import numpy as np + from phono3py import Phono3py cwd = Path(__file__).parent @@ -40,3 +44,125 @@ def test_displacements_setter_Si(si_pbesol_111_222_fd: Phono3py): ) ph3.displacements = displacements ph3.phonon_displacements = phonon_displacements + + +def test_type1_forces_energies_setter_Si(si_111_222_fd: Phono3py): + """Test type1 supercell_energies, phonon_supercell_energies attributes.""" + ph3_in = si_111_222_fd + ref_ph_supercell_energies = [-346.85204143] + ref_supercell_energies = [ + -43.3509760, + -43.33608775, + -43.35352904, + -43.34370672, + -43.34590849, + -43.34540162, + -43.34421408, + -43.34481089, + -43.34703607, + -43.34241924, + -43.34786243, + -43.34168203, + -43.34274245, + -43.34703607, + -43.34786243, + -43.34184454, + ] + np.testing.assert_allclose(ph3_in.supercell_energies, ref_supercell_energies) + np.testing.assert_allclose( + ph3_in.phonon_supercell_energies, ref_ph_supercell_energies + ) + + ref_force00 = [-0.4109520800000000, 0.0000000100000000, 0.0000000300000000] + ref_force_last = [0.1521426300000000, 0.0715600600000000, -0.0715600700000000] + ref_ph_force00 = [-0.4027479600000000, 0.0000000200000000, 0.0000001000000000] + np.testing.assert_allclose(ph3_in.forces[0, 0], ref_force00) + np.testing.assert_allclose(ph3_in.forces[-1, -1], ref_force_last) + np.testing.assert_allclose(ph3_in.phonon_forces[0, 0], ref_ph_force00) + + ph3 = Phono3py( + ph3_in.unitcell, + supercell_matrix=ph3_in.supercell_matrix, + phonon_supercell_matrix=ph3_in.phonon_supercell_matrix, + primitive_matrix=ph3_in.primitive_matrix, + ) + ph3.dataset = ph3_in.dataset + ph3.phonon_dataset = ph3_in.phonon_dataset + + ph3.supercell_energies = ph3_in.supercell_energies + 1 + ph3.phonon_supercell_energies = ph3_in.phonon_supercell_energies + 1 + np.testing.assert_allclose(ph3_in.supercell_energies + 1, ph3.supercell_energies) + np.testing.assert_allclose( + ph3_in.phonon_supercell_energies + 1, ph3.phonon_supercell_energies + ) + + ph3.forces = ph3_in.forces + 1 + ph3.phonon_forces = ph3_in.phonon_forces + 1 + np.testing.assert_allclose(ph3_in.forces + 1, ph3.forces) + np.testing.assert_allclose(ph3_in.phonon_forces + 1, ph3.phonon_forces) + + +def test_type2_forces_energies_setter_Si(si_111_222_rd: Phono3py): + """Test type2 supercell_energies, phonon_supercell_energies attributes.""" + ph3_in = si_111_222_rd + ref_ph_supercell_energies = [ + -346.81061270, # 1 + -346.81263617, # 2 + ] + ref_supercell_energies = [ + -43.35270268, # 1 + -43.35211687, # 2 + -43.35122776, # 3 + -43.35226673, # 4 + -43.35146358, # 5 + -43.35133209, # 6 + -43.35042212, # 7 + -43.35008442, # 8 + -43.34968796, # 9 + -43.35348999, # 10 + -43.35134937, # 11 + -43.35335251, # 12 + -43.35160892, # 13 + -43.35009115, # 14 + -43.35202797, # 15 + -43.35076370, # 16 + -43.35174477, # 17 + -43.35107001, # 18 + -43.35037949, # 19 + -43.35126123, # 20 + ] + np.testing.assert_allclose(ph3_in.supercell_energies, ref_supercell_energies) + np.testing.assert_allclose( + ph3_in.phonon_supercell_energies, ref_ph_supercell_energies + ) + + ref_force00 = [0.0445647800000000, 0.1702929900000000, 0.0913398200000000] + ref_force_last = [-0.1749668700000000, 0.0146997300000000, -0.1336066300000000] + ref_ph_force00 = [-0.0161598900000000, -0.1161657500000000, 0.1399128100000000] + ref_ph_force_last = [0.1049486700000000, 0.0795870900000000, 0.1062164600000000] + + np.testing.assert_allclose(ph3_in.forces[0, 0], ref_force00) + np.testing.assert_allclose(ph3_in.forces[-1, -1], ref_force_last) + np.testing.assert_allclose(ph3_in.phonon_forces[0, 0], ref_ph_force00) + np.testing.assert_allclose(ph3_in.phonon_forces[-1, -1], ref_ph_force_last) + + ph3 = Phono3py( + ph3_in.unitcell, + supercell_matrix=ph3_in.supercell_matrix, + phonon_supercell_matrix=ph3_in.phonon_supercell_matrix, + primitive_matrix=ph3_in.primitive_matrix, + ) + ph3.dataset = ph3_in.dataset + ph3.phonon_dataset = ph3_in.phonon_dataset + ph3.supercell_energies = ph3_in.supercell_energies + 1 + ph3.phonon_supercell_energies = ph3_in.phonon_supercell_energies + 1 + + np.testing.assert_allclose(ph3_in.supercell_energies + 1, ph3.supercell_energies) + np.testing.assert_allclose( + ph3_in.phonon_supercell_energies + 1, ph3.phonon_supercell_energies + ) + + ph3.forces = ph3_in.forces + 1 + ph3.phonon_forces = ph3_in.phonon_forces + 1 + np.testing.assert_allclose(ph3_in.forces + 1, ph3.forces) + np.testing.assert_allclose(ph3_in.phonon_forces + 1, ph3.phonon_forces) diff --git a/test/conftest.py b/test/conftest.py index bfb10096..9b7403f1 100644 --- a/test/conftest.py +++ b/test/conftest.py @@ -1,5 +1,6 @@ """Pytest conftest.py.""" +import tarfile from pathlib import Path import numpy as np @@ -482,6 +483,20 @@ def aln_lda(request) -> Phono3py: ) +@pytest.fixture(scope="session") +def si_111_222_fd() -> Phono3py: + """Return Phono3py class instance of Si-1x1x1-2x2x2 FD.""" + yaml_filename = cwd / "phono3py_params_Si-111-222-fd.yaml.xz" + return phono3py.load(yaml_filename, produce_fc=False, log_level=1) + + +@pytest.fixture(scope="session") +def si_111_222_rd() -> Phono3py: + """Return Phono3py class instance of Si-1x1x1-2x2x2 RD.""" + yaml_filename = cwd / "phono3py_params_Si-111-222-rd.yaml.xz" + return phono3py.load(yaml_filename, produce_fc=False, log_level=1) + + @pytest.fixture(scope="session") def ph_nacl() -> Phonopy: """Return Phonopy class instance of NaCl 2x2x2.""" @@ -508,3 +523,76 @@ def ph_si() -> Phonopy: log_level=1, produce_fc=True, ) + + +@pytest.fixture(scope="session") +def si_111_222_fd_raw_data() -> tarfile.TarFile: + """Return Si fc3 111 fc2 222 vasp inputs. + + tar.getnames() + ['Si-111-222-fd', + 'Si-111-222-fd/vasprun-001.xml', + 'Si-111-222-fd/vasprun-015.xml', + 'Si-111-222-fd/vasprun-014.xml', + 'Si-111-222-fd/vasprun-000.xml', + 'Si-111-222-fd/vasprun-016.xml', + 'Si-111-222-fd/vasprun-002.xml', + 'Si-111-222-fd/vasprun-003.xml', + 'Si-111-222-fd/vasprun-013.xml', + 'Si-111-222-fd/vasprun-007.xml', + 'Si-111-222-fd/vasprun-fc2-000.xml', + 'Si-111-222-fd/vasprun-fc2-001.xml', + 'Si-111-222-fd/vasprun-006.xml', + 'Si-111-222-fd/vasprun-012.xml', + 'Si-111-222-fd/vasprun-004.xml', + 'Si-111-222-fd/vasprun-010.xml', + 'Si-111-222-fd/vasprun-011.xml', + 'Si-111-222-fd/vasprun-005.xml', + 'Si-111-222-fd/phono3py_disp.yaml', + 'Si-111-222-fd/vasprun-008.xml', + 'Si-111-222-fd/vasprun-009.xml'] + member = tar.getmember("Si-111-222-fd/phono3py_disp.yaml") + tar.extractfile(member) -> byte file object + + """ + tar = tarfile.open(cwd / "Si-111-222-fd.tar.xz") + return tar + + +@pytest.fixture(scope="session") +def si_111_222_rd_raw_data() -> tarfile.TarFile: + """Return Si fc3 111 fc2 222 vasp inputs. + + tar.getnames() + ['Si-111-222-rd', + 'Si-111-222-rd/vasprun-001.xml', + 'Si-111-222-rd/vasprun-015.xml', + 'Si-111-222-rd/vasprun-014.xml', + 'Si-111-222-rd/vasprun-000.xml', + 'Si-111-222-rd/vasprun-016.xml', + 'Si-111-222-rd/vasprun-002.xml', + 'Si-111-222-rd/vasprun-003.xml', + 'Si-111-222-rd/vasprun-017.xml', + 'Si-111-222-rd/vasprun-013.xml', + 'Si-111-222-rd/vasprun-007.xml', + 'Si-111-222-rd/vasprun-fc2-000.xml', + 'Si-111-222-rd/vasprun-fc2-001.xml', + 'Si-111-222-rd/vasprun-006.xml', + 'Si-111-222-rd/vasprun-012.xml', + 'Si-111-222-rd/vasprun-004.xml', + 'Si-111-222-rd/vasprun-010.xml', + 'Si-111-222-rd/vasprun-fc2-002.xml', + 'Si-111-222-rd/vasprun-011.xml', + 'Si-111-222-rd/vasprun-005.xml', + 'Si-111-222-rd/phono3py_disp.yaml', + 'Si-111-222-rd/vasprun-020.xml', + 'Si-111-222-rd/vasprun-008.xml', + 'Si-111-222-rd/vasprun-009.xml', + 'Si-111-222-rd/vasprun-019.xml', + 'Si-111-222-rd/vasprun-018.xml'] + member = tar.getmember("Si-111-222-rd/phono3py_disp.yaml") + tar.extractfile(member) -> byte file object + + """ + tar = tarfile.open(cwd / "Si-111-222-fd.tar.xz") + return tar diff --git a/test/phono3py_params_Si-111-222-fd.yaml.xz b/test/phono3py_params_Si-111-222-fd.yaml.xz new file mode 100644 index 00000000..26b05e18 Binary files /dev/null and b/test/phono3py_params_Si-111-222-fd.yaml.xz differ diff --git a/test/phono3py_params_Si-111-222-rd.yaml.xz b/test/phono3py_params_Si-111-222-rd.yaml.xz new file mode 100644 index 00000000..626dfbc8 Binary files /dev/null and b/test/phono3py_params_Si-111-222-rd.yaml.xz differ