Skip to content

Commit

Permalink
Merge pull request #224 from pyiron/qh_rescale_fix
Browse files Browse the repository at this point in the history 
QuasiHarmonic: Fix strain
  • Loading branch information
@jan-janssen
jan-janssen authored Mar 18, 2024
2 parents 362c429 + 096caa0 commit aa3db9a
Show file tree
Hide file tree
Showing 3 changed files with 63 additions and 66 deletions.
19 changes: 11 additions & 8 deletions atomistics/workflows/evcurve/workflow.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -164,14 +164,7 @@ def generate_structures(self) -> dict:
Returns:
(dict)
"""
strains = self.strains
if strains is None:
strains = np.linspace(
-self.vol_range,
self.vol_range,
int(self.num_points),
)
for strain in strains:
for strain in self._get_strains():
basis = _strain_axes(
structure=self.structure, axes=self.axes, volume_strain=strain
)
Expand Down Expand Up @@ -218,3 +211,13 @@ def get_thermal_properties(
constant_volume=constant_volume,
output_keys=output_keys,
)

def _get_strains(self):
strains = self.strains
if strains is None:
strains = np.linspace(
-self.vol_range,
self.vol_range,
int(self.num_points),
)
return strains
98 changes: 46 additions & 52 deletions atomistics/workflows/quasiharmonic.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
from atomistics.workflows.evcurve.workflow import (
EnergyVolumeCurveWorkflow,
fit_ev_curve,
_strain_axes,
)
from atomistics.workflows.phonons.workflow import PhonopyWorkflow
from atomistics.workflows.phonons.helper import get_supercell_matrix
Expand All @@ -27,7 +28,6 @@ def get_thermal_properties(
eng_internal_dict: dict,
phonopy_dict: dict,
volume_lst: np.ndarray,
volume_rescale_factor: float,
fit_type: str,
fit_order: int,
t_min: float = 1.0,
Expand Down Expand Up @@ -58,7 +58,6 @@ def get_thermal_properties(
if quantum_mechanical:
tp_collect_dict = _get_thermal_properties_quantum_mechanical(
phonopy_dict=phonopy_dict,
volume_rescale_factor=volume_rescale_factor,
t_min=t_min,
t_max=t_max,
t_step=t_step,
Expand All @@ -84,7 +83,6 @@ def get_thermal_properties(
)
tp_collect_dict = _get_thermal_properties_classical(
phonopy_dict=phonopy_dict,
volume_rescale_factor=volume_rescale_factor,
t_min=t_min,
t_max=t_max,
t_step=t_step,
Expand All @@ -94,7 +92,7 @@ def get_thermal_properties(

temperatures = tp_collect_dict[1.0]["temperatures"]
strain_lst = eng_internal_dict.keys()
eng_int_lst = np.array(list(eng_internal_dict.values())) / volume_rescale_factor
eng_int_lst = np.array(list(eng_internal_dict.values()))

vol_lst, eng_lst = [], []
for i, temp in enumerate(temperatures):
Expand Down Expand Up @@ -129,7 +127,6 @@ def get_thermal_properties(

def _get_thermal_properties_quantum_mechanical(
phonopy_dict: dict,
volume_rescale_factor: float,
t_min: float = 1.0,
t_max: float = 1500.0,
t_step: float = 50.0,
Expand All @@ -154,28 +151,24 @@ def _get_thermal_properties_quantum_mechanical(
Returns:
:class:`Thermal`: thermal properties as returned by Phonopy
"""
tp_collect_dict = {}
for strain, phono in phonopy_dict.items():
tp_collect_dict[strain] = {
k: v if k == "temperatures" else v / 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,
output_keys=output_keys,
).items()
}
return tp_collect_dict
return {
strain: 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,
output_keys=output_keys,
)
for strain, phono in phonopy_dict.items()
}


def _get_thermal_properties_classical(
phonopy_dict: dict,
volume_rescale_factor: float,
t_min: float = 1.0,
t_max: float = 1500.0,
t_step: float = 50.0,
Expand Down Expand Up @@ -223,9 +216,7 @@ def _get_thermal_properties_classical(
)
tp_collect_dict[strain] = {
"temperatures": temperatures,
"free_energy": np.array(t_property_lst)
* kJ_mol_to_eV
/ volume_rescale_factor,
"free_energy": np.array(t_property_lst) * kJ_mol_to_eV,
}
return tp_collect_dict

Expand Down Expand Up @@ -292,23 +283,8 @@ def __init__(
primitive_matrix: np.ndarray = None,
number_of_snapshots: int = None,
):
repeat_vector = np.array(
np.diag(
get_supercell_matrix(
interaction_range=interaction_range,
cell=structure.cell.array,
)
),
dtype=int,
)
# Phonopy internally repeats structures that are "too small"
# Here we manually guarantee that all structures passed are big enough
# 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_repeated,
structure=structure,
num_points=num_points,
fit_type=fit_type,
fit_order=fit_order,
Expand All @@ -323,31 +299,50 @@ 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
# Phonopy internally repeats structures that are "too small"
# Here we manually guarantee that all structures passed are big enough
# 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
self._repeat_vector = np.array(
np.diag(
get_supercell_matrix(
interaction_range=interaction_range,
cell=structure.cell.array,
)
),
dtype=int,
)

def generate_structures(self) -> dict:
task_dict = super().generate_structures()
task_dict["calc_forces"] = {}
for strain, structure in task_dict["calc_energy"].items():
self._phonopy_dict[strain] = PhonopyWorkflow(
structure=structure,
task_dict = {"calc_forces": {}}
for strain in self._get_strains():
strain_ind = 1 + np.round(strain, 7)
basis = _strain_axes(
structure=self.structure, axes=self.axes, volume_strain=strain
)
structure_ev = basis.repeat(self._repeat_vector)
self._structure_dict[strain_ind] = structure_ev
self._phonopy_dict[strain_ind] = PhonopyWorkflow(
structure=basis,
interaction_range=self._interaction_range,
factor=self._factor,
displacement=self._displacement,
dos_mesh=self._dos_mesh,
primitive_matrix=self._primitive_matrix,
number_of_snapshots=self._number_of_snapshots,
)
structure_task_dict = self._phonopy_dict[strain].generate_structures()
structure_task_dict = self._phonopy_dict[strain_ind].generate_structures()
task_dict["calc_forces"].update(
{
(strain, key): structure_phono
(strain_ind, key): structure_phono
for key, structure_phono in structure_task_dict[
"calc_forces"
].items()
}
)
task_dict["calc_energy"] = self._structure_dict
return task_dict

def analyse_structures(
Expand Down Expand Up @@ -401,8 +396,7 @@ def get_thermal_properties(
return get_thermal_properties(
eng_internal_dict=self._eng_internal_dict,
phonopy_dict=self._phonopy_dict,
volume_lst=np.array(self.get_volume_lst()) / self._volume_rescale_factor,
volume_rescale_factor=self._volume_rescale_factor,
volume_lst=np.array(self.get_volume_lst()) / np.prod(self._repeat_vector),
fit_type=self.fit_type,
fit_order=self.fit_order,
t_min=t_min,
Expand Down
12 changes: 6 additions & 6 deletions tests/test_quasiharmonic_lammps.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,16 +59,16 @@ def test_calc_phonons(self):
self.assertTrue(tp_collect_dict["free_energy"][-1] > -2.7)
self.assertTrue(tp_collect_dict["entropy"][0] < 0.1)
self.assertTrue(tp_collect_dict["entropy"][0] > 0.0)
self.assertTrue(tp_collect_dict["entropy"][-1] < 273)
self.assertTrue(tp_collect_dict["entropy"][-1] > 272)
self.assertTrue(tp_collect_dict["entropy"][-1] < 270)
self.assertTrue(tp_collect_dict["entropy"][-1] > 269)
self.assertTrue(tp_collect_dict["heat_capacity"][0] < 0.1)
self.assertTrue(tp_collect_dict["heat_capacity"][0] > 0.0)
self.assertTrue(tp_collect_dict["heat_capacity"][-1] < 100)
self.assertTrue(tp_collect_dict["heat_capacity"][-1] > 99)
self.assertTrue(tp_collect_dict["volumes"][-1] < 68.6)
self.assertTrue(tp_collect_dict["volumes"][-1] > 68.5)
self.assertTrue(tp_collect_dict["volumes"][0] < 66.8)
self.assertTrue(tp_collect_dict["volumes"][0] > 66.7)
self.assertTrue(tp_collect_dict["volumes"][-1] < 66.6)
self.assertTrue(tp_collect_dict["volumes"][-1] > 66.5)
self.assertTrue(tp_collect_dict["volumes"][0] < 66.5)
self.assertTrue(tp_collect_dict["volumes"][0] > 66.4)
thermal_properties_dict = workflow.get_thermal_properties(
temperatures=[100, 1000],
output_keys=["temperatures", "volumes"],
Expand Down

0 comments on commit aa3db9a

Please sign in to comment.