diff --git a/abipy/examples/flows/run_qha_vzsisa.py b/abipy/examples/flows/run_qha_vzsisa.py index b0442e209..82190690e 100755 --- a/abipy/examples/flows/run_qha_vzsisa.py +++ b/abipy/examples/flows/run_qha_vzsisa.py @@ -22,9 +22,11 @@ def build_flow(options): __file__ = os.path.join(os.getcwd(), "run_qha_vzsisa.py") options.workdir = os.path.basename(__file__).replace(".py", "").replace("run_", "flow_") + from abipy.flowtk.psrepos import get_oncvpsp_pseudos + pseudos = get_oncvpsp_pseudos(xc_name="PBE", version="0.4") + # Initialize structure and pseudos structure = abilab.Structure.from_file(abidata.cif_file("si.cif")) - pseudos = abidata.pseudos("14si.pspnc") # Select k-mesh for electrons and q-mesh for phonons. #ngkpt = [2, 2, 2]; ngqpt = [2, 2, 2] @@ -41,7 +43,22 @@ def build_flow(options): ph_scales = [1, 1.02, 1.04] # EinfVib2(D) scf_input = abilab.AbinitInput(structure, pseudos) - scf_input.set_vars(ecut=8, nband=4, tolvrs=1e-8, nstep=50) + #scf_input.set_vars(ecut=8, nband=4, tolvrs=1e-8, nstep=50) + + # Set other important variables (consistent with tutorial) + # Aall the other DFPT runs will inherit these parameters. + scf_input.set_vars( + nband=4, + nline=10, + nbdbuf=0, + nstep=100, + ecut=8.0, + ecutsm=1.0, + #chksymbreak=1, + occopt=1, + tolvrs=1.0e-18, # SCF stopping criterion (modify default) + ) + scf_input.set_kmesh(ngkpt=ngkpt, shiftk=[0, 0, 0]) return VzsisaFlow.from_scf_input(options.workdir, scf_input, bo_scales, ph_scales, ngqpt, diff --git a/abipy/flowtk/zsisa.py b/abipy/flowtk/zsisa.py new file mode 100644 index 000000000..613365a48 --- /dev/null +++ b/abipy/flowtk/zsisa.py @@ -0,0 +1,261 @@ +# 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.flowtk.works import Work, PhononWork +from abipy.flowtk.tasks import RelaxTask +from abipy.flowtk.flows import Flow + + +class ZsisaFlow(Flow): + """ + Flow for QHA calculations with the VZSISA approach. + + .. rubric:: Inheritance Diagram + .. inheritance-diagram:: ZsisaFlow + """ + + @classmethod + def from_scf_input(cls, workdir, scf_input, eps, ngqpt, with_becs, with_quad, + edos_ngkpt=None, manager=None) -> ZsisaFlow: + """ + Build a flow for QHA calculations from an |AbinitInput| for GS-SCF calculation. + + Args: + workdir: Working directory of the flow. + eps: + scf_input: |AbinitInput| for GS-SCF run used as template to generate the other inputs. + strains_a: List of scaling factors for the a and b lattice vectors. + strains_c: List of scaling factors for the c lattice vectors. + 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) + + # optcell = 2: full optimization of cell geometry ( + work = ZsisaWork.from_scf_input(scf_input, eps, ngqpt, with_becs, with_quad, + ionmov=2, edos_ngkpt=edos_ngkpt) + + flow.register_work(work) + 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 = {} + + ## 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_vol] + + #entries, gsr_relax_volumes = [], [] + #for task in work.relax_tasks_vol: + # with task.open_gsr() as gsr: + # 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_entries"] = entries + #data["gsr_relax_volumes"] = gsr_relax_volumes + + #data["ddb_paths"] = [ph_work.outdir.has_abiext("DDB") for ph_work in work.ph_works] + #data["ddb_relax_volumes"] = [ph_work[0].input.structure.volume for ph_work in work.ph_works] + + #data["gsr_edos_path"] = [] if not work.edos_work else [task.gsr_path for task in work.edos_work] + + #mjson_write(data, self.outdir.path_in("zsisa.json"), indent=4) + + return super().finalize() + + +class ZsisaWork(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:: ZsisaWork + """ + + @classmethod + def from_scf_input(cls, scf_input, eps, ngqpt, + with_becs: bool, with_quad: bool, + ionmov: int, edos_ngkpt=None) -> ZsisaWork: + """ + Build the work from an |AbinitInput| representing a GS-SCF calculation. + + Args: + scf_input: |AbinitInput| for GS-SCF used as template to generate the other inputs. + with_becs: Activate calculation of Electric field and Born effective charges. + ionmov: Abinit input variables. + edos_ngkpt: Three integers defining the the k-sampling for the computation of the + electron DOS with the relaxed structures. Useful for metals + in which the electronic contribution should be included. + None disables the computation of the e-DOS. + """ + work = cls() + + # Save attributes in work + work.initial_scf_input = scf_input + work.eps = float(eps) + 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 = 2: full optimization of cell geometry ( + relax_template.set_vars(optcell=2, ionmov=ionmov, tolvrs=1e-8, toldff=1.e-6) + #if optcell is not None and optcell != 0: + relax_template.set_vars_ifnotin(ecutsm=0.5, 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() + v0 = relaxed_structure.volume + self.deformed_structures_dict = generate_deformations(relaxed_structure, eps=self.eps) + + relax_template = self.relax_template + self.relax_tasks_vol = [] + for structure in self.deformed_structures_dict.values(): + # Relax at fixed unit cell. + new_input = relax_template.new_with_structure(structure, optcell=0) + task = self.register_relax_task(new_input) + self.relax_tasks_vol.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 = None if self.edos_ngkpt is None else Work() + + # Build phonon works for the different relaxed structures. + self.ph_works = [] + + start = 1 + for task, def_name in zip(self[start:], self.deformed_structures_dict.keys(), strict=True): + relaxed_structure = task.get_final_structure() + scf_input = self.initial_scf_input.new_with_structure(relaxed_structure) + ph_work = PhononWork.from_scf_input(scf_input, self.ngqpt, is_ngqpt=True, tolerance=None, + with_becs=self.with_becs, ddk_tolerance=None) + + ph_work.set_name(def_name) + self.ph_works.append(ph_work) + self.flow.register_work(ph_work) + + # Add electron DOS calculation. + 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"}) + + if self.edos_ngkpt is not None: self.flow.register_work(self.edos_work) + self.flow.allocate(build=True) + + # Build phonon works for the different relaxed structures. + #for task, bo_scale in zip(self.relax_tasks_vol, self.strains_a): + # if all(abs(bo_scale - self.strains_c)) > 1e-3: continue + # relaxed_structure = task.get_final_structure() + # scf_input = self.initial_scf_input.new_with_structure(relaxed_structure) + + # ph_work = PhononWork.from_scf_input(scf_input, self.ngqpt, is_ngqpt=True, tolerance=None, + # with_becs=self.with_becs, with_quad=self.with_quad, + # ddk_tolerance=None) + + # # Reduce the number of files produced in the DFPT tasks to avoid possible disk quota issues. + # prtvars = dict(prtwf=-1, prtden=0, prtpot=0) + # for task in ph_work[1:]: + # task.input.set_vars(**prtvars) + + # self.flow.register_work(ph_work) + # self.ph_works.append(ph_work) + + # # Add electron DOS calculation. + # if self.edos_work 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"}) + + #if self.edos_ngkpt is not None: self.flow.register_work(self.edos_work) + self.flow.allocate(build=True) + + return super().on_all_ok() + + +#class ThermalRelaxTask(RelaxTask): +# +# def _on_ok(self): +# results = super()._on_ok() +# # Check for convergence. +# #if not self.collinear_done: +# # self.input.set_vars(strtarget=strtarget) +# # self.finalized = False +# # self.restart() +# +# return results +# +# +#class ThermalRelaxWork(Work): +# """ +# .. rubric:: Inheritance Diagram +# .. inheritance-diagram:: ThermalRelaxWork +# """ +# +# @classmethod +# def from_relax_input(cls, relax_input, qha, temperatures, pressures): +# """ +# """ +# work = cls() +# +# work.temperatures = temperatures +# work.pressures = pressures +# #work.qha = qha +# +# for pressure in pressures: +# for temperature in temperatures: +# strtarget = qha.get_strtarget(temperature, pressure) +# new_input = relax_input.new_with_vars(strtarget=strtarget) +# work.register_relax_task(new_input) +# +# return work +# +# def on_all_ok(self): +# """ +# Implement the post-processing step at the end of the Work. +# """ +# return super().on_all_ok()