diff --git a/abipy/dfpt/deformation_utils.py b/abipy/dfpt/deformation_utils.py index 0561d76b0..b0521fdc0 100644 --- a/abipy/dfpt/deformation_utils.py +++ b/abipy/dfpt/deformation_utils.py @@ -51,7 +51,7 @@ def _add(name, new_rprim, i, j, k, l, m, n) -> None: structures_new[name] = new_structure strain_inds.append([i, j, k, l, m, n]) - i, j, k, l, m, n = 6 * [None] + i, j, k, l, m, n = 6 * [0] if 1 <= spgrp_number <= 2: disp=[[1,1,1,1,1,1], [0,1,1,1,1,1], [2,1,1,1,1,1], [1,0,1,1,1,1], [1,2,1,1,1,1], [1,1,0,1,1,1], @@ -193,5 +193,5 @@ def _add(name, new_rprim, i, j, k, l, m, n) -> None: else: raise ValueError(f"Invalid {spgrp_number=}") - return structures_new, np.array(strain_inds, dtype=object), spgrp_number + return structures_new, np.array(strain_inds, dtype=int), spgrp_number diff --git a/abipy/dfpt/qha_2D.py b/abipy/dfpt/qha_2D.py index 25bcf369a..054e6ed00 100644 --- a/abipy/dfpt/qha_2D.py +++ b/abipy/dfpt/qha_2D.py @@ -12,44 +12,96 @@ import numpy as np import abipy.core.abinit_units as abu +from scipy.interpolate import RectBivariateSpline #, RegularGridInterpolator #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.tools.typing import Figure +from abipy.tools.serialization import HasPickleIO, mjson_load from abipy.electrons.gsr import GsrFile from abipy.dfpt.ddb import DdbFile from abipy.dfpt.phonons import PhdosFile # PhononBandsPlotter, PhononDos, -from scipy.interpolate import RectBivariateSpline #, RegularGridInterpolator +from abipy.dfpt.vzsisa import anaget_phdoses_with_gauss + -class QHA_2D: +class QHA_2D(HasPickleIO): """ Quasi-Harmonic Approximation (QHA) analysis in 2D. Provides methods for calculating and visualizing energy, free energy, and thermal expansion. """ - #@classmethod - #def from_ddb_files(cls, ddb_files): - # for row in ddb_files: - # for gp in row: - # if os.path.exists(gp): + @classmethod + def from_json_file(cls, + filepath: PathLike, + nqsmall_or_qppa: int, + anaget_kwargs: dict | None = None, + smearing_ev: float | None = None, + verbose: int = 0) -> Vzsisa: + """ + Build an instance from a json file `filepath` typically produced by an Abipy flow. + For the meaning of the other arguments see from_gsr_ddb_paths. + """ + data = mjson_load(filepath) + + bo_strains_ac = [data["strains_a"], data["strains_c"]] + phdos_strains_ac = [data["strains_a"], data["strains_c"]] - # return cls.from_files(ddb_files, phdos_paths_2D, gsr_file="DDB") + return cls.from_gsr_ddb_paths(nqsmall_or_qppa, + data["gsr_relax_paths"], data["ddb_relax_paths"], + bo_strains_ac, phdos_strains_ac, + anaget_kwargs=anaget_kwargs, smearing_ev=smearing_ev, verbose=verbose) @classmethod - def from_files(cls, gsr_paths_2D, phdos_paths_2D, gsr_file="GSR.nc") -> QHA_2D: + def from_gsr_ddb_paths(cls, + nqsmall_or_qppa: int, + gsr_paths, + ddb_paths, + bo_strains_ac, + phdos_strains_ac, + anaget_kwargs: dict | None = None, + smearing_ev: float | None = None, + verbose: int = 0) -> QHA_2D: + """ + Creates an instance from a list of GSR files and a list of DDB files. + This is a simplified interface that computes the PHDOS.nc files automatically + from the DDB files by invoking anaddb + + Args: + nqsmall_or_qppa: Define the q-mesh for the computation of the PHDOS. + if > 0, it is interpreted as nqsmall + if < 0, it is interpreted as qppa. + gsr_paths: list of paths to GSR files. + ddb_paths: list of paths to DDB files. + bo_strains_ac: List of strains for the a and the c lattice vector. + phdos_strains_ac: List of strains for the a and the c lattice vector. + anaget_kwargs: dict with extra arguments passed to anaget_phdoses_with_gauss. + smearing_ev: Gaussian smearing in eV. + verbose: Verbosity level. + """ + phdos_paths, phbands_paths = anaget_phdoses_with_gauss(nqsmall_or_qppa, smearing_ev, ddb_paths, anaget_kwargs, verbose) + + new = cls.from_files(ddb_files, phdos_paths_2D, bo_strains_ac, phdos_strains_ac, gsr_file="GSR.nc") + #new.pickle_dump(workdir, basename=None) + return new + + @classmethod + def from_files(cls, gsr_paths_2D, phdos_paths_2D, bo_strains_ac, phdos_strains_ac, gsr_file="GSR.nc") -> QHA_2D: """ Creates an instance of QHA from 2D lists of GSR and PHDOS files. Args: gsr_paths_2D: 2D list of paths to GSR files. phdos_paths_2D: 2D list of paths to PHDOS.nc files. - - Returns: - A new instance of QHA. + bo_strains_ac: List of strains for the a and the c lattice vector. + phdos_strains_ac: List of strains for the a and the c lattice vector. """ energies, structures, phdoses , structures_from_phdos = [], [], [],[] + #shape = (len(strains_a), len(strains_c)) + #gsr_paths_2d = np.reshape(gsr_paths_2D, shape) + #phdos_paths_2d = np.reshape(phdos_paths_2D, shape) + if gsr_file == "GSR.nc": # Process GSR files for row in gsr_paths_2D: @@ -98,15 +150,18 @@ def from_files(cls, gsr_paths_2D, phdos_paths_2D, gsr_file="GSR.nc") -> QHA_2D: phdoses.append(row_doses) structures_from_phdos.append(row_structures) - return cls(structures, phdoses, energies , structures_from_phdos) + return cls(structures, phdoses, energies, structures_from_phdos, bo_strains_ac, phdos_strains_ac) def __init__(self, structures, phdoses, energies, structures_from_phdos, + bo_strains_ac, phdos_strains_ac, eos_name: str='vinet', pressure: float=0.0): """ Args: structures (list): List of structures at different volumes. phdoses: List of density of states (DOS) data for phonon calculations. energies (list): SCF energies for the structures in eV. + bo_strains_ac: List of strains for the a and the c lattice vector. + phdos_strains_ac: List of strains for the a and the c lattice vector. eos_name (str): Expression used to fit the energies (e.g., 'vinet'). pressure (float): External pressure in GPa to include in p*V term. """ @@ -114,6 +169,10 @@ def __init__(self, structures, phdoses, energies, structures_from_phdos, self.structures = structures self.structures_from_phdos = structures_from_phdos self.energies = np.array(energies, dtype=np.float64) + + self.bo_strains_ac = bo_strains_ac + self.phdos_strains_ac = phdos_strains_ac + self.eos_name = eos_name self.pressure = pressure self.volumes = np.array([[s.volume if s else np.nan for s in row] for row in structures]) @@ -509,7 +568,7 @@ def get_vib_free_energies(self, tstart: float, tstop: float, num: int) -> np.nda tstop: Stop temperature. num: Number of temperature points. - Return: A 3D array of vibrational free energies. + Return: A 3D array of vibrational free energies of shape (num_c, num_a, num_temp) """ f = np.zeros((len(self.lattice_c_from_phdos[0]), len(self.lattice_a_from_phdos[:, 0]), num)) for i in range(len(self.lattice_a_from_phdos[:, 0])): diff --git a/abipy/dfpt/qha_general_stress.py b/abipy/dfpt/qha_general_stress.py index 50f70847f..2dc208999 100644 --- a/abipy/dfpt/qha_general_stress.py +++ b/abipy/dfpt/qha_general_stress.py @@ -46,20 +46,21 @@ def from_json_file(cls, verbose: Verbosity level """ data = mjson_load(filepath) - ddb_paths = data["ddb_relax_paths"] - gsr_path = data["gsr_relax_paths"] + ddb_relax_paths = data["ddb_relax_paths"] + gsr_relax_paths = data["gsr_relax_paths"] phdos_paths, phbands_paths = anaget_phdoses_with_gauss(nqsmall_or_qppa, smearing_ev, - ddb_paths, anaget_kwargs, verbose) + ddb_relax_paths, anaget_kwargs, verbose) # Create a 6D array with shape (3, 3, 3, 3, 3, 3) initialized with None. gsr_paths_6d = np.full((3, 3, 3, 3, 3, 3), None, dtype=object) phdos_paths_6d = np.full((3, 3, 3, 3, 3, 3), None, dtype=object) - cell6_inds = np.array(data["cell6_inds"], dtype=np.int) - for gsr_path, phdos_path, cell_ind in zip(gsr_paths, phdos_paths, cell6_inds, strict=True): - gsr_paths_6d[cell_ind] = gsr_path - phdos_paths_6d[cell_ind] = phdos_path + strain_inds = data["strain_inds"] + #print(strain_inds) + for gsr_path, phdos_path, inds in zip(gsr_relax_paths, phdos_paths, strain_inds, strict=True): + gsr_paths_6d[inds] = gsr_path + phdos_paths_6d[inds] = phdos_path new = cls.from_files(gsr_paths_6d, phdos_paths_6d, gsr_guess, model='ZSISA') #new.pickle_dump(workdir, basename=None) diff --git a/abipy/dfpt/tests/test_deformation_utils.py b/abipy/dfpt/tests/test_deformation_utils.py new file mode 100644 index 000000000..2e12dc8ab --- /dev/null +++ b/abipy/dfpt/tests/test_deformation_utils.py @@ -0,0 +1,20 @@ +"""Tests for deformation_utils module""" +import numpy as np +import abipy.data as abidata + +from abipy.core.testing import AbipyTest +from abipy.dfpt.deformation_utils import generate_deformations + + +class DeformationUtilsTest(AbipyTest): + + def test_generate_deformations(self): + """Testing generate_deformations""" + eps = 0.005 + si = self.get_structure("Si") + structures_dict, strain_inds, spgrp_number = generate_deformations(si, eps) + assert len(structures_dict) == 3 + assert spgrp_number == 227 + assert strain_inds.shape == (3, 6) + self.assert_equal(strain_inds[:, 0], [0, 1, 2]) + assert np.all(strain_inds[1:, 1] == 0) diff --git a/abipy/dfpt/tests/test_qha_2d.py b/abipy/dfpt/tests/test_qha_2d.py index a1c3f524c..f7944270b 100644 --- a/abipy/dfpt/tests/test_qha_2d.py +++ b/abipy/dfpt/tests/test_qha_2d.py @@ -1,5 +1,6 @@ """Tests for QHA_2D""" import os +import numpy as np import abipy.data as abidata from abipy.dfpt.qha_2D import QHA_2D @@ -8,20 +9,25 @@ class Qha2dTest(AbipyTest): - def test_qha_2d(self): + def test_zsisa_approximation(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 ] + bo_strains_a = [995, 1000, 1005, 1010, 1015] + bo_strains_c = [995, 1000, 1005, 1010, 1015] + phdos_strains_a = [1000, 1005, 1010] + phdos_strains_c = [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] + gsr_paths = [[os.path.join(root, f"scale_{s1}_{s3}/out_GSR_DDB") for s3 in bo_strains_c] for s1 in bo_strains_a] + dos_paths = [[os.path.join(root, f"scale_{s1}_{s3}/out_PHDOS.nc") for s3 in phdos_strains_c] for s1 in phdos_strains_a] + + bo_strains_ac = [bo_strains_a, bo_strains_c] + bo_strains_ac = (np.array(bo_strains_ac) - 1000) / 100 + phdos_strains_ac = [phdos_strains_a, phdos_strains_c] + phdos_strains_ac = (np.array(phdos_strains_ac) - 1000) / 100 - qha = QHA_2D.from_files(gsr_paths, dos_paths, gsr_file="DDB") + qha = QHA_2D.from_files(gsr_paths, dos_paths, bo_strains_ac, phdos_strains_ac, gsr_file="DDB") # Test properties and methods. assert qha.pressure == 0 @@ -48,8 +54,78 @@ def test_qha_2d(self): 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) + + def test_qha_2d(self): + + bo_strains_a = [995, 1000, 1005, 1010, 1015] + bo_strains_c = [995, 1000, 1005, 1010, 1015] + #bo_strains_a = [1000, 1005, 1010, 1015 , 1020] + #bo_strains_c = [1000, 1005, 1010, 1015 , 1020] + + # 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_QHA") + + #gsr_paths = [[f"scale_{s1}_{s3}/out_GSR.nc" for s3 in bo_strains_c] for s1 in bo_strains_a] + gsr_paths = [[os.path.join(root, f"scale_{s1}_{s3}/out_GSR_DDB") for s3 in bo_strains_c] for s1 in bo_strains_a] + phdos_paths = [[os.path.join(root, f"scale_{s1}_{s3}/out_PHDOS.nc") for s3 in bo_strains_c] for s1 in bo_strains_a] + + bo_strains_ac = [bo_strains_a, bo_strains_c] + bo_strains_ac = (np.array(bo_strains_ac) - 1000) / 100 + + #qha = QHA_2D.from_files(gsr_paths, phdos_paths, bo_strains_ac, phdos_strains_ac, gsr_file="GSR.nc") + qha = QHA_2D.from_files(gsr_paths, phdos_paths, bo_strains_ac, bo_strains_ac, gsr_file="DDB") + + # Test properties and methods. + assert qha.pressure == 0 + + f_mat = qha.get_vib_free_energies(0, 1000, 2) + #print(f_mat) + ref_mat = [ + [ + [ 0.2364734 , -1.0138777 ], + [ 0.23334375, -1.02191455], + [ 0.23021055, -1.03048734], + [ 0.22707712, -1.03958781], + [ 0.22394442, -1.04917831] + ], + [ + [ 0.23500674, -1.0175096 ], + [ 0.23190046, -1.0257329 ], + [ 0.2287909 , -1.03448119], + [ 0.22568001, -1.0437223 ], + [ 0.22257111, -1.05344296] + ], + [ + [ 0.23352563, -1.02131569], + [ 0.23044276, -1.02971896], + [ 0.22735644, -1.03862873], + [ 0.22427218, -1.04801921], + [ 0.2211813 , -1.05787544] + ], + [ + [ 0.23203085, -1.02529729], + [ 0.22897121, -1.0338748 ], + [ 0.22590762, -1.04294394], + [ 0.22284219, -1.05248871], + [ 0.21977746, -1.06247964] + ], + [ + [ 0.2305236 , -1.02944732], + [ 0.22748688, -1.03820052], + [ 0.22444611, -1.04742487], + [ 0.22140302, -1.05711238], + [ 0.21836046, -1.06723406] + ] + ] + + self.assert_almost_equal(f_mat, ref_mat) + + if self.has_matplotlib(): + assert qha.plot_energies(show=False) + assert qha.plot_free_energies(tstop=500, tstart=0, num=6, show=False) + assert qha.plot_thermal_expansion(tstop=1000, tstart=0, num=101, show=False) + assert qha.plot_lattice(tstop=1000, tstart=00, num=101, show=False) diff --git a/abipy/dfpt/vzsisa.py b/abipy/dfpt/vzsisa.py index 1c288959a..d0bfdf5f4 100644 --- a/abipy/dfpt/vzsisa.py +++ b/abipy/dfpt/vzsisa.py @@ -119,9 +119,9 @@ def from_gsr_ddb_paths(cls, if < 0, it is interpreted as qppa. gsr_paths: list of paths to GSR files. ddb_paths: list of paths to DDB files. - anaget_kwargs: dict with extra arguments passed to anaget_phdoses_with_gauss - smearing_ev: Gaussian smearing in eV - verbose: + anaget_kwargs: dict with extra arguments passed to anaget_phdoses_with_gauss. + smearing_ev: Gaussian smearing in eV. + verbose: Verbosity level. """ phdos_paths, phbands_paths = anaget_phdoses_with_gauss(nqsmall_or_qppa, smearing_ev, ddb_paths, anaget_kwargs, verbose) new = cls.from_gsr_phdos_files(gsr_paths, phdos_paths) diff --git a/abipy/examples/flows/run_qha_2d.py b/abipy/examples/flows/run_qha_2d.py new file mode 100755 index 000000000..b4168762c --- /dev/null +++ b/abipy/examples/flows/run_qha_2d.py @@ -0,0 +1,117 @@ +#!/usr/bin/env python +r""" +Flow for quasi-harmonic calculations under development +===================================================== +Warning: This code is still under development. +""" +import sys +import os +import numpy as np +import abipy.abilab as abilab +import abipy.data as abidata + +from abipy import flowtk +from abipy.flowtk.qha_2d import Qha2dFlow + + +def build_flow(options): + """ + Create a `QhaFlow` for quasi-harmonic calculations. + """ + # Working directory (default is the name of the script with '.py' removed and "run_" replaced by "flow_") + if not options.workdir: + __file__ = os.path.join(os.getcwd(), "run_qha_2d.py") + options.workdir = os.path.basename(__file__).replace(".py", "").replace("run_", "flow_") + + # Initialize structure and pseudos for ZnO. + structure = abilab.Structure.from_abistring(""" +natom 4 +ntypat 2 +typat +1 1 2 +2 +znucl 30 8 +xred + 0.0000000000 0.0000000000 -0.0025137620 + 0.3333333333 0.6666666667 0.4974862380 + 0.0000000000 0.0000000000 0.3835201241 + 0.3333333333 0.6666666667 0.8835201241 +acell 1.0 1.0 1.0 +rprim + 6.3016720624 0.0000000000 0.0000000000 + -3.1508360312 5.4574080923 0.0000000000 + 0.0000000000 0.0000000000 9.7234377918 +""") + + # Use NC PBE pseudos from pseudodojo v0.4 + from abipy.flowtk.psrepos import get_oncvpsp_pseudos + pseudos = get_oncvpsp_pseudos(xc_name="PBE", version="0.4", accuracy="standard").get_pseudos_for_structure(structure) + + #ecut = max(p.hint_for_accuracy("normal").ecut for p in pseudos) + #print(ecut) + + # Select k-mesh for electrons and q-mesh for phonons. + #ngkpt = [6, 6, 4]; ngqpt = [1, 1, 1] + ngkpt = [2, 2, 2]; ngqpt = [1, 1, 1] + + with_becs = True + with_quad = False + #with_quad = not structure.has_zero_dynamical_quadrupoles + + scf_input = abilab.AbinitInput(structure, pseudos) + #scf_input.set_cutoffs_for_accuracy("standard") + + # Set other important variables (consistent with tutorial) + # All the other DFPT runs will inherit these parameters. + scf_input.set_vars( + nline=10, + nbdbuf=0, + nstep=100, + ecut=42.0, + ecutsm=1.0, + occopt=1, + nband=26, + #tolvrs=1.0e-18, # SCF stopping criterion (modify default) + tolvrs=1.0e-6, # SCF stopping criterion (modify default) + ) + + #scf_input.set_scf_nband_semicond() + scf_input.set_kmesh(ngkpt=ngkpt, shiftk=[0, 0, 0]) + + bo_strains_a = [-5, 0, 5, 10, 15] + bo_strains_c = [-5, 0, 5, 10, 15] + #bo_strains_a = [0, 5, 10, 15, 20] + #bo_strains_c = [0, 5, 10, 15, 20] + bo_strains_a = [0, 5] + bo_strains_c = [0, 5] + bo_strains_a = np.array(bo_strains_a) / 100 + bo_strains_c = np.array(bo_strains_c) / 100 + + bo_strains_ac = [bo_strains_a, bo_strains_c] + phdos_strains_ac = bo_strains_ac + + return Qha2dFlow.from_scf_input(options.workdir, scf_input, bo_strains_ac, phdos_strains_ac, ngqpt, + with_becs, with_quad, edos_ngkpt=None) + + +# This block generates the thumbnails in the Abipy gallery. +# You can safely REMOVE this part if you are using this script for production runs. +if os.getenv("READTHEDOCS", False): + __name__ = None + import tempfile + options = flowtk.build_flow_main_parser().parse_args(["-w", tempfile.mkdtemp()]) + build_flow(options).graphviz_imshow() + + +@flowtk.flow_main +def main(options): + """ + This is our main function that will be invoked by the script. + flow_main is a decorator implementing the command line interface. + Command line args are stored in `options`. + """ + return build_flow(options) + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/abipy/flowtk/qha_2d.py b/abipy/flowtk/qha_2d.py new file mode 100644 index 000000000..c68cee192 --- /dev/null +++ b/abipy/flowtk/qha_2d.py @@ -0,0 +1,218 @@ +# coding: utf-8 +""" +Workflows for calculations within the quasi-harmonic approximation. +""" +from __future__ import annotations + +import numpy as np + +from abipy.tools.serialization import mjson_write +from abipy.dfpt.deformation_utils import generate_deformations +from abipy.abio.inputs import AbinitInput +from abipy.flowtk.works import Work, PhononWork +from abipy.flowtk.tasks import RelaxTask +from abipy.flowtk.flows import Flow + + +class Qha2dFlow(Flow): + """ + Flow for QHA calculations with two degrees of freedom. + + .. rubric:: Inheritance Diagram + .. inheritance-diagram:: Qha2dFlow + """ + + @classmethod + def from_scf_input(cls, + workdir: PathLike, + scf_input: AbinitInput, + bo_strains_ac, + phdos_strains_ac, + ngqpt, + with_becs: bool, + with_quad: bool, + edos_ngkpt=None, + manager=None) -> Qha2dFlow: + """ + Build a flow for QHA calculations from an |AbinitInput| for GS-SCF calculation. + + Args: + workdir: Working directory of the flow. + scf_input: |AbinitInput| for GS-SCF run used as template to generate the other inputs. + bo_strains_ac + phdos_strains_ac + ngqpt: Three integers defining the q-mesh for phonon calculation. + with_becs: Activate calculation of Electric field and Born effective charges. + with_quad: Activate calculation of dynamical quadrupoles. Require `with_becs` + Note that only selected features are compatible with dynamical quadrupoles. + Please consult + edos_ngkpt: Three integers defining the the k-sampling for the computation of the + electron DOS with the relaxed structures. Useful for metals or small gap semiconductors + in which the electronic contribution should be included. + None disables the computation of the e-DOS. + manager: |TaskManager| instance. Use default if None. + """ + flow = cls(workdir=workdir, manager=manager) + flow.register_work(Qha2dWork.from_scf_input(scf_input, bo_strains_ac, phdos_strains_ac, + ngqpt, with_becs, with_quad, + ionmov=2, edos_ngkpt=edos_ngkpt)) + return flow + + def finalize(self): + """ + This method is called when the flow is completed. + It performs some basic post-processing of the results to facilitate further analysis. + """ + work = self[0] + data = {"bo_strains_ac": work.bo_strains_ac, "phdos_strains_ac": work.phdos_strains_ac} + + # Build list of strings with path to the relevant output files ordered by V. + data["gsr_relax_paths"] = [task.gsr_path for task in work.relax_tasks_deformed] + + gsr_relax_entries, gsr_relax_volumes = [], [] + for task in work.relax_tasks_deformed: + with task.open_gsr() as gsr: + gsr_relax_entries.append(dict( + volume=gsr.structure.volume, + energy_eV=float(gsr.energy), + pressure_GPa=float(gsr.pressure), + #structure=gsr.structure, + )) + gsr_relax_volumes.append(gsr.structure.volume) + + data["gsr_relax_entries"] = gsr_relax_entries + + data["ddb_relax_paths"] = [ph_work.outdir.has_abiext("DDB") for ph_work in work.ph_works] + data["gsr_relax_edos_paths"] = [] if not work.edos_work else [task.gsr_path for task in work.edos_work] + + # Write json file + mjson_write(data, self.outdir.path_in("qha_2d.json"), indent=4) + + return super().finalize() + + +class Qha2dWork(Work): + """ + This work performs a structural relaxation of the initial structure, then a set of distorted + structures is genenerated and the relaxed structures are used + to compute phonons, BECS and the dielectric tensor with DFPT. + + .. rubric:: Inheritance Diagram + .. inheritance-diagram:: Qha2dWork + """ + + @classmethod + def from_scf_input(cls, + scf_input: AbinitInput, + bo_strains_ac, + phdos_strains_ac, + ngqpt, + with_becs: bool, + with_quad: bool, + ionmov: int, + edos_ngkpt=None) -> Qha2dWork: + """ + Build the work from an |AbinitInput| representing a GS-SCF calculation. + See Qha2dFlow for the meaning of the arguments. + """ + work = cls() + + # Save attributes in work + work.initial_scf_input = scf_input + + work.bo_strains_ac = bo_strains_ac + work.phdos_strains_ac = phdos_strains_ac + # Make sure ench row is a numpy array. + for i in range(2): + work.bo_strains_ac[i] = np.array(bo_strains_ac[i]) + work.phdos_strains_ac[i] = np.array(phdos_strains_ac[i]) + + work.ngqpt = ngqpt + work.with_becs = with_becs + work.with_quad = with_quad + work.edos_ngkpt = edos_ngkpt if edos_ngkpt is None else np.reshape(edos_ngkpt, (3,)) + + # Create input for relaxation and register the relaxation task. + work.relax_template = relax_template = scf_input.deepcopy() + relax_template.pop_tolerances() + # optcell = 3: constant-volume optimization of cell geometry + relax_template.set_vars(optcell=3, ionmov=ionmov, tolvrs=1e-8, toldff=1.e-6) + #if optcell is not None and optcell != 0: + #relax_template.set_vars_ifnotin(ecutsm=1.0, dilatmx=1.05) + + work.initial_relax_task = work.register_relax_task(relax_template) + + return work + + def on_ok(self, sender): + """ + This method is called when one task reaches status `S_OK`. + It executes on_all_ok when all tasks in self have reached `S_OK`. + """ + if sender == self.initial_relax_task: + # Get relaxed structure and build new task for structural relaxation at fixed volume. + relaxed_structure = sender.get_final_structure() + + relax_template = self.relax_template + self.relax_tasks_deformed = [] + + import itertools + for s1, s3 in itertools.product(self.bo_strains_ac[0], self.bo_strains_ac[1]): + deformation_name = f"{s1=}, {s3=}" + # Apply strain to the structure + strain_tensor = np.diag([s1, s1, s3]) + strained_structure = relaxed_structure.apply_strain(strain_tensor, inplace=False) + #print("strained_structure:", strained_structure) + + # Relax deformed structure with fixed unit cell. + task = self.register_relax_task(relax_template.new_with_structure(strained_structure, optcell=0)) + task.bo_strain = (s1, s3) + task.in_phdos_strains = np.any(np.abs(s1 - self.phdos_strains_ac[0]) < 1e-3) and \ + np.any(np.abs(s3 - self.phdos_strains_ac[1]) < 1e-3) + self.relax_tasks_deformed.append(task) + + self.flow.allocate(build=True) + + return super().on_ok(sender) + + def on_all_ok(self): + """ + This callback is called when all tasks in the Work reach status `S_OK`. + Here we add a new PhononWork for each volume using the relaxed structure. + """ + self.edos_work = Work() + + # Build phonon works for the different relaxed structures. + self.ph_works = [] + + for task in self.relax_tasks_deformed: + s1, s3 = task.bo_strain + deformation_name = f"{s1=}, {s3=}" + + relaxed_structure = task.get_final_structure() + scf_input = self.initial_scf_input.new_with_structure(relaxed_structure) + + if task.in_phdos_strains: + ph_work = PhononWork.from_scf_input(scf_input, self.ngqpt, is_ngqpt=True, tolerance=None, + with_becs=self.with_becs, ddk_tolerance=None) + + # Reduce the number of files produced in the DFPT tasks to avoid possible disk quota issues. + prtvars = dict(prtden=0, prtpot=0) + for task in ph_work[1:]: + task.input.set_vars(**prtvars) + + ph_work.set_name(deformation_name) + self.ph_works.append(ph_work) + self.flow.register_work(ph_work) + + # Add task for electron DOS calculation to edos_work + if self.edos_ngkpt is not None: + edos_input = scf_input.make_edos_input(self.edos_ngkpt) + self.edos_work.register_nscf_task(edos_input, deps={ph_work[0]: "DEN"}).set_name(deformationname) + + if self.edos_ngkpt is not None: + self.flow.register_work(self.edos_work) + + self.flow.allocate(build=True) + + return super().on_all_ok() diff --git a/abipy/flowtk/zsisa.py b/abipy/flowtk/zsisa.py index 03b7dec1e..4029c81f9 100644 --- a/abipy/flowtk/zsisa.py +++ b/abipy/flowtk/zsisa.py @@ -152,6 +152,11 @@ def on_ok(self, sender): task = self.register_relax_task(relax_template.new_with_structure(structure, optcell=0)) self.relax_tasks_deformed.append(task) + # Build work for elastic properties (clamped-ions) + # activate internal strain and piezoelectric part. + #from abipy.flowtk.dfpt import ElasticWork + #elastic_work = ElasticWork.from_scf_input(scf_input, with_relaxed_ion=True, with_piezo=True) + self.flow.allocate(build=True) return super().on_ok(sender) @@ -250,5 +255,3 @@ def _on_ok(self): # self.restart() return results - -