From e57e1c56bfd7f41f102d4cd1f60ebe553357813f Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Sat, 29 Jun 2024 16:33:15 +0900 Subject: [PATCH 1/2] Update pypolymlp feature --- phono3py/api_phono3py.py | 218 ++++++++++++++++++------- phono3py/cui/create_force_constants.py | 154 +++++++++-------- phono3py/cui/load.py | 54 +++--- phono3py/cui/phono3py_script.py | 9 +- phono3py/interface/phono3py_yaml.py | 2 +- 5 files changed, 272 insertions(+), 165 deletions(-) diff --git a/phono3py/api_phono3py.py b/phono3py/api_phono3py.py index c6190549..42b83762 100644 --- a/phono3py/api_phono3py.py +++ b/phono3py/api_phono3py.py @@ -303,6 +303,8 @@ def __init__( # MLP self._mlp = None self._mlp_dataset = None + self._phonon_mlp = None + self._phonon_mlp_dataset = None # Setup interaction self._interaction = None @@ -700,22 +702,25 @@ def mlp_dataset(self) -> Optional[dict]: @mlp_dataset.setter def mlp_dataset(self, mlp_dataset: dict): - if not isinstance(mlp_dataset, dict): - raise TypeError("mlp_dataset has to be a dictionary.") - if "displacements" not in mlp_dataset: - raise RuntimeError("Displacements have to be given.") - if "forces" not in mlp_dataset: - raise RuntimeError("Forces have to be given.") - if "supercell_energy" in mlp_dataset: - raise RuntimeError("Supercell energies have to be given.") - if len(mlp_dataset["displacements"]) != len(mlp_dataset["forces"]): - raise RuntimeError("Length of displacements and forces are different.") - if len(mlp_dataset["displacements"]) != len(mlp_dataset["supercell_energies"]): - raise RuntimeError( - "Length of displacements and supercell_energies are different." - ) + self._check_mlp_dataset(mlp_dataset) self._mlp_dataset = mlp_dataset + @property + def phonon_mlp_dataset(self) -> Optional[dict]: + """Return phonon displacement-force dataset. + + The phonon 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. + + """ + return self._phonon_mlp_dataset + + @phonon_mlp_dataset.setter + def phonon_mlp_dataset(self, mlp_dataset: dict): + self._check_mlp_dataset(mlp_dataset) + self._phonon_mlp_dataset = mlp_dataset + @property def band_indices(self) -> list[np.ndarray]: """Setter and getter of band indices. @@ -950,18 +955,18 @@ def phonon_displacements(self): shape=(supercells, natom, 3), dtype='double', order='C' """ - if self._phonon_supercell_matrix is None: - raise RuntimeError("phonon_supercell_matrix is not set.") - - dataset = self._phonon_dataset - if "first_atoms" in dataset: - num_scells = len(dataset["first_atoms"]) + if self._phonon_dataset is None: + raise RuntimeError("phonon_dataset is not set.") + if "first_atoms" in self._phonon_dataset: + num_scells = len(self._phonon_dataset["first_atoms"]) natom = len(self._phonon_supercell) displacements = np.zeros((num_scells, natom, 3), dtype="double", order="C") - for i, disp1 in enumerate(dataset["first_atoms"]): + for i, disp1 in enumerate(self._phonon_dataset["first_atoms"]): displacements[i, disp1["number"]] = disp1["displacement"] - elif "forces" in dataset or "displacements" in dataset: - displacements = dataset["displacements"] + elif ( + "forces" in self._phonon_dataset or "displacements" in self._phonon_dataset + ): + displacements = self._phonon_dataset["displacements"] else: raise RuntimeError("displacement dataset has wrong format.") @@ -1045,6 +1050,16 @@ def grid(self): """ return self._bz_grid + @property + def mlp(self): + """Return MLP instance.""" + return self._mlp + + @property + def phonon_mlp(self): + """Return MLP instance for fc2.""" + return self._phonon_mlp + def init_phph_interaction( self, nac_q_direction=None, @@ -1214,48 +1229,56 @@ def generate_displacements( ): """Generate displacement dataset in supercell for fc3. - This systematically generates single and pair atomic displacements - in supercells to calculate fc3 considering crystal symmetry. - When this method is called, existing cache of supercells with - displacements for fc3 are removed. - - For fc3, two atoms are displaced for each configuration - considering crystal symmetry. The first displacement is chosen - in the perfect supercell, and the second displacement in the - displaced supercell. The first displacements are taken along - the basis vectors of the supercell. This is because the - symmetry is expected to be less broken by the introduced first - displacement, and as the result, the number of second - displacements may become smaller than the case that the first - atom is displaced not along the basis vectors. + Systematic displacements + ------------------------ + + Unless number_of_snapshots is specified, this method systematically + generates single and pair atomic displacements in supercells to + calculate fc3 considering crystal symmetry. + + For fc3, two atoms are displaced for each configuration considering + crystal symmetry. The first displacement is chosen in the perfect + supercell, and the second displacement in the displaced supercell. The + first displacements are taken along the basis vectors of the supercell. + This is because the symmetry is expected to be less broken by the + introduced first displacement, and as the result, the number of second + displacements may become smaller than the case that the first atom is + displaced not along the basis vectors. + + Random displacements + -------------------- + Unless number_of_snapshots is specified, displacements are generated + randomly. There are several options how the random displacements are + generated. Note ---- - When phonon_supercell_matrix is not given, fc2 is also - computed from the same set of the displacements for fc3 and - respective supercell forces. When phonon_supercell_matrix is - set, the displacements in phonon_supercell are generated unless - those already exist. + When phonon_supercell_matrix is not given, fc2 is also computed from the + same set of the displacements for fc3 and respective supercell forces. + When phonon_supercell_matrix is set, the displacements in + phonon_supercell are generated unless those already exist. If a specific + set of displacements for fc2 is expected, generate_fc2_displacements + should be called. Parameters ---------- distance : float, optional Constant displacement Euclidean distance. Default is 0.03. cutoff_pair_distance : float, optional - This is used as a cutoff Euclidean distance to determine if - each pair of displacements is considered to calculate fc3 or not. - Default is None, which means cutoff is not used. + This is used as a cutoff Euclidean distance to determine if each + pair of displacements is considered to calculate fc3 or not. Default + is None, which means cutoff is not used. is_plusminus : True, False, or 'auto', optional With True, atomis are displaced in both positive and negative - directions. With False, only one direction. With 'auto', - mostly equivalent to is_plusminus=True, but only one direction - is chosen when the displacements in both directions are - symmetrically equivalent. Default is 'auto'. + directions. With False, only one direction. With 'auto', mostly + equivalent to is_plusminus=True, but only one direction is chosen + when the displacements in both directions are symmetrically + equivalent. Default is 'auto'. is_diagonal : Bool, optional With False, the second displacements are made along the basis vectors of the supercell. With True, direction not along the basis - vectors can be chosen when the number of the displacements - may be reduced. + 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 @@ -1362,14 +1385,6 @@ def generate_fc2_displacements( value. """ - if self._phonon_supercell_matrix is None: - msg = ( - "phonon_supercell_matrix is not set. " - "This method is used to generate displacements to " - "calculate phonon_fc2." - ) - raise RuntimeError(msg) - if number_of_snapshots is not None and number_of_snapshots > 0: self._phonon_dataset = self._generate_random_displacements( number_of_snapshots, @@ -2214,6 +2229,75 @@ def evaluate_mlp(self): self.supercell_energies = energies self.forces = forces + def develop_phonon_mlp( + self, params: Optional[Union[PypolymlpParams, dict, str]] = None + ): + """Develop MLP of pypolymlp for fc2. + + Parameters + ---------- + params : PypolymlpParams or dict, optional + Parameters for developing MLP. Default is None. When dict is given, + PypolymlpParams instance is created from the dict. + + """ + if self._phonon_mlp_dataset is None: + raise RuntimeError("MLP dataset is not set.") + + if params is not None: + _params = parse_mlp_params(params) + else: + _params = params + + disps = self._phonon_mlp_dataset["displacements"] + forces = self._phonon_mlp_dataset["forces"] + energies = self._phonon_mlp_dataset["supercell_energies"] + n = int(len(disps) * 0.9) + train_data = PypolymlpData( + displacements=disps[:n], forces=forces[:n], supercell_energies=energies[:n] + ) + test_data = PypolymlpData( + displacements=disps[n:], forces=forces[n:], supercell_energies=energies[n:] + ) + self._phonon_mlp = develop_polymlp( + self._phonon_supercell, + train_data, + test_data, + params=_params, + verbose=self._log_level - 1 > 0, + ) + + def evaluate_phonon_mlp(self): + """Evaluate the machine learning potential of pypolymlp. + + This method calculates the supercell energies and forces from the MLP + for the displacements in self._dataset of type 2. The results are stored + in self._dataset. + + The displacements may be generated by the produce_force_constants method + with number_of_snapshots > 0. With MLP, a small distance parameter, such + as 0.001, can be numerically stable for the computation of force + constants. + + """ + if self._mlp is None and self._phonon_mlp is None: + raise RuntimeError("MLP is not developed yet.") + + if self.phonon_supercells_with_displacements is None: + raise RuntimeError( + "Displacements are not set. Run generate_fc2_displacements." + ) + + if self._phonon_mlp is None: + mlp = self._mlp + else: + mlp = self._phonon_mlp + energies, forces, _ = evalulate_polymlp( + mlp, self.phonon_supercells_with_displacements + ) + self.phonon_supercell_energies = energies + self.phonon_forces = forces + ################### # private methods # ################### @@ -2569,3 +2653,19 @@ def _generate_random_displacements( ) dataset["displacements"] = d return dataset + + def _check_mlp_dataset(self, mlp_dataset: dict): + if not isinstance(mlp_dataset, dict): + raise TypeError("mlp_dataset has to be a dictionary.") + if "displacements" not in mlp_dataset: + raise RuntimeError("Displacements have to be given.") + if "forces" not in mlp_dataset: + raise RuntimeError("Forces have to be given.") + if "supercell_energy" in mlp_dataset: + raise RuntimeError("Supercell energies have to be given.") + if len(mlp_dataset["displacements"]) != len(mlp_dataset["forces"]): + raise RuntimeError("Length of displacements and forces are different.") + if len(mlp_dataset["displacements"]) != len(mlp_dataset["supercell_energies"]): + raise RuntimeError( + "Length of displacements and supercell_energies are different." + ) diff --git a/phono3py/cui/create_force_constants.py b/phono3py/cui/create_force_constants.py index ee40bd8c..9f6f778a 100644 --- a/phono3py/cui/create_force_constants.py +++ b/phono3py/cui/create_force_constants.py @@ -167,29 +167,19 @@ def create_phono3py_force_constants( _read_phono3py_fc2(phono3py, symmetrize_fc2, input_filename, log_level) else: if phono3py.phonon_supercell_matrix is None: - if fc_calculator == "alm" and phono3py.fc2 is not None: - if log_level: - print("fc2 that was fit simultaneously with fc3 " "by ALM is used.") - else: - _create_phono3py_fc2( - phono3py, - ph3py_yaml, - symmetrize_fc2, - settings.is_compact_fc, - fc_calculator, - fc_calculator_options, - log_level, - ) + force_filename = "FORCES_FC3" else: - _create_phono3py_phonon_fc2( - phono3py, - ph3py_yaml, - symmetrize_fc2, - settings.is_compact_fc, - fc_calculator, - fc_calculator_options, - log_level, - ) + force_filename = "FORCES_FC2" + _create_phono3py_fc2( + phono3py, + ph3py_yaml, + force_filename, + symmetrize_fc2, + settings.is_compact_fc, + fc_calculator, + fc_calculator_options, + log_level, + ) if output_filename is None: filename = "fc2.hdf5" else: @@ -500,7 +490,7 @@ def _create_phono3py_fc3( phono3py.mlp_dataset = dataset run_pypolymlp_to_compute_forces( phono3py, - mlp_params=mlp_params, + mlp_params, displacement_distance=displacement_distance, number_of_snapshots=number_of_snapshots, random_seed=random_seed, @@ -519,10 +509,10 @@ def _create_phono3py_fc3( def run_pypolymlp_to_compute_forces( ph3py: Phono3py, mlp_params: Union[str, dict, PypolymlpParams], - displacement_distance: float, - number_of_snapshots: int, - random_seed: int, - log_level: int, + displacement_distance: Optional[float] = None, + number_of_snapshots: Optional[int] = None, + random_seed: Optional[int] = None, + log_level: int = 0, ): """Run pypolymlp to compute forces.""" if log_level: @@ -542,13 +532,14 @@ def run_pypolymlp_to_compute_forces( ph3py.develop_mlp(params=mlp_params) + if log_level: + print("-" * 30 + " pypolymlp end " + "-" * 31, flush=True) + if displacement_distance is None: _displacement_distance = 0.001 else: _displacement_distance = displacement_distance - if log_level > 1: - print("") if log_level: if number_of_snapshots: print("Generate random displacements") @@ -572,10 +563,10 @@ def run_pypolymlp_to_compute_forces( if log_level: print( - " Number of supercells for computing forces: " - f"{ph3py.displacements.shape[0]}" + f"Evaluate forces in {ph3py.displacements.shape[0]} supercells " + "by pypolymlp", + flush=True, ) - print("Evaluate forces by pypolymlp", flush=True) if ph3py.mlp_dataset is None: msg = "mlp_dataset has to be set before calling this method." @@ -585,26 +576,85 @@ def run_pypolymlp_to_compute_forces( ph3py.evaluate_mlp() + +def run_pypolymlp_to_compute_phonon_forces( + ph3py: Phono3py, + mlp_params: Optional[Union[str, dict, PypolymlpParams]] = None, + displacement_distance: Optional[float] = None, + number_of_snapshots: Optional[int] = None, + random_seed: Optional[int] = None, + log_level: int = 0, +): + """Run pypolymlp to compute phonon forces.""" + if ph3py.phonon_mlp_dataset is not None: + if log_level: + print("-" * 29 + " pypolymlp start " + "-" * 30) + print("Pypolymlp is a generator of polynomial machine learning potentials.") + print("Please cite the paper: A. Seko, J. Appl. Phys. 133, 011101 (2023).") + print("Pypolymlp is developed at https://github.com/sekocha/pypolymlp.") + if mlp_params: + print("Parameters:") + for k, v in asdict(parse_mlp_params(mlp_params)).items(): + if v is not None: + print(f" {k}: {v}") + if log_level > 1: + print("") + if log_level: + print("Developing MLPs by pypolymlp...", flush=True) + + ph3py.develop_phonon_mlp(params=mlp_params) + + if log_level: + print("-" * 30 + " pypolymlp end " + "-" * 31, flush=True) + + if displacement_distance is None: + _displacement_distance = 0.001 + else: + _displacement_distance = displacement_distance if log_level: - print("-" * 30 + " pypolymlp end " + "-" * 31, flush=True) + print("Generate random displacements for fc2") + print( + f" Displacement distance: {_displacement_distance:.5f}".rstrip("0").rstrip( + "." + ) + ) + ph3py.generate_fc2_displacements( + distance=_displacement_distance, + is_plusminus=True, + number_of_snapshots=number_of_snapshots, + random_seed=random_seed, + ) + if log_level: + print( + f"Evaluate forces in {ph3py.phonon_displacements.shape[0]} " + "supercells by pypolymlp", + flush=True, + ) + ph3py.evaluate_phonon_mlp() def _create_phono3py_fc2( phono3py: Phono3py, ph3py_yaml: Optional[Phono3pyYaml], + force_filename, symmetrize_fc2, is_compact_fc, fc_calculator, fc_calculator_options, log_level, ): + """Read forces and produce fc2. + + force_filename is either "FORCES_FC2" or "FORCES_FC3". + + """ _ph3py_yaml = _get_default_ph3py_yaml(ph3py_yaml) try: dataset = parse_forces( phono3py, ph3py_yaml=_ph3py_yaml, - force_filename="FORCES_FC3", + force_filename=force_filename, fc_type="fc2", log_level=log_level, ) @@ -633,42 +683,6 @@ def _get_default_ph3py_yaml(ph3py_yaml: Optional[Phono3pyYaml]): return _ph3py_yaml -def _create_phono3py_phonon_fc2( - phono3py: Phono3py, - ph3py_yaml: Optional[Phono3pyYaml], - symmetrize_fc2: bool, - is_compact_fc: bool, - fc_calculator: Optional[str], - fc_calculator_options: Optional[str], - log_level: int, -): - _ph3py_yaml = _get_default_ph3py_yaml(ph3py_yaml) - - try: - dataset = parse_forces( - phono3py, - ph3py_yaml=_ph3py_yaml, - force_filename="FORCES_FC2", - fc_type="phonon_fc2", - log_level=log_level, - ) - except RuntimeError as e: - if log_level: - print(str(e)) - print_error() - sys.exit(1) - except FileNotFoundError as e: - file_exists(e.filename, log_level) - - phono3py.phonon_dataset = dataset - phono3py.produce_fc2( - symmetrize_fc2=symmetrize_fc2, - is_compact_fc=is_compact_fc, - fc_calculator=extract_fc2_fc3_calculators(fc_calculator, 2), - fc_calculator_options=extract_fc2_fc3_calculators(fc_calculator_options, 2), - ) - - def _convert_unit_in_dataset( dataset: dict, distance_to_A: Optional[float] = None, diff --git a/phono3py/cui/load.py b/phono3py/cui/load.py index 238cfa31..14730c85 100644 --- a/phono3py/cui/load.py +++ b/phono3py/cui/load.py @@ -450,47 +450,33 @@ def compute_force_constants_from_datasets( fc2 : bool """ - if not read_fc["fc3"]: - if ph3py.dataset or ph3py.mlp_dataset: - if ph3py.mlp_dataset and use_pypolymlp: - run_pypolymlp_to_compute_forces( - ph3py, - mlp_params=mlp_params, - displacement_distance=displacement_distance, - number_of_snapshots=number_of_snapshots, - random_seed=random_seed, - log_level=log_level, - ) - - ph3py.produce_fc3( - symmetrize_fc3r=symmetrize_fc, - is_compact_fc=is_compact_fc, - fc_calculator=extract_fc2_fc3_calculators(fc_calculator, 3), - fc_calculator_options=extract_fc2_fc3_calculators( - fc_calculator_options, 3 - ), + fc3_calculator = extract_fc2_fc3_calculators(fc_calculator, 3) + fc2_calculator = extract_fc2_fc3_calculators(fc_calculator, 2) + if not read_fc["fc3"] and (ph3py.dataset or ph3py.mlp_dataset): + if ph3py.mlp_dataset and use_pypolymlp: + run_pypolymlp_to_compute_forces( + ph3py, + mlp_params=mlp_params, + displacement_distance=displacement_distance, + number_of_snapshots=number_of_snapshots, + random_seed=random_seed, + log_level=log_level, ) + ph3py.produce_fc3( + symmetrize_fc3r=symmetrize_fc, + is_compact_fc=is_compact_fc, + fc_calculator=fc3_calculator, + fc_calculator_options=extract_fc2_fc3_calculators(fc_calculator_options, 3), + ) - if log_level and symmetrize_fc and fc_calculator is None: - print("fc3 was symmetrized.") + if log_level and symmetrize_fc and fc_calculator is None: + print("fc3 was symmetrized.") if not read_fc["fc2"] and (ph3py.dataset or ph3py.phonon_dataset): - # if ( - # ph3py.phonon_supercell_matrix is None - # and fc_calculator in ("alm", "symfc") - # and ph3py.fc2 is not None - # ): - # if log_level: - # if fc_calculator == "alm": - # print("fc2 that was fit simultaneously with fc3 by ALM is used.") - # elif fc_calculator == "symfc": - # print( - # "fc2 that was fit simultaneously with fc3 by symfc is used.") - # else: ph3py.produce_fc2( symmetrize_fc2=symmetrize_fc, is_compact_fc=is_compact_fc, - fc_calculator=extract_fc2_fc3_calculators(fc_calculator, 2), + fc_calculator=fc2_calculator, fc_calculator_options=extract_fc2_fc3_calculators(fc_calculator_options, 2), ) if log_level and symmetrize_fc and fc_calculator is None: diff --git a/phono3py/cui/phono3py_script.py b/phono3py/cui/phono3py_script.py index 2fda8b0c..a1fc34be 100644 --- a/phono3py/cui/phono3py_script.py +++ b/phono3py/cui/phono3py_script.py @@ -160,11 +160,18 @@ def finalize_phono3py( yaml_filename = filename _physical_units = get_default_physical_units(phono3py.calculator) + + write_force_sets = phono3py.mlp is not None + _write_displacements = write_displacements or phono3py.mlp is not None + ph3py_yaml = Phono3pyYaml( configuration=confs_dict, calculator=phono3py.calculator, physical_units=_physical_units, - settings={"force_sets": False, "displacements": write_displacements}, + settings={ + "force_sets": write_force_sets, + "displacements": _write_displacements, + }, ) ph3py_yaml.set_phonon_info(phono3py) with open(yaml_filename, "w") as w: diff --git a/phono3py/interface/phono3py_yaml.py b/phono3py/interface/phono3py_yaml.py index 2decda01..caf9c493 100644 --- a/phono3py/interface/phono3py_yaml.py +++ b/phono3py/interface/phono3py_yaml.py @@ -326,7 +326,7 @@ def _displacements_yaml_lines(self, with_forces: bool = False) -> list: """ lines = [] - if self._data.phonon_supercell_matrix is not None: + if self._data.phonon_dataset is not None: lines += self._displacements_yaml_lines_2types( self._data.phonon_dataset, with_forces=with_forces, From 4f95be3e830c29b904c649b2fd3746d0614b7c03 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Sat, 29 Jun 2024 16:47:11 +0900 Subject: [PATCH 2/2] Update phonopy version dependency --- requirements.txt | 2 +- setup.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/requirements.txt b/requirements.txt index b7fdc84d..146789f8 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.24,<2.25 +phonopy >=2.25,<2.26 diff --git a/setup.py b/setup.py index d08cc8b5..1784f761 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.24,<2.25", + "phonopy>=2.25,<2.26", ], provides=["phono3py"], scripts=scripts_phono3py,