From a7dc13be4a740e19a9c04ac13a000d74f6954e8b Mon Sep 17 00:00:00 2001 From: Matteo Giantomassi Date: Thu, 12 Dec 2024 14:46:47 +0100 Subject: [PATCH] Add plot_qha_v-ZSISA.py and preliminary tests in test_qha_approximation.py --- abipy/dfpt/qha_aproximation.py | 33 ++++--- abipy/dfpt/tests/test_qha_approximation.py | 103 +++++++++++++++++++++ abipy/examples/plot/plot_qha_v-ZSISA.py | 67 ++++++++++++++ 3 files changed, 191 insertions(+), 12 deletions(-) create mode 100644 abipy/dfpt/tests/test_qha_approximation.py create mode 100755 abipy/examples/plot/plot_qha_v-ZSISA.py diff --git a/abipy/dfpt/qha_aproximation.py b/abipy/dfpt/qha_aproximation.py index 859b56895..237d84c28 100644 --- 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 new file mode 100644 index 000000000..583ab1568 --- /dev/null +++ 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(' 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 new file mode 100755 index 000000000..508d6e071 --- /dev/null +++ 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")