From cab90ca7160fd0f16cd9c607569849d1c618e8b8 Mon Sep 17 00:00:00 2001 From: Matteo Giantomassi Date: Thu, 12 Dec 2024 15:51:20 +0100 Subject: [PATCH] Add test_qha_2d.py plot_qha_2d.py --- .github/workflows/gh-pages.yml | 5 + .github/workflows/test.yml | 5 + abipy/dfpt/qha_2D.py | 36 +-- abipy/dfpt/tests/test_qha_2d.py | 55 +++++ abipy/dfpt/tests/test_qha_approximation.py | 4 +- abipy/examples/plot/plot_qha_2d.py | 34 +++ docs/ghp_import.py | 248 --------------------- tasks.py | 17 ++ 8 files changed, 130 insertions(+), 274 deletions(-) create mode 100644 abipy/dfpt/tests/test_qha_2d.py create mode 100755 abipy/examples/plot/plot_qha_2d.py delete mode 100755 docs/ghp_import.py diff --git a/.github/workflows/gh-pages.yml b/.github/workflows/gh-pages.yml index 80ffb9028..ef2132417 100644 --- a/.github/workflows/gh-pages.yml +++ b/.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 diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index e958ac960..5f8ee23c1 100644 --- a/.github/workflows/test.yml +++ b/.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 diff --git a/abipy/dfpt/qha_2D.py b/abipy/dfpt/qha_2D.py index a5fd63cb4..c4e7aed8b 100644 --- a/abipy/dfpt/qha_2D.py +++ b/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 - - diff --git a/abipy/dfpt/tests/test_qha_2d.py b/abipy/dfpt/tests/test_qha_2d.py new file mode 100644 index 000000000..a1c3f524c --- /dev/null +++ b/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) diff --git a/abipy/dfpt/tests/test_qha_approximation.py b/abipy/dfpt/tests/test_qha_approximation.py index 583ab1568..59ce5092a 100644 --- a/abipy/dfpt/tests/test_qha_approximation.py +++ b/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) diff --git a/abipy/examples/plot/plot_qha_2d.py b/abipy/examples/plot/plot_qha_2d.py new file mode 100755 index 000000000..98b0b4db5 --- /dev/null +++ b/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) diff --git a/docs/ghp_import.py b/docs/ghp_import.py deleted file mode 100755 index 41f437658..000000000 --- a/docs/ghp_import.py +++ /dev/null @@ -1,248 +0,0 @@ -#! /usr/bin/env python -# -# This file is part of the ghp-import package released under -# the Tumbolia Public License. See the LICENSE file for more -# information. - -import errno -import optparse as op -import os -import subprocess as sp -import sys -import time -import unicodedata - -__usage__ = "%prog [OPTIONS] DIRECTORY" - - -if sys.version_info[0] == 3: - def enc(text): - if isinstance(text, bytes): - return text - return text.encode() - - def dec(text): - if isinstance(text, bytes): - return text.decode('utf-8') - return text - - def write(pipe, data): - try: - pipe.stdin.write(data) - except IOError as e: - if e.errno != errno.EPIPE: - raise -else: - def enc(text): - if isinstance(text, unicode): - return text.encode('utf-8') - return text - - def dec(text): - if isinstance(text, unicode): - return text - return text.decode('utf-8') - - def write(pipe, data): - pipe.stdin.write(data) - - -class Git(object): - def __init__(self, use_shell=False): - self.use_shell = use_shell - - self.cmd = None - self.pipe = None - self.stderr = None - self.stdout = None - - def check_repo(self, parser): - if self.call('rev-parse') != 0: - error = self.stderr - if not error: - error = "Unknown Git error" - error = dec(error) - if error.startswith("fatal: "): - error = error[len("fatal: "):] - parser.error(error) - - def try_rebase(self, remote, branch): - rc = self.call('rev-list', '--max-count=1', '%s/%s' % (remote, branch)) - if rc != 0: - return True - rev = dec(self.stdout.strip()) - rc = self.call('update-ref', 'refs/heads/%s' % branch, rev) - if rc != 0: - return False - return True - - def get_config(self, key): - self.call('config', key) - return self.stdout.strip() - - def get_prev_commit(self, branch): - rc = self.call('rev-list', '--max-count=1', branch, '--') - if rc != 0: - return None - return dec(self.stdout).strip() - - def open(self, *args, **kwargs): - self.cmd = ['git'] + list(args) - if sys.version_info >= (3, 2, 0): - kwargs['universal_newlines'] = False - for k in 'stdin stdout stderr'.split(): - kwargs[k] = sp.PIPE - kwargs['shell'] = self.use_shell - self.pipe = sp.Popen(self.cmd, **kwargs) - return self.pipe - - def call(self, *args, **kwargs): - self.open(*args, **kwargs) - (self.stdout, self.stderr) = self.pipe.communicate() - return self.pipe.wait() - - def check_call(self, *args, **kwargs): - kwargs["shell"] = self.use_shell - sp.check_call(['git'] + list(args), **kwargs) - - -def normalize_path(path): - # Fix unicode pathnames on OS X - # See: http://stackoverflow.com/a/5582439/44289 - if sys.platform == "darwin": - return unicodedata.normalize("NFKC", dec(path)) - return path - - -def mk_when(timestamp=None): - if timestamp is None: - timestamp = int(time.time()) - currtz = time.strftime('%z') - return "%s %s" % (timestamp, currtz) - - -def start_commit(pipe, git, branch, message): - uname = dec(git.get_config("user.name")) - email = dec(git.get_config("user.email")) - write(pipe, enc('commit refs/heads/%s\n' % branch)) - write(pipe, enc('committer %s <%s> %s\n' % (uname, email, mk_when()))) - write(pipe, enc('data %d\n%s\n' % (len(enc(message)), message))) - head = git.get_prev_commit(branch) - if head: - write(pipe, enc('from %s\n' % head)) - write(pipe, enc('deleteall\n')) - - -def add_file(pipe, srcpath, tgtpath): - with open(srcpath, "rb") as handle: - if os.access(srcpath, os.X_OK): - write(pipe, enc('M 100755 inline %s\n' % tgtpath)) - else: - write(pipe, enc('M 100644 inline %s\n' % tgtpath)) - data = handle.read() - write(pipe, enc('data %d\n' % len(data))) - write(pipe, enc(data)) - write(pipe, enc('\n')) - - -def add_nojekyll(pipe): - write(pipe, enc('M 100644 inline .nojekyll\n')) - write(pipe, enc('data 0\n')) - write(pipe, enc('\n')) - - -def add_cname(pipe, cname): - write(pipe, enc('M 100644 inline CNAME\n')) - write(pipe, enc('data %d\n%s\n' % (len(enc(cname)), cname))) - - -def gitpath(fname): - norm = os.path.normpath(fname) - return "/".join(norm.split(os.path.sep)) - - -def run_import(git, srcdir, opts): - cmd = ['git', 'fast-import', '--date-format=raw', '--quiet'] - kwargs = { - "stdin": sp.PIPE, - "shell": opts.use_shell - } - if sys.version_info >= (3, 2, 0): - kwargs["universal_newlines"] = False - pipe = sp.Popen(cmd, **kwargs) - start_commit(pipe, git, opts.branch, opts.mesg) - for path, dnames, fnames in os.walk(srcdir, followlinks=opts.followlinks): - for fn in fnames: - fpath = os.path.join(path, fn) - fpath = normalize_path(fpath) - gpath = gitpath(os.path.relpath(fpath, start=srcdir)) - add_file(pipe, fpath, gpath) - if opts.nojekyll: - add_nojekyll(pipe) - if opts.cname is not None: - add_cname(pipe, opts.cname) - write(pipe, enc('\n')) - pipe.stdin.close() - if pipe.wait() != 0: - sys.stdout.write(enc("Failed to process commit.\n")) - - -def options(): - return [ - op.make_option('-n', '--no-jekyll', dest='nojekyll', default=False, - action="store_true", - help='Include a .nojekyll file in the branch.'), - op.make_option('-c', '--cname', dest='cname', default=None, - help='Write a CNAME file with the given CNAME.'), - op.make_option('-m', '--message', dest='mesg', - default='Update documentation', - help='The commit message to use on the target branch.'), - op.make_option('-p', '--push', dest='push', default=False, - action='store_true', - help='Push the branch to origin/{branch} after committing.'), - op.make_option('-f', '--force', dest='force', - default=False, action='store_true', - help='Force the push to the repository'), - op.make_option('-r', '--remote', dest='remote', default='origin', - help='The name of the remote to push to. [%default]'), - op.make_option('-b', '--branch', dest='branch', default='gh-pages', - help='Name of the branch to write to. [%default]'), - op.make_option('-s', '--shell', dest='use_shell', default=False, - action='store_true', - help='Use the shell when invoking Git. [%default]'), - op.make_option('-l', '--follow-links', dest='followlinks', - default=False, action='store_true', - help='Follow symlinks when adding files. [%default]') - ] - - -def main(): - parser = op.OptionParser(usage=__usage__, option_list=options()) - opts, args = parser.parse_args() - - if len(args) == 0: - parser.error("No import directory specified.") - - if len(args) > 1: - parser.error("Unknown arguments specified: %s" % ', '.join(args[1:])) - - if not os.path.isdir(args[0]): - parser.error("Not a directory: %s" % args[0]) - - git = Git(use_shell=opts.use_shell) - git.check_repo(parser) - - if not git.try_rebase(opts.remote, opts.branch): - parser.error("Failed to rebase %s branch." % opts.branch) - - run_import(git, args[0], opts) - - if opts.push: - if opts.force: - git.check_call('push', opts.remote, opts.branch, '--force') - else: - git.check_call('push', opts.remote, opts.branch) - - -if __name__ == '__main__': - main() diff --git a/tasks.py b/tasks.py index dbdc31f33..7d378b987 100644 --- a/tasks.py +++ b/tasks.py @@ -15,6 +15,23 @@ DOCS_DIR = os.path.join(ABIPY_ROOTDIR, "docs") +@task +def pull(ctx): + """"Execute `git stash && git pull --recurse-submodules && git stash apply && makemake`""" + ctx.run("git stash") + ctx.run("git pull --recurse-submodules") + ctx.run("git stash apply") + + +@task +def submodules(ctx): + """Update submodules.""" + with cd(ABIPY_ROOTDIR): + # https://stackoverflow.com/questions/1030169/easy-way-to-pull-latest-of-all-git-submodules + ctx.run("git submodule update --remote --init", pty=True) + ctx.run("git submodule update --recursive --remote", pty=True) + + @task def make_doc(ctx): with cd(DOCS_DIR):