Implement thermal expansion for quasi-harmonic and energy-volume curve

atomistics/shared/thermo/thermalexpansion.py

@@ -0,0 +1,17 @@
+import numpy as np
+
+from atomistics.shared.thermo.debye import get_debye_model
+from atomistics.shared.thermo.thermo import get_thermo_bulk_model
+
+
+def get_thermal_expansion_with_evcurve(
+    fit_dict, masses, t_min=1, t_max=1500, t_step=50, temperatures=None
+):
+    if temperatures is None:
+        temperatures = np.arange(t_min, t_max + t_step, t_step)
+    debye_model = get_debye_model(fit_dict=fit_dict, masses=masses, num_steps=50)
+    pes = get_thermo_bulk_model(
+        temperatures=temperatures,
+        debye_model=debye_model,
+    )
+    return temperatures, pes.get_minimum_energy_path()

atomistics/workflows/evcurve/workflow.py

@@ -4,6 +4,7 @@
 
 from atomistics.workflows.evcurve.fit import EnergyVolumeFit
 from atomistics.workflows.shared.workflow import Workflow
+from atomistics.shared.thermo.thermalexpansion import get_thermal_expansion_with_evcurve
 
 
 def _strain_axes(
@@ -150,3 +151,16 @@ def analyse_structures(self, output_dict):
 
     def get_volume_lst(self):
         return get_volume_lst(structure_dict=self._structure_dict)
+
+    def get_thermal_expansion(
+        self, output_dict, t_min=1, t_max=1500, t_step=50, temperatures=None
+    ):
+        fit_dict = self.analyse_structures(output_dict=output_dict)
+        return get_thermal_expansion_with_evcurve(
+            fit_dict=fit_dict,
+            masses=self.structure.get_masses(),
+            t_min=t_min,
+            t_max=t_max,
+            t_step=t_step,
+            temperatures=temperatures,
+        )

atomistics/workflows/phonons/workflow.py

@@ -218,3 +218,8 @@ def plot_dos(self, *args, axis=None, **kwargs):
             axis=axis,
             **kwargs,
         )
+
+    def get_thermal_expansion(
+        self, output_dict, t_min=1, t_max=1500, t_step=50, temperatures=None
+    ):
+        raise NotImplementedError()

atomistics/workflows/quasiharmonic/workflow.py

@@ -1,3 +1,5 @@
+import numpy as np
+
 from atomistics.workflows.evcurve.workflow import EnergyVolumeCurveWorkflow
 from atomistics.workflows.phonons.workflow import PhonopyWorkflow
 from atomistics.workflows.phonons.units import VaspToTHz
@@ -91,3 +93,35 @@ def get_thermal_properties(self, t_min=1, t_max=1500, t_step=50, temperatures=No
                 t_step=t_step, t_max=t_max, t_min=t_min, temperatures=temperatures
             )
         return tp_collect_dict
+
+    def get_thermal_expansion(
+        self, output_dict, t_min=1, t_max=1500, t_step=50, temperatures=None
+    ):
+        (
+            eng_internal_dict,
+            mesh_collect_dict,
+            dos_collect_dict,
+        ) = 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
+        )
+
+        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()))
+
+        fit_deg = 4
+        vol_best = volume_lst[int(len(volume_lst) / 2)]
+        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
+            )
+            p = np.polyfit(volume_lst, free_eng_lst, deg=fit_deg)
+            extrema = np.roots(np.polyder(p, m=1)).real
+            vol_select = extrema[np.argmin(np.abs(extrema - vol_best))]
+            eng_lst.append(np.poly1d(p)(vol_select))
+            vol_lst.append(vol_select)
+        return temperatures, vol_lst