diff --git a/phono3py/api_phono3py.py b/phono3py/api_phono3py.py index 24dba4af..e758b027 100644 --- a/phono3py/api_phono3py.py +++ b/phono3py/api_phono3py.py @@ -40,10 +40,12 @@ from typing import Literal, Optional, Union import numpy as np +from phonopy import Phonopy from phonopy.exception import ForceCalculatorRequiredError from phonopy.harmonic.displacement import ( directions_to_displacement_dataset, get_least_displacements, + get_random_displacements_dataset, ) from phonopy.harmonic.dynamical_matrix import DynamicalMatrix from phonopy.harmonic.force_constants import get_fc2 as get_phonopy_fc2 @@ -54,6 +56,7 @@ symmetrize_force_constants, ) from phonopy.interface.fc_calculator import get_fc2 +from phonopy.interface.pypolymlp import PypolymlpParams from phonopy.structure.atoms import PhonopyAtoms from phonopy.structure.cells import ( Primitive, @@ -288,10 +291,13 @@ def __init__( self._frequency_points = None self._temperatures = None - # Other variables + # Force constants self._fc2 = None self._fc3 = None + # MLP container + self._mlp_phonopy = None + # Setup interaction self._interaction = None self._band_indices = None @@ -984,6 +990,30 @@ def phonon_supercell_energies(self): def phonon_supercell_energies(self, values): self._set_phonon_forces_energies(values, target="supercell_energies") + @property + def mlp_dataset(self) -> Optional[dict]: + """Return displacement-force dataset. + + The supercell matrix is equal to that of usual displacement-force + dataset. Only type 2 format is supported. "displacements", + "forces", and "supercell_energies" should be contained. + + """ + if self._mlp_phonopy is None: + return None + return self._mlp_phonopy.mlp_dataset + + @mlp_dataset.setter + def mlp_dataset(self, mlp_dataset: dict): + if self._mlp_phonopy is None: + self._mlp_phonopy = Phonopy( + self._unitcell, + supercell_matrix=self._supercell_matrix, + symprec=self._symprec, + log_level=self._log_level, + ) + self._mlp_phonopy.mlp_dataset = mlp_dataset + @property def phph_interaction(self): """Return Interaction instance.""" @@ -1166,6 +1196,10 @@ def generate_displacements( cutoff_pair_distance=None, is_plusminus="auto", is_diagonal=True, + number_of_snapshots: Optional[int] = None, + random_seed: Optional[int] = None, + is_random_distance: bool = False, + min_distance: Optional[float] = None, ): """Generate displacement dataset in supercell for fc3. @@ -1211,20 +1245,48 @@ def generate_displacements( vectors of the supercell. With True, direction not along the basis vectors can be chosen when the number of the displacements may be reduced. + number_of_snapshots : int or None, optional + Number of snapshots of supercells with random displacements. Random + displacements are generated displacing all atoms in random + directions with a fixed displacement distance specified by + 'distance' parameter, i.e., all atoms in supercell are displaced + with the same displacement distance in direct space. Default is + None. + random_seed : int or None, optional + Random seed for random displacements generation. Default is None. + is_random_distance : bool, optional + Random direction displacements are generated also with random + amplitudes. The maximum value is defined by `distance` and minimum + value is given by `min_distance`. Default is False. + min_distance : float or None, optional + In random direction displacements generation with random distance + (`is_random_distance=True`), the minimum distance is given by this + value. """ - direction_dataset = get_third_order_displacements( - self._supercell, - self._symmetry, - is_plusminus=is_plusminus, - is_diagonal=is_diagonal, - ) - self._dataset = direction_to_displacement( - direction_dataset, - distance, - self._supercell, - cutoff_distance=cutoff_pair_distance, - ) + if number_of_snapshots is not None and number_of_snapshots > 0: + self._dataset = self._generate_random_displacements( + number_of_snapshots, + len(self._supercell), + distance=distance, + is_plusminus=is_plusminus is True, + random_seed=random_seed, + is_random_distance=is_random_distance, + min_distance=min_distance, + ) + else: + direction_dataset = get_third_order_displacements( + self._supercell, + self._symmetry, + is_plusminus=is_plusminus, + is_diagonal=is_diagonal, + ) + self._dataset = direction_to_displacement( + direction_dataset, + distance, + self._supercell, + cutoff_distance=cutoff_pair_distance, + ) self._supercells_with_displacements = None if self._phonon_supercell_matrix is not None and self._phonon_dataset is None: @@ -1233,7 +1295,14 @@ def generate_displacements( ) def generate_fc2_displacements( - self, distance=0.03, is_plusminus="auto", is_diagonal=False + self, + distance: float = 0.03, + is_plusminus: str = "auto", + is_diagonal: bool = False, + number_of_snapshots: Optional[int] = None, + random_seed: Optional[int] = None, + is_random_distance: bool = False, + min_distance: Optional[float] = None, ): """Generate displacement dataset in phonon supercell for fc2. @@ -1263,6 +1332,23 @@ def generate_fc2_displacements( vectors of the supercell. With True, direction not along the basis vectors can be chosen when the number of the displacements may be reduced. Default is False. + number_of_snapshots : int or None, optional + Number of snapshots of supercells with random displacements. Random + displacements are generated displacing all atoms in random + directions with a fixed displacement distance specified by + 'distance' parameter, i.e., all atoms in supercell are displaced + with the same displacement distance in direct space. Default is + None. + random_seed : int or None, optional + Random seed for random displacements generation. Default is None. + is_random_distance : bool, optional + Random direction displacements are generated also with random + amplitudes. The maximum value is defined by `distance` and minimum + value is given by `min_distance`. Default is False. + min_distance : float or None, optional + In random direction displacements generation with random distance + (`is_random_distance=True`), the minimum distance is given by this + value. """ if self._phonon_supercell_matrix is None: @@ -1273,54 +1359,85 @@ def generate_fc2_displacements( ) raise RuntimeError(msg) - phonon_displacement_directions = get_least_displacements( - self._phonon_supercell_symmetry, - is_plusminus=is_plusminus, - is_diagonal=is_diagonal, - ) - self._phonon_dataset = directions_to_displacement_dataset( - phonon_displacement_directions, distance, self._phonon_supercell - ) + if number_of_snapshots is not None and number_of_snapshots > 0: + self._phonon_dataset = self._generate_random_displacements( + number_of_snapshots, + len(self._phonon_supercell), + distance=distance, + is_plusminus=is_plusminus is True, + random_seed=random_seed, + is_random_distance=is_random_distance, + min_distance=min_distance, + ) + else: + phonon_displacement_directions = get_least_displacements( + self._phonon_supercell_symmetry, + is_plusminus=is_plusminus, + is_diagonal=is_diagonal, + ) + self._phonon_dataset = directions_to_displacement_dataset( + phonon_displacement_directions, distance, self._phonon_supercell + ) self._phonon_supercells_with_displacements = None def produce_fc3( self, - symmetrize_fc3r=False, - is_compact_fc=False, - fc_calculator=None, - fc_calculator_options=None, + symmetrize_fc3r: bool = False, + is_compact_fc: bool = False, + fc_calculator: Optional[str] = None, + fc_calculator_options: Optional[str] = None, + use_pypolymlp: bool = False, + mlp_params: Optional[Union[PypolymlpParams, dict]] = None, ): """Calculate fc3 from displacements and forces. Parameters ---------- symmetrize_fc3r : bool - Only for type 1 displacement_dataset, translational and - permutation symmetries are applied after creating fc3. This - symmetrization is not very sophisticated and can break space - group symmetry, but often useful. If better symmetrization is - expected, it is recommended to use external force constants - calculator such as ALM. Default is False. + Only for type 1 displacement_dataset, translational and permutation + symmetries are applied after creating fc3. This symmetrization is + not very sophisticated and can break space group symmetry, but often + useful. If better symmetrization is expected, it is recommended to + use external force constants calculator such as ALM. Default is + False. is_compact_fc : bool fc3 shape is - False: (supercell, supercell, supecell, 3, 3, 3) - True: (primitive, supercell, supecell, 3, 3, 3) + False: (supercell, supercell, supecell, 3, 3, 3) True: + (primitive, supercell, supecell, 3, 3, 3) where 'supercell' and 'primitive' indicate number of atoms in these cells. Default is False. fc_calculator : str or None Force constants calculator given by str. fc_calculator_options : dict Options for external force constants calculator. + use_pypolymlp : bool + Use MLP of pypolymlp to calculate fc3. Default is False. mlp_dataset + has to be set before calling this method. + mlp_params : PypolymlpParams or dict, optional + Parameters for developing MLP. Default is None. When dict is given, + PypolymlpParams instance is created from the dict. """ - disp_dataset = self._dataset + if use_pypolymlp: + if self._mlp_phonopy is None: + msg = "mlp_dataset has to be set before calling this method." + raise RuntimeError(msg) + if self._dataset is None or "displacements" not in self._dataset: + raise RuntimeError( + "Type 2 displacements are not set. Run generate_displacements." + ) + self._mlp_phonopy.develop_mlp(params=mlp_params) + self._mlp_phonopy.displacements = self._dataset["displacements"] + self._mlp_phonopy.evaluate_mlp() + self._dataset["forces"] = self._mlp_phonopy.forces + self._dataset["supercell_energies"] = self._mlp_phonopy.supercell_energies fc3_calculator, fc3_calculator_options = self._extract_fc2_fc3_calculators( fc_calculator, fc_calculator_options, 3 ) if fc3_calculator is not None: - disps, forces = get_displacements_and_forces_fc3(disp_dataset) + disps, forces = get_displacements_and_forces_fc3(self._dataset) fc2, fc3 = get_fc3( self._supercell, self._primitive, @@ -1333,7 +1450,7 @@ def produce_fc3( log_level=self._log_level, ) else: - if "displacements" in disp_dataset: + if "displacements" in self._dataset: msg = ( "fc_calculator has to be set to produce force " "constans from this dataset." @@ -1342,7 +1459,7 @@ def produce_fc3( fc2, fc3 = get_phono3py_fc3( self._supercell, self._primitive, - disp_dataset, + self._dataset, self._symmetry, is_compact_fc=is_compact_fc, verbose=self._log_level, @@ -2404,3 +2521,31 @@ def _set_phonon_forces_energies( self._phonon_dataset[target] = _values else: raise RuntimeError("Set of FC2 displacements is not available.") + + def _generate_random_displacements( + self, + number_of_snapshots: int, + number_of_atoms: int, + distance: float = 0.03, + is_plusminus: bool = False, + random_seed: Optional[int] = None, + is_random_distance: bool = False, + min_distance: Optional[float] = None, + ): + if random_seed is not None and random_seed >= 0 and random_seed < 2**32: + _random_seed = random_seed + dataset = {"random_seed": _random_seed} + else: + _random_seed = None + dataset = {} + d = get_random_displacements_dataset( + number_of_snapshots, + number_of_atoms, + distance, + random_seed=_random_seed, + is_plusminus=is_plusminus, + is_random_distance=is_random_distance, + min_distance=min_distance, + ) + dataset["displacements"] = d + return dataset diff --git a/phono3py/cui/phono3py_argparse.py b/phono3py/cui/phono3py_argparse.py index 78e734e5..48bf9820 100644 --- a/phono3py/cui/phono3py_argparse.py +++ b/phono3py/cui/phono3py_argparse.py @@ -618,6 +618,13 @@ def get_parser(fc_symmetry=False, is_nac=False, load_phono3py_yaml=False): default=None, help="Conversion factor for ph-ph interaction", ) + parser.add_argument( + "--pypolymlp", + dest="use_pypolymlp", + action="store_true", + default=None, + help="Use pypolymlp and symfc for generating force constants", + ) parser.add_argument( "--qpoints", nargs="+", diff --git a/test/api/test_api_phono3py.py b/test/api/test_api_phono3py.py index cbed29eb..35a1b2fb 100644 --- a/test/api/test_api_phono3py.py +++ b/test/api/test_api_phono3py.py @@ -5,6 +5,8 @@ from pathlib import Path import numpy as np +import pytest +from phonopy.interface.pypolymlp import PypolymlpParams from phono3py import Phono3py @@ -166,3 +168,37 @@ def test_type2_forces_energies_setter_Si(si_111_222_rd: Phono3py): 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_use_pypolymlp_mgs(mgo_222rd_444rd: Phono3py): + """Test use_pypolymlp in produce_fc3.""" + pytest.importorskip("pypolymlp") + + ph3_in = mgo_222rd_444rd + 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, + log_level=2, + ) + ph3.mlp_dataset = { + "displacements": ph3_in.displacements[:10], + "forces": ph3_in.forces[:10], + "supercell_energies": ph3_in.supercell_energies[:10], + } + ph3.phonon_dataset = ph3_in.phonon_dataset + ph3.nac_params = ph3_in.nac_params + ph3.displacements = ph3_in.displacements[:100] + + atom_energies = {"Mg": -0.00896717, "O": -0.95743902} + params = PypolymlpParams(gtinv_maxl=(4, 4), atom_energies=atom_energies) + ph3.produce_fc3(use_pypolymlp=True, mlp_params=params, fc_calculator="symfc") + ph3.produce_fc2(fc_calculator="symfc") + + ph3.mesh_numbers = 30 + ph3.init_phph_interaction() + ph3.run_thermal_conductivity(temperatures=[300]) + assert ( + pytest.approx(63.6001137, abs=1e-3) == ph3.thermal_conductivity.kappa[0, 0, 0] + ) diff --git a/test/conftest.py b/test/conftest.py index 5d552de1..91fda478 100644 --- a/test/conftest.py +++ b/test/conftest.py @@ -497,6 +497,17 @@ def si_111_222_rd() -> Phono3py: return phono3py.load(yaml_filename, produce_fc=False, log_level=1) +@pytest.fixture(scope="session") +def mgo_222rd_444rd() -> Phono3py: + """Return Phono3py class instance of MgO-2x2x2-4x4x4 RD-RD. + + 4 and 400 supercells for fc2 and fc3, respectively. + + """ + yaml_filename = cwd / "phono3py_params_MgO-222rd-444fd.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.""" diff --git a/test/phono3py_params_MgO-222rd-444fd.yaml.xz b/test/phono3py_params_MgO-222rd-444fd.yaml.xz new file mode 100644 index 00000000..59148cb0 Binary files /dev/null and b/test/phono3py_params_MgO-222rd-444fd.yaml.xz differ