From adc76bf736448e198487237447ce6db6ee804f2a Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Wed, 13 Dec 2023 16:49:56 +0100 Subject: [PATCH 1/3] Implement get_thermal_properties() for quasi-harmonic analog to the harmonic implementation --- .../workflows/quasiharmonic/workflow.py | 205 ++++++++++++------ 1 file changed, 137 insertions(+), 68 deletions(-) diff --git a/atomistics/workflows/quasiharmonic/workflow.py b/atomistics/workflows/quasiharmonic/workflow.py index 3fbf41f3..3af8d0bb 100644 --- a/atomistics/workflows/quasiharmonic/workflow.py +++ b/atomistics/workflows/quasiharmonic/workflow.py @@ -48,8 +48,9 @@ def __init__( # This provides some computational efficiency for classical calculations # And for quantum calculations _ensures_ that force matrices and energy/atom # get treated with the same kmesh + structure_repeated = structure.repeat(repeat_vector) super().__init__( - structure=structure.repeat(repeat_vector), + structure=structure_repeated, num_points=num_points, fit_type=fit_type, fit_order=fit_order, @@ -64,6 +65,8 @@ def __init__( self._factor = factor self._primitive_matrix = primitive_matrix self._phonopy_dict = {} + self._volume_rescale_factor = len(structure_repeated) / len(structure) + self._eng_internal_dict = None def generate_structures(self): task_dict = super().generate_structures() @@ -90,7 +93,7 @@ def generate_structures(self): return task_dict def analyse_structures(self, output_dict): - eng_internal_dict = output_dict["energy"] + self._eng_internal_dict = output_dict["energy"] mesh_collect_dict, dos_collect_dict = {}, {} for strain, phono in self._phonopy_dict.items(): mesh_dict, dos_dict = phono.analyse_structures( @@ -100,7 +103,7 @@ def analyse_structures(self, output_dict): ) mesh_collect_dict[strain] = mesh_dict dos_collect_dict[strain] = dos_dict - return eng_internal_dict, mesh_collect_dict, dos_collect_dict + return self._eng_internal_dict, mesh_collect_dict, dos_collect_dict def get_thermal_properties( self, @@ -112,6 +115,7 @@ def get_thermal_properties( pretend_real=False, band_indices=None, is_projection=False, + quantum_mechanical=True, ): """ Returns thermal properties at constant volume in the given temperature range. Can only be called after job @@ -127,21 +131,117 @@ def get_thermal_properties( Returns: :class:`Thermal`: thermal properties as returned by Phonopy """ - tp_collect_dict = {} - for strain, phono in self._phonopy_dict.items(): - tp_collect_dict[strain] = phono.get_thermal_properties( - t_step=t_step, - t_max=t_max, + if self._eng_internal_dict is None: + raise ValueError("Please first execute analyse_output() before calling get_thermal_properties().") + if quantum_mechanical: + tp_collect_dict = self._get_thermal_properties_quantum_mechanical( t_min=t_min, + t_max=t_max, + t_step=t_step, temperatures=temperatures, cutoff_frequency=cutoff_frequency, pretend_real=pretend_real, band_indices=band_indices, is_projection=is_projection, ) + else: + if is_projection: + raise ValueError( + "is_projection!=False is incompatible to quantum_mechanical=False." + ) + if pretend_real: + raise ValueError( + "pretend_real!=False is incompatible to quantum_mechanical=False." + ) + if band_indices is not None: + raise ValueError( + "band_indices!=None is incompatible to quantum_mechanical=False." + ) + tp_collect_dict = self._get_thermal_properties_classical( + t_min=t_min, + t_max=t_max, + t_step=t_step, + temperatures=temperatures, + cutoff_frequency=cutoff_frequency, + ) + + temperatures = tp_collect_dict[1.0]["temperatures"] + quantities = tp_collect_dict[1.0].keys() + strain_lst = self._eng_internal_dict.keys() + volume_lst = np.array(self.get_volume_lst()) / self._volume_rescale_factor + eng_int_lst = np.array(list(self._eng_internal_dict.values())) / self._volume_rescale_factor + + vol_lst, eng_lst = [], [] + for i, temp in enumerate(temperatures): + free_eng_lst = ( + np.array([tp_collect_dict[s]["free_energy"][i] for s in strain_lst]) + + eng_int_lst + ) + fit_dict = fit_ev_curve( + volume_lst=volume_lst, + energy_lst=free_eng_lst, + fit_type=self.fit_type, + fit_order=self.fit_order, + ) + eng_lst.append(fit_dict["energy_eq"]) + vol_lst.append(fit_dict["volume_eq"]) + + mat_dict = {} + for k in quantities: + if k != "temperatures": + mat_dict[k] = np.array([ + np.poly1d(np.polyfit(volume_lst, q_over_v, 1))(vol_opt) + for q_over_v, vol_opt in zip( + np.array([tp_collect_dict[s][k] for s in strain_lst]).T, + vol_lst + ) + ]) + mat_dict["volumes"] = vol_lst + mat_dict["temperatures"] = temperatures + return mat_dict + + def _get_thermal_properties_quantum_mechanical( + self, + t_min=1, + t_max=1500, + t_step=50, + temperatures=None, + cutoff_frequency=None, + pretend_real=False, + band_indices=None, + is_projection=False, + ): + """ + Returns thermal properties at constant volume in the given temperature range. Can only be called after job + successfully ran. + + Args: + t_min (float): minimum sample temperature + t_max (float): maximum sample temperature + t_step (int): temperature sample interval + temperatures (array_like, float): custom array of temperature samples, if given t_min, t_max, t_step are + ignored. + + Returns: + :class:`Thermal`: thermal properties as returned by Phonopy + """ + tp_collect_dict = {} + for strain, phono in self._phonopy_dict.items(): + tp_collect_dict[strain] = { + k: v if k == "temperatures" else v / self._volume_rescale_factor + for k, v in phono.get_thermal_properties( + t_step=t_step, + t_max=t_max, + t_min=t_min, + temperatures=temperatures, + cutoff_frequency=cutoff_frequency, + pretend_real=pretend_real, + band_indices=band_indices, + is_projection=is_projection, + ).items()} return tp_collect_dict - def get_thermal_properties_classical( + def _get_thermal_properties_classical( self, t_min=1, t_max=1500, @@ -149,6 +249,20 @@ def get_thermal_properties_classical( temperatures=None, cutoff_frequency=None, ): + """ + Returns thermal properties at constant volume in the given temperature range. Can only be called after job + successfully ran. + + Args: + t_min (float): minimum sample temperature + t_max (float): maximum sample temperature + t_step (int): temperature sample interval + temperatures (array_like, float): custom array of temperature samples, if given t_min, t_max, t_step are + ignored. + + Returns: + :class:`Thermal`: thermal properties as returned by Phonopy + """ if temperatures is None: temperatures = np.arange(t_min, t_max, t_step) if cutoff_frequency is None or cutoff_frequency < 0: @@ -174,7 +288,7 @@ def get_thermal_properties_classical( ) tp_collect_dict[strain] = { "temperatures": temperatures, - "free_energy": np.array(t_property_lst) * kJ_mol_to_eV, + "free_energy": np.array(t_property_lst) * kJ_mol_to_eV / self._volume_rescale_factor, } return tp_collect_dict @@ -191,61 +305,16 @@ def get_thermal_expansion( is_projection=False, quantum_mechanical=True, ): - ( - eng_internal_dict, - mesh_collect_dict, - dos_collect_dict, - ) = self.analyse_structures(output_dict=output_dict) - if quantum_mechanical: - tp_collect_dict = self.get_thermal_properties( - t_min=t_min, - t_max=t_max, - t_step=t_step, - temperatures=temperatures, - cutoff_frequency=cutoff_frequency, - pretend_real=pretend_real, - band_indices=band_indices, - is_projection=is_projection, - ) - else: - if is_projection: - raise ValueError( - "is_projection!=False is incompatible to quantum_mechanical=False." - ) - if pretend_real: - raise ValueError( - "pretend_real!=False is incompatible to quantum_mechanical=False." - ) - if band_indices is not None: - raise ValueError( - "band_indices!=None is incompatible to quantum_mechanical=False." - ) - - tp_collect_dict = self.get_thermal_properties_classical( - t_min=t_min, - t_max=t_max, - t_step=t_step, - temperatures=temperatures, - cutoff_frequency=cutoff_frequency, - ) - - temperatures = tp_collect_dict[1.0]["temperatures"] - strain_lst = eng_internal_dict.keys() - volume_lst = self.get_volume_lst() - eng_int_lst = np.array(list(eng_internal_dict.values())) - - vol_lst, eng_lst = [], [] - for i, temp in enumerate(temperatures): - free_eng_lst = ( - np.array([tp_collect_dict[s]["free_energy"][i] for s in strain_lst]) - + eng_int_lst - ) - fit_dict = fit_ev_curve( - volume_lst=volume_lst, - energy_lst=free_eng_lst, - fit_type=self.fit_type, - fit_order=self.fit_order, - ) - eng_lst.append(fit_dict["energy_eq"]) - vol_lst.append(fit_dict["volume_eq"]) - return temperatures, vol_lst + self.analyse_structures(output_dict=output_dict) + tp_collect_dict = self.get_thermal_properties( + t_min=t_min, + t_max=t_max, + t_step=t_step, + temperatures=temperatures, + cutoff_frequency=cutoff_frequency, + pretend_real=pretend_real, + band_indices=band_indices, + is_projection=is_projection, + quantum_mechanical=quantum_mechanical, + ) + return tp_collect_dict["temperatures"], tp_collect_dict["volumes"] From f110b6c477caf9b74bba7b659f82c0b95615a5ac Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Wed, 13 Dec 2023 16:50:36 +0100 Subject: [PATCH 2/3] fix black formatting --- atomistics/calculators/ase.py | 14 ++++---- atomistics/calculators/lammps/calculator.py | 14 ++++---- atomistics/calculators/qe.py | 22 ++++++------- .../workflows/quasiharmonic/workflow.py | 32 ++++++++++++------- 4 files changed, 46 insertions(+), 36 deletions(-) diff --git a/atomistics/calculators/ase.py b/atomistics/calculators/ase.py index 5f3658f0..f22a3178 100644 --- a/atomistics/calculators/ase.py +++ b/atomistics/calculators/ase.py @@ -29,13 +29,13 @@ def evaluate_with_ase( ase_optimizer_kwargs=ase_optimizer_kwargs, ) elif "optimize_positions_and_volume" in tasks: - results["structure_with_optimized_positions_and_volume"] = ( - optimize_positions_and_volume_with_ase( - structure=structure, - ase_calculator=ase_calculator, - ase_optimizer=ase_optimizer, - ase_optimizer_kwargs=ase_optimizer_kwargs, - ) + results[ + "structure_with_optimized_positions_and_volume" + ] = optimize_positions_and_volume_with_ase( + structure=structure, + ase_calculator=ase_calculator, + ase_optimizer=ase_optimizer, + ase_optimizer_kwargs=ase_optimizer_kwargs, ) elif "calc_energy" in tasks or "calc_forces" in tasks or "calc_stress" in tasks: if "calc_energy" in tasks and "calc_forces" in tasks and "calc_stress" in tasks: diff --git a/atomistics/calculators/lammps/calculator.py b/atomistics/calculators/lammps/calculator.py index 6c9d5278..6fee3952 100644 --- a/atomistics/calculators/lammps/calculator.py +++ b/atomistics/calculators/lammps/calculator.py @@ -362,13 +362,13 @@ def evaluate_with_lammps_library( ): results = {} if "optimize_positions_and_volume" in tasks: - results["structure_with_optimized_positions_and_volume"] = ( - optimize_positions_and_volume_with_lammps( - structure=structure, - potential_dataframe=potential_dataframe, - lmp=lmp, - **lmp_optimizer_kwargs, - ) + results[ + "structure_with_optimized_positions_and_volume" + ] = optimize_positions_and_volume_with_lammps( + structure=structure, + potential_dataframe=potential_dataframe, + lmp=lmp, + **lmp_optimizer_kwargs, ) elif "optimize_positions" in tasks: results["structure_with_optimized_positions"] = optimize_positions_with_lammps( diff --git a/atomistics/calculators/qe.py b/atomistics/calculators/qe.py index 47b2fd62..e4c14c22 100644 --- a/atomistics/calculators/qe.py +++ b/atomistics/calculators/qe.py @@ -275,17 +275,17 @@ def evaluate_with_qe( ): results = {} if "optimize_positions_and_volume" in tasks: - results["structure_with_optimized_positions_and_volume"] = ( - optimize_positions_and_volume_with_qe( - structure=structure, - calculation_name=calculation_name, - working_directory=working_directory, - kpts=kpts, - pseudopotentials=pseudopotentials, - tstress=tstress, - tprnfor=tprnfor, - **kwargs, - ) + results[ + "structure_with_optimized_positions_and_volume" + ] = optimize_positions_and_volume_with_qe( + structure=structure, + calculation_name=calculation_name, + working_directory=working_directory, + kpts=kpts, + pseudopotentials=pseudopotentials, + tstress=tstress, + tprnfor=tprnfor, + **kwargs, ) elif "calc_energy" in tasks or "calc_forces" in tasks: if "calc_energy" in tasks and "calc_forces" in tasks: diff --git a/atomistics/workflows/quasiharmonic/workflow.py b/atomistics/workflows/quasiharmonic/workflow.py index 3af8d0bb..3322b6bf 100644 --- a/atomistics/workflows/quasiharmonic/workflow.py +++ b/atomistics/workflows/quasiharmonic/workflow.py @@ -132,7 +132,9 @@ def get_thermal_properties( :class:`Thermal`: thermal properties as returned by Phonopy """ if self._eng_internal_dict is None: - raise ValueError("Please first execute analyse_output() before calling get_thermal_properties().") + raise ValueError( + "Please first execute analyse_output() before calling get_thermal_properties()." + ) if quantum_mechanical: tp_collect_dict = self._get_thermal_properties_quantum_mechanical( t_min=t_min, @@ -169,7 +171,10 @@ def get_thermal_properties( quantities = tp_collect_dict[1.0].keys() strain_lst = self._eng_internal_dict.keys() volume_lst = np.array(self.get_volume_lst()) / self._volume_rescale_factor - eng_int_lst = np.array(list(self._eng_internal_dict.values())) / self._volume_rescale_factor + eng_int_lst = ( + np.array(list(self._eng_internal_dict.values())) + / self._volume_rescale_factor + ) vol_lst, eng_lst = [], [] for i, temp in enumerate(temperatures): @@ -189,13 +194,15 @@ def get_thermal_properties( mat_dict = {} for k in quantities: if k != "temperatures": - mat_dict[k] = np.array([ - np.poly1d(np.polyfit(volume_lst, q_over_v, 1))(vol_opt) - for q_over_v, vol_opt in zip( - np.array([tp_collect_dict[s][k] for s in strain_lst]).T, - vol_lst - ) - ]) + mat_dict[k] = np.array( + [ + np.poly1d(np.polyfit(volume_lst, q_over_v, 1))(vol_opt) + for q_over_v, vol_opt in zip( + np.array([tp_collect_dict[s][k] for s in strain_lst]).T, + vol_lst, + ) + ] + ) mat_dict["volumes"] = vol_lst mat_dict["temperatures"] = temperatures return mat_dict @@ -238,7 +245,8 @@ def _get_thermal_properties_quantum_mechanical( pretend_real=pretend_real, band_indices=band_indices, is_projection=is_projection, - ).items()} + ).items() + } return tp_collect_dict def _get_thermal_properties_classical( @@ -288,7 +296,9 @@ def _get_thermal_properties_classical( ) tp_collect_dict[strain] = { "temperatures": temperatures, - "free_energy": np.array(t_property_lst) * kJ_mol_to_eV / self._volume_rescale_factor, + "free_energy": np.array(t_property_lst) + * kJ_mol_to_eV + / self._volume_rescale_factor, } return tp_collect_dict From cca83c7923936a604a3a6b8159767289a20768a3 Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Wed, 13 Dec 2023 17:19:03 +0100 Subject: [PATCH 3/3] fix tests --- tests/test_quasiharmonic_ase_emt.py | 2 +- tests/test_quasiharmonic_lammps.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_quasiharmonic_ase_emt.py b/tests/test_quasiharmonic_ase_emt.py index de2b4f42..30c304e0 100644 --- a/tests/test_quasiharmonic_ase_emt.py +++ b/tests/test_quasiharmonic_ase_emt.py @@ -47,7 +47,7 @@ def test_calc_phonons(self): quantum_mechanical=False ) self.assertEqual(len(eng_internal_dict.keys()), 11) - self.assertEqual(len(tp_collect_dict.keys()), 11) + self.assertEqual(len(tp_collect_dict.keys()), 5) self.assertEqual(len(temperatures_qh_qm), 2) self.assertEqual(len(volumes_qh_qm), 2) self.assertTrue(volumes_qh_qm[0] < volumes_qh_qm[-1]) diff --git a/tests/test_quasiharmonic_lammps.py b/tests/test_quasiharmonic_lammps.py index f8fbbdb7..a3c53c28 100644 --- a/tests/test_quasiharmonic_lammps.py +++ b/tests/test_quasiharmonic_lammps.py @@ -60,7 +60,7 @@ def test_calc_phonons(self): quantum_mechanical=False ) self.assertEqual(len(eng_internal_dict.keys()), 11) - self.assertEqual(len(tp_collect_dict.keys()), 11) + self.assertEqual(len(tp_collect_dict.keys()), 5) self.assertEqual(len(temperatures_qh_qm), 2) self.assertEqual(len(volumes_qh_qm), 2) self.assertTrue(volumes_qh_qm[0] < volumes_qh_qm[-1])