Add test_qha_2d.py plot_qha_2d.py

.github/workflows/gh-pages.yml

@@ -57,6 +57,11 @@ jobs:
             conda install abinit -c conda-forge --yes
             mpirun -n 1 abinit --version
             mpirun -n 1 abinit --build
+
+            # Update submodules with data.
+            git submodule update --remote --init
+            git submodule update --recursive --remote
+
             pip install --editable .
             mkdir -p $HOME/.abinit/abipy/
             cp abipy/data/managers/gh_manager.yml $HOME/.abinit/abipy/manager.yml

.github/workflows/test.yml

@@ -51,6 +51,11 @@ jobs:
             conda install abinit -c conda-forge --yes
             mpirun -n 1 abinit --version
             mpirun -n 1 abinit --build
+
+            # Update submodules with data.
+            git submodule update --remote --init
+            git submodule update --recursive --remote
+
             pip install -r requirements.txt
             pip install -r requirements-optional.txt
             pip install -r requirements-panel.txt

abipy/dfpt/qha_2D.py

@@ -3,20 +3,15 @@
 import numpy as np
 import abipy.core.abinit_units as abu
 
-from scipy.interpolate import UnivariateSpline
-from monty.collections import dict2namedtuple
-from monty.functools import lazy_property
-from pymatgen.analysis.eos import EOS
-from abipy.core.func1d import Function1D
+#from scipy.interpolate import UnivariateSpline
+#from monty.collections import dict2namedtuple
+#from monty.functools import lazy_property
 from abipy.tools.plotting import add_fig_kwargs, get_ax_fig_plt, get_axarray_fig_plt
 from abipy.electrons.gsr import GsrFile
 from abipy.dfpt.ddb import DdbFile
-from abipy.dfpt.phonons import PhononBandsPlotter, PhononDos, PhdosFile
-from abipy.dfpt.gruneisen import GrunsNcFile
-from scipy.interpolate import RectBivariateSpline, RegularGridInterpolator
+from abipy.dfpt.phonons import PhdosFile # PhononBandsPlotter, PhononDos,
+from scipy.interpolate import RectBivariateSpline #, RegularGridInterpolator
 
-import matplotlib.pyplot as plt
-from abipy.tools.plotting import get_ax_fig_plt
 
 class QHA_2D:
     """
@@ -56,14 +51,14 @@ def __init__(self, structures, doses, energies, structures_from_phdos, eos_name=
         # Find index of minimum energy
         self.min_energy_idx = np.unravel_index(np.nanargmin(self.energies), self.energies.shape)
 
-    def plot_energies(self, ax=None):
+    @add_fig_kwargs
+    def plot_energies(self, ax=None, **kwargs):
         """
         Plot energy surface and visualize minima in a 3D plot.
 
         Args:
             ax: Matplotlib axis for the plot. If None, creates a new figure.
         """
-
         a0 = self.lattice_a[:, 0]
         c0 = self.lattice_c[0, :]
 
@@ -75,7 +70,6 @@ def plot_energies(self, ax=None):
 
         X, Y = np.meshgrid(c0,a0)
 
-
         # Plot the surface
         ax.plot_wireframe(X, Y, self.energies, cmap='viridis')
         ax.scatter(self.lattice_c[0,self.iy0], self.lattice_a[self.ix0,0], self.energies[self.ix0, self.iy0], color='red', s=100)
@@ -100,7 +94,6 @@ def plot_energies(self, ax=None):
         ax.set_ylabel('Lattice Parameter A (Å)')
         ax.set_zlabel('Energy (eV)')
         ax.set_title('Energy Surface in 3D')
-        plt.show()
 
         return fig
 
@@ -133,6 +126,7 @@ def find_minimum(self, f_interp, xy_init, tol=1e-6, max_iter=1000, step_size=0.0
         min_energy = f_interp(xy[0], xy[1])
         return xy[0], xy[1], min_energy
 
+    @add_fig_kwargs
     def plot_free_energies(self, tstart=800 , tstop=0 ,num=5, ax=None, **kwargs):
         """
         Plot free energy as a function of temperature in a 3D plot.
@@ -224,11 +218,10 @@ def plot_free_energies(self, tstart=800 , tstop=0 ,num=5, ax=None, **kwargs):
         ax.set_zlabel('Energy (eV)')
         #ax.set_title('Energies as a 3D Plot')
         plt.savefig("energy.pdf", format="pdf", bbox_inches="tight")
-        plt.show()  # Display the plot
-
 
         return fig
 
+    @add_fig_kwargs
     def plot_thermal_expansion(self, tstart=800, tstop=0, num=81, ax=None, **kwargs):
         """
         Plots thermal expansion coefficients along the a-axis, c-axis, and volumetric alpha.
@@ -327,16 +320,15 @@ def plot_thermal_expansion(self, tstart=800, tstop=0, num=81, ax=None, **kwargs)
         file_path = 'thermal-expansion_data.txt'
         np.savetxt(file_path, data_to_save, fmt='%4.6e', delimiter='\t\t',  header='\t\t\t'.join(columns), comments='')
 
-
         ax.grid(True)
         ax.legend()
         ax.set_xlabel('Temperature (K)')
         ax.set_ylabel(r'Thermal Expansion Coefficients ($\alpha$)')
         plt.savefig("thermal_expansion.pdf", format="pdf", bbox_inches="tight")
-        plt.show()
 
         return fig
 
+    @add_fig_kwargs
     def plot_lattice(self, tstart=800, tstop=0, num=81, ax=None, **kwargs):
         """
         Plots thermal expansion coefficients along the a-axis, c-axis, and volumetric alpha.
@@ -353,6 +345,7 @@ def plot_lattice(self, tstart=800, tstop=0, num=81, ax=None, **kwargs):
         ph_energies = self.get_vib_free_energies(tstart, tstop, num)
         min_x, min_y = np.zeros(num), np.zeros(num)
         min_tot_energy = np.zeros(num)
+        import matplotlib.pyplot as plt
         fig, axs = plt.subplots(1, 3, figsize=(18, 6), sharex=True)
 
         if (len(self.lattice_a_from_phdos)==len(self.lattice_a) or  len(self.lattice_c_from_phdos)==len(self.lattice_c)):
@@ -417,12 +410,11 @@ def plot_lattice(self, tstart=800, tstop=0, num=81, ax=None, **kwargs):
             min_volumes = min_x**2 * min_y * scale
 
             axs[0].plot(tmesh, min_x, color='c', label=r"$a$ (E$\infty$Vib2)", linewidth=2)
-            axs[1].set_title("Plots of a, c, and V (E$\infty$Vib2)")
+            axs[1].set_title(r"Plots of a, c, and V (E$\infty$Vib2)")
             axs[1].plot(tmesh, min_y, color='r', label=r"$c$ (E$\infty$Vib2)", linewidth=2)
             axs[2].plot(tmesh, min_volumes, color='b', label=r"$V$ (E$\infty$Vib2)", linewidth=2)
 
 
-
         axs[0].set_ylabel("a")
         axs[0].legend()
         axs[0].grid(True)
@@ -438,8 +430,6 @@ def plot_lattice(self, tstart=800, tstop=0, num=81, ax=None, **kwargs):
 
         # Adjust layout and show the figure
         plt.tight_layout()
-        plt.show()
-
 
         return fig
 
@@ -534,5 +524,3 @@ def get_vib_free_energies(self, tstart, tstop, num):
                 if dos is not None:
                     f[j][i] = dos.get_free_energy(tstart, tstop, num).values
         return f
-
-

abipy/dfpt/tests/test_qha_2d.py

@@ -0,0 +1,55 @@
+"""Tests for QHA_2D"""
+import os
+import abipy.data as abidata
+
+from abipy.dfpt.qha_2D import QHA_2D
+from abipy.core.testing import AbipyTest
+
+
+class Qha2dTest(AbipyTest):
+
+    def test_qha_2d(self):
+
+        strains_a = [ 995,1000, 1005, 1010, 1015  ]
+        strains_c = [ 995,1000, 1005, 1010, 1015  ]
+        strains_a1 = [ 1000, 1005, 1010 ]
+        strains_c1 = [ 1000, 1005, 1010 ]
+
+        # Root points to the directory in the git submodule with the output results.
+        root = os.path.join(abidata.dirpath, "data_v-ZSISA-QHA.git", "ZnO_ZSISA_approximation")
+
+        gsr_paths = [[os.path.join(root, f"scale_{s1}_{s3}/out_GSR_DDB") for s3 in strains_c] for s1 in strains_a]
+        dos_paths = [[os.path.join(root, f"scale_{s1}_{s3}/out_PHDOS.nc") for s3 in strains_c1] for s1 in strains_a1]
+
+        qha = QHA_2D.from_files(gsr_paths, dos_paths, gsr_file="DDB")
+
+        # Test properties and methods.
+        assert qha.pressure == 0
+
+        f_mat = qha.get_vib_free_energies(0, 1000, 2)
+        ref_mat = [
+            [
+                [0.23190045769242368, -1.0257329004129736],
+                [0.22879089684005882, -1.0344811855487408],
+                [0.0, 0.0],
+            ],
+            [
+                [0.23044275786150215, -1.0297189635203232],
+                [0.22735644425413698, -1.038628734783889],
+                [0.22427217946933886, -1.0480192052628534],
+            ],
+            [
+                [0.0, 0.0],
+                [0.2259076206066778, -1.0429439417530417],
+                [0.0, 0.0],
+            ],
+        ]
+
+        self.assert_almost_equal(f_mat, ref_mat)
+
+        if self.has_matplotlib():
+        #if False:
+            qha.plot_energies(show=False)
+            qha.plot_free_energies(tstop=500, tstart=0, num=6, show=False)
+            qha.plot_thermal_expansion(tstop=1000, tstart=0, num=101, show=False)
+            qha.plot_lattice(tstop=1000, tstart=0, num=101, show=False)

abipy/dfpt/tests/test_qha_approximation.py

@@ -84,8 +84,8 @@ def test_v_ZSIZA(self):
         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:
+        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)

abipy/examples/plot/plot_qha_2d.py

@@ -0,0 +1,34 @@
+#!/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_2D import QHA_2D
+
+strains_a = [ 995,1000, 1005, 1010, 1015  ]
+strains_c = [ 995,1000, 1005, 1010, 1015  ]
+strains_a1 = [ 1000, 1005, 1010 ]
+strains_c1 = [ 1000, 1005, 1010 ]
+
+# Root points to the directory in the git submodule with the output results.
+root = os.path.join(abidata.dirpath, "data_v-ZSISA-QHA.git", "ZnO_ZSISA_approximation")
+
+#gsr_paths = [[f"scale_{s1}_{s3}/out_GSR.nc" for s3 in strains_c] for s1 in strains_a]
+gsr_paths = [[os.path.join(root, f"scale_{s1}_{s3}/out_GSR_DDB") for s3 in strains_c] for s1 in strains_a]
+dos_paths = [[os.path.join(root, f"scale_{s1}_{s3}/out_PHDOS.nc") for s3 in strains_c1] for s1 in strains_a1]
+
+qha = QHA_2D.from_files(gsr_paths, dos_paths, gsr_file="DDB")
+
+qha.plot_energies()
+
+qha.plot_free_energies(tstop=500, tstart=0, num=6)
+
+qha.plot_thermal_expansion(tstop=1000, tstart=0, num=101)
+
+qha.plot_lattice(tstop=1000, tstart=0, num=101)