Add plot_qha_v-ZSISA.py and preliminary tests in test_qha_approximati…

…on.py
abinit · Dec 12, 2024 · a7dc13b · a7dc13b
1 parent 2e41eac
commit a7dc13b
Show file tree

Hide file tree

Showing 3 changed files with 191 additions and 12 deletions.
diff --git a/abipy/dfpt/qha_aproximation.py b/abipy/dfpt/qha_aproximation.py
@@ -630,7 +630,8 @@ def plot_thermal_expansion_coeff(self, tstart=0, tstop=1000, num=101, tref=None,
         ax.get_yaxis().get_major_formatter().set_powerlimits((0, 0))
 
         return fig
-    def get_abc(self, tstart=0, tstop=1000, num=101,volumes="volumes"):
+
+    def get_abc(self, tstart=0, tstop=1000, num=101, volumes="volumes"):
         """
         Plots the thermal expansion coefficient as a function of the temperature.
 
@@ -655,9 +656,9 @@ def get_abc(self, tstart=0, tstop=1000, num=101,volumes="volumes"):
         pc = np.poly1d(param)
         cc_qha=pc(volumes)
 
-        return aa_qha,bb_qha,cc_qha
+        return aa_qha, bb_qha, cc_qha
 
-    def get_angles(self, tstart=0, tstop=1000, num=101,volumes="volumes"):
+    def get_angles(self, tstart=0, tstop=1000, num=101, volumes="volumes"):
         """
         Plots the thermal expansion coefficient as a function of the temperature.
 
@@ -682,7 +683,7 @@ def get_angles(self, tstart=0, tstop=1000, num=101,volumes="volumes"):
         pc = np.poly1d(param)
         alpha=pc(volumes)
 
-        return alpha,beta,gamma
+        return alpha, beta, gamma
 
 ###################################################################################################
     @add_fig_kwargs
@@ -738,7 +739,7 @@ def plot_thermal_expansion_coeff_abc(self, tstart=0, tstop=1000, num=101, tref=N
         tmesh = np.linspace(tstart, tstop, num)
         dt= tmesh[1] - tmesh[0]
 
-        aa,bb,cc = self.get_abc(tstart, tstop, num,vol)
+        aa,bb,cc = self.get_abc(tstart, tstop, num, vol)
         if (tref!=None):
             aa_tref,bb_tref,cc_tref = self.get_abc(tref, tref, 1,vol_tref)
 
@@ -794,6 +795,7 @@ def plot_thermal_expansion_coeff_abc(self, tstart=0, tstop=1000, num=101, tref=N
         ax.get_yaxis().get_major_formatter().set_powerlimits((0, 0))
 
         return fig
+
     @add_fig_kwargs
     def plot_thermal_expansion_coeff_angles(self, tstart=0, tstop=1000, num=101, tref=None, ax=None, **kwargs):
         """
@@ -976,6 +978,7 @@ def plot_abc_vs_t(self, tstart=0, tstop=1000, num=101, lattice=None, tref=None,
         ax.get_yaxis().get_major_formatter().set_powerlimits((0, 0))
 
         return fig
+
     @add_fig_kwargs
     def plot_angles_vs_t(self, tstart=0, tstop=1000, num=101, angle=None, tref=None, ax=None, **kwargs):
         """
@@ -1049,8 +1052,9 @@ def plot_angles_vs_t(self, tstart=0, tstop=1000, num=101, angle=None, tref=None,
         ax.get_yaxis().get_major_formatter().set_powerlimits((0, 0))
 
         return fig
+
     #&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
-    def fit_forth(self, tstart=0, tstop=1000, num=1,energy="energy",volumes="volumes"):
+    def fit_forth(self, tstart=0, tstop=1000, num=1, energy="energy", volumes="volumes"):
         """
         Performs a fit of the energies as a function of the volume at different temperatures.
 
@@ -1145,6 +1149,7 @@ def vol_EinfVib1_forth(self, tstart=0, tstop=1000, num=101):
         vol=f.min_vol
 
         return vol
+
     def vol_Einf_Vib2_forth(self, tstart=0, tstop=1000, num=101):
         """
         Args:
@@ -1181,6 +1186,7 @@ def vol_Einf_Vib2_forth(self, tstart=0, tstop=1000, num=101):
         vol=f.min_vol
 
         return vol
+
     def vol_Einf_Vib4_forth(self, tstart=0, tstop=1000, num=101):
         """
 
@@ -1191,7 +1197,6 @@ def vol_Einf_Vib4_forth(self, tstart=0, tstop=1000, num=101):
 
         Returns: |Vol|
         """
-
         tmesh = np.linspace(tstart, tstop, num)
         ph_energies = self.get_vib_free_energies(tstart, tstop, num)
         energies = self.energies
@@ -1223,7 +1228,7 @@ def vol_Einf_Vib4_forth(self, tstart=0, tstop=1000, num=101):
         vol=f.min_vol
 
         return vol
-#$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+
     @add_fig_kwargs
     def plot_vol_vs_t_4th(self, tstart=0, tstop=1000, num=101, ax=None, **kwargs):
         """
@@ -1332,7 +1337,8 @@ def get_thermal_expansion_coeff_4th(self, tstart=0, tstop=1000, num=101 , tref=N
             #alpha= - 1/f0.min_vol * d2f_t_v[1:-1] / F2D[1:-1]
             alpha= - 1/f0.min_vol * d2f_t_v / F2D
 
-        return  alpha
+        return alpha
+
     @add_fig_kwargs
     def plot_thermal_expansion_coeff_4th(self, tstart=0, tstop=1000, num=101, tref=None, ax=None, **kwargs):
         """
@@ -1452,6 +1458,7 @@ def plot_thermal_expansion_coeff_4th(self, tstart=0, tstop=1000, num=101, tref=N
         ax.get_yaxis().get_major_formatter().set_powerlimits((0, 0))
 
         return fig
+
     @add_fig_kwargs
     def plot_abc_vs_t_4th(self, tstart=0, tstop=1000, num=101, lattice=None,tref=None, ax=None, **kwargs):
         """
@@ -1525,6 +1532,7 @@ def plot_abc_vs_t_4th(self, tstart=0, tstop=1000, num=101, lattice=None,tref=Non
         ax.get_yaxis().get_major_formatter().set_powerlimits((0, 0))
 
         return fig
+
     @add_fig_kwargs
     def plot_angles_vs_t_4th(self, tstart=0, tstop=1000, num=101,angle=None,  tref=None, ax=None, **kwargs):
         """
@@ -1598,7 +1606,7 @@ def plot_angles_vs_t_4th(self, tstart=0, tstop=1000, num=101,angle=None,  tref=N
         ax.get_yaxis().get_major_formatter().set_powerlimits((0, 0))
 
         return fig
-###################################################################################################
+
     @add_fig_kwargs
     def plot_thermal_expansion_coeff_abc_4th(self, tstart=0, tstop=1000, num=101, tref=None, ax=None, **kwargs):
         """
@@ -1708,6 +1716,7 @@ def plot_thermal_expansion_coeff_abc_4th(self, tstart=0, tstop=1000, num=101, tr
         ax.get_yaxis().get_major_formatter().set_powerlimits((0, 0))
 
         return fig
+
     @add_fig_kwargs
     def plot_thermal_expansion_coeff_angles_4th(self, tstart=0, tstop=1000, num=101, tref=None, ax=None, **kwargs):
         """
@@ -1816,7 +1825,7 @@ def plot_thermal_expansion_coeff_angles_4th(self, tstart=0, tstop=1000, num=101,
         ax.get_yaxis().get_major_formatter().set_powerlimits((0, 0))
 
         return fig
-#*********************************************************************************************
+
     @classmethod
     def from_files_app(cls, gsr_paths, phdos_paths):
         """
@@ -1916,7 +1925,7 @@ def get_thermodynamic_properties(self, tstart=0, tstop=1000, num=101):
             zpe[i] = d.zero_point_energy
 
         return dict2namedtuple(tmesh=tmesh, cv=cv, free_energy=free_energy, entropy=entropy, zpe=zpe)
-#=======================================================================================================
+
     @classmethod
     def from_files_app_ddb(cls, ddb_paths, phdos_paths):
         """

diff --git a/abipy/dfpt/tests/test_qha_approximation.py b/abipy/dfpt/tests/test_qha_approximation.py
@@ -0,0 +1,103 @@
+"""Tests for frozen_phonons"""
+import os
+import abipy.data as abidata
+
+from abipy.dfpt.qha_aproximation import QHA_App
+from abipy.core.testing import AbipyTest
+
+
+class QhaTest(AbipyTest):
+
+    def test_v_ZSIZA(self):
+
+        # Root points to the directory in the git submodule with the output results.
+        root = os.path.join(abidata.dirpath, "data_v-ZSISA-QHA.git", "Si_v_ZSISA_approximation")
+
+        strains = [96, 98, 100, 102, 104, 106]
+        strains2 = [98, 100, 102, 104, 106] # EinfVib4(D)
+        #strains2 = [96, 98, 100, 102, 104] # EinfVib4(S)
+        #strains2 = [100, 102, 104] # EinfVib2(D)
+
+        gsr_paths = [os.path.join(root, "scale_{:d}_GSR.nc".format(s)) for s in strains]
+        ddb_paths = [os.path.join(root, "scale_{:d}_GSR_DDB".format(s)) for s in strains]
+        dos_paths = [os.path.join(root, "scale_{:d}_PHDOS.nc".format(s)) for s in strains2]
+
+        qha = QHA_App.from_files_app_ddb(ddb_paths, dos_paths)
+        tstart, tstop = 0, 800
+
+        # Test basic properties and get methods of qha
+        assert qha.nvols == 5
+        f = qha.get_thermal_expansion_coeff(tstart=0, tstop=1000, num=7, tref=None)
+        self.assert_almost_equal(f.values, [-0.0000000e+00,  5.1944952e-07,  7.0431284e-06,  9.4401141e-06,
+                                             1.0630764e-05,  1.1414897e-05,  1.2035432e-05])
+
+        # FIXME: Calling these functions directly raises:
+        #numpy.core._exceptions._UFuncNoLoopError: ufunc 'multiply' did not contain a loop with signature matching types (dtype('<U7'), dtype('<U7')) -> None
+
+        #aas, bbs, ccs = qha.get_abc(tstart=0, tstop=1000, num=7, volumes="volumes")
+        #self.assert_almost_equal(aas, [0])
+        #self.assert_almost_equal(bbs, [0])
+        #self.assert_almost_equal(ccs, [0])
+
+        #alphas, betas, gammas = qha.get_angles(tstart=0, tstop=1000, num=7, volumes="volumes")
+        #self.assert_almost_equal(alphas, [0])
+        #self.assert_almost_equal(betas, [0])
+        #self.assert_almost_equal(gammas, [0])
+
+        #fit = qha.fit_forth(tstart=0, tstop=1000, num=1, energy="energy", volumes="volumes")
+        #self.assert_almost_equal(fit.min_vol, 0)
+        #self.assert_almost_equal(fit.min_en, 0)
+        #self.assert_almost_equal(fit.F2D_V, 0)
+
+        vols = qha.vol_E2Vib1_forth(tstart=0, tstop=1000, num=5)
+        self.assert_almost_equal(vols, [40.2380458, 40.2411873, 40.3202008, 40.4207685, 40.5276663])
+
+        vols = qha.vol_EinfVib1_forth(tstart=0, tstop=1000, num=5)
+        self.assert_almost_equal(vols, [40.2444849, 40.247707 , 40.3291996, 40.4341991, 40.5474182])
+
+        vols = qha.vol_Einf_Vib4_forth(tstart=0, tstop=1000, num=5)
+        self.assert_almost_equal(vols, [40.2456751, 40.2432144, 40.3196922, 40.4241022, 40.5411743])
+
+        alphas = qha.get_thermal_expansion_coeff_4th(tstart=0, tstop=1000, num=5, tref=None)
+        self.assert_almost_equal(alphas, [-0.0000000e+00,  4.6240312e-06,  9.4381513e-06,  1.1048456e-05, 1.2028183e-05])
+
+        free_energies = qha.get_vib_free_energies(tstart=0, tstop=100, num=2)
+        ref_energies = [
+           [0.12305556, 0.11944085],
+           [0.12140287, 0.11786901],
+           [0.11976321, 0.11629327],
+           [0.11813874, 0.11471864],
+           [0.11653109, 0.11314891]
+        ]
+        self.assert_almost_equal(free_energies, ref_energies)
+
+        data = qha.get_thermodynamic_properties(tstart=0, tstop=1000, num=2)
+        self.assert_almost_equal(data.tmesh, [   0., 1000.])
+        self.assert_almost_equal(data.cv,
+            [[0.        , 0.00050313],
+             [0.        , 0.0005035 ],
+             [0.        , 0.00050386],
+             [0.        , 0.0005042 ],
+             [0.        , 0.00050453]])
+
+        self.assert_almost_equal(data.free_energy[0], [ 0.12305556, -0.45237198])
+        self.assert_almost_equal(data.entropy[0], [0., 0.0009788])
+        self.assert_almost_equal(data.zpe, [0.1230556, 0.1214029, 0.1197632, 0.1181387, 0.1165311])
+
+        #if self.has_matplotlib():
+        if False:
+            assert qha.plot_energies(tstop=tstop, tstart=tstart, num=11, show=False)
+            # Vinet
+            assert qha.plot_vol_vs_t(tstop=tstop, tstart=tstart, num=101, show=False)
+            assert qha.plot_abc_vs_t(tstop=tstop, tstart=tstart, num=101, show=False)
+            assert qha.plot_abc_vs_t(tstop=tstop, tstart=tstart, num=101, lattice="b", show=False)
+            assert qha.plot_thermal_expansion_coeff(tstop=tstop, tstart=tstart ,num=101, show=False)
+            assert qha.plot_thermal_expansion_coeff_abc(tstop=tstop, tstart=tstart ,num=101, show=False)
+            assert qha.plot_angles_vs_t(tstop=tstop, tstart=tstart, num=101, show=False)
+            # 4th order polinomial
+            assert qha.plot_vol_vs_t_4th(tstop=tstop, tstart=tstart, num=101, show=False)
+            assert qha.plot_abc_vs_t_4th(tstop=tstop, tstart=tstart, num=101, lattice="a", show=False)
+            assert qha.plot_abc_vs_t_4th(tstop=tstop, tstart=tstart, show=False)
+            assert qha.plot_thermal_expansion_coeff_4th(tref=293, show=False)
+            assert qha.plot_thermal_expansion_coeff_abc_4th(tstop=tstop, tstart=tstart ,num=101, tref=293, show=False)
+            assert qha.plot_angles_vs_t_4th(tstop=tstop, tstart=tstart, num=101, angle=3, show=False)
diff --git a/abipy/examples/plot/plot_qha_v-ZSISA.py b/abipy/examples/plot/plot_qha_v-ZSISA.py
@@ -0,0 +1,67 @@
+#!/usr/bin/env python
+"""
+Quasi-harmonic approximation
+============================
+
+This example shows how to use the GSR.nc and PHDOS.nc files computed with different volumes
+to compute thermodynamic properties within the quasi-harmonic approximation.
+"""
+import os
+import abipy.data as abidata
+
+from abipy.dfpt.qha_aproximation import QHA_App
+
+# Root points to the directory in the git submodule with the output results.
+root = os.path.join(abidata.dirpath, "data_v-ZSISA-QHA.git", "Si_v_ZSISA_approximation")
+
+strains = [96, 98, 100, 102, 104, 106]
+strains2 = [98, 100, 102, 104, 106] # EinfVib4(D)
+#strains2 = [96, 98, 100, 102, 104] # EinfVib4(S)
+#strains2 = [100, 102, 104] # EinfVib2(D)
+
+gsr_paths = [os.path.join(root, "scale_{:d}_GSR.nc".format(s)) for s in strains]
+ddb_paths = [os.path.join(root, "scale_{:d}_GSR_DDB".format(s)) for s in strains]
+dos_paths = [os.path.join(root, "scale_{:d}_PHDOS.nc".format(s)) for s in strains2]
+
+qha = QHA_App.from_files_app_ddb(ddb_paths, dos_paths)
+tstart, tstop = 0, 800
+qha.plot_energies(tstop=tstop, tstart=tstart, num=11,
+                  title="Energies as a function of volume for different T")
+
+# Vinet
+qha.plot_vol_vs_t(tstop=tstop, tstart=tstart, num=101,
+                  title="Volume as a function of T")
+
+qha.plot_abc_vs_t(tstop=tstop, tstart=tstart, num=101,
+                  title="Lattice as a function of T")
+
+qha.plot_abc_vs_t(tstop=tstop, tstart=tstart, num=101, lattice="b",
+                  title="Lattice as a function of T")
+
+qha.plot_thermal_expansion_coeff(tstop=tstop, tstart=tstart ,num=101,
+                                 title="Volumetric thermal expansion coefficient as a function of T")
+
+qha.plot_thermal_expansion_coeff_abc(tstop=tstop, tstart=tstart ,num=101,
+                                     title="Thermal expansion coefficient as a function of T")
+
+qha.plot_angles_vs_t(tstop=tstop, tstart=tstart, num=101,
+                     title="Angles as a function of T")
+
+# 4th order polinomial
+qha.plot_vol_vs_t_4th(tstop=tstop, tstart=tstart, num=101,
+                      title="Volume as a function of T")
+
+qha.plot_abc_vs_t_4th(tstop=tstop, tstart=tstart, num=101, lattice="a",
+                      title="Lattice as a function of T")
+
+qha.plot_abc_vs_t_4th(tstop=tstop, tstart=tstart,
+                      title="Lattice as a function of T")
+
+qha.plot_thermal_expansion_coeff_4th(tref=293,
+                                     title="Volumetric thermal expansion coefficient as a function of T")
+
+qha.plot_thermal_expansion_coeff_abc_4th(tstop=tstop, tstart=tstart ,num=101, tref=293,
+                                         title="Thermal expansion coefficient as a function of T")
+
+qha.plot_angles_vs_t_4th(tstop=tstop, tstart=tstart, num=101, angle=3,
+                         title="Angles as a function of T")