diff --git a/phono3py/api_phono3py.py b/phono3py/api_phono3py.py index 1cb98b7f..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 @@ -610,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 @@ -643,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): @@ -822,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): @@ -928,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 @@ -942,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): @@ -2280,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/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/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