diff --git a/abipy/core/structure.py b/abipy/core/structure.py index 63e449ea5..2196cd0ee 100644 --- a/abipy/core/structure.py +++ b/abipy/core/structure.py @@ -172,6 +172,7 @@ def display_structure(obj, **kwargs): def get_structures_from_file(filepath: PathLike, index) -> list[Structure]: """ + Read and return list of structures from filepath """ #if index is None: # index = -1 @@ -2409,16 +2410,31 @@ def calc_ngkpt(self, nksmall) -> np.ndarray: return ngkpt - #def any_to_ngkpt(self, any_var) -> np.ndarray: - # any_var = np.array(any_var) - # if len(any_var) == 3: - # return any_var - # if (nksmall := float(any_var)) > 0: - # return self.calc_ngkpt(nksmall) + def as_ngkpt(self, ngkpt) -> np.ndarray: + """ + Flexible API to compute the ABINIT variable ``ngkpt`` using different approaches. + + The following cases are supported: + - If `ngkpt` is a 1D vector with 3 items, return `ngkpt` as-is. + - If `ngkpt` is a positive float, interpret it as `nksmall` (a scaling factor for the k-point grid). + - If `ngkpt` is a negative float, interpret it as the desired number of k-points per atom. + + This method allows users to flexibly specify the k-point grid based on their preferred input format. + """ + ngkpt = np.array(ngkpt) + + if ngkpt.ndim == 1 and len(ngkpt) == 3: + return ngkpt + + if (nksmall := float(ngkpt)) > 0: + return self.calc_ngkpt(nksmall) + + if (kppa := -float(ngkpt)) > 0: + import pymatgen.io.abinit.abiobjects as aobj + ksampling = aobj.KSampling.automatic_density(self, kppa, chksymbreak=0, shifts=(0,0,0)) + return ksampling.to_abivars()["ngkpt"] - # import pymatgen.io.abinit.abiobjects as aobj - # ksampling = aobj.KSampling.automatic_density(self, kppa, chksymbreak=0, shifts=(0,0,0)) - # return ksampling.to_abivars()["ngkpt"] + raise ValueError(f"Don't know how to convert {type(ngkpt)=}, {ngkpt=} to k-mesh!") def calc_shiftk(self, symprec=0.01, angle_tolerance=5) -> np.ndarray: """ diff --git a/abipy/core/tests/test_structure.py b/abipy/core/tests/test_structure.py index 82fa04f26..aa06c2f0f 100644 --- a/abipy/core/tests/test_structure.py +++ b/abipy/core/tests/test_structure.py @@ -159,6 +159,7 @@ def test_utils(self): assert len(si.abi_string) assert si.reciprocal_lattice == si.lattice.reciprocal_lattice + # Test k-sampling methods. kptbounds = si.calc_kptbounds() ksamp = si.calc_ksampling(nksmall=10) @@ -168,6 +169,13 @@ def test_utils(self): self.assert_equal(ksamp.ngkpt, [10, 10, 10]) self.assert_equal(ksamp.shiftk, shiftk) + # Test as_ngkpt API + self.assert_equal(si.as_ngkpt([1, 2, 3]), [1, 2, 3]) + self.assert_equal(si.as_ngkpt(3), [3, 3, 3]) + self.assert_equal(si.as_ngkpt(-1000), [8, 8, 8]) + with self.assertRaises(ValueError): + si.as_ngkpt("foo") + lif = Structure.from_abistring(""" acell 7.7030079150 7.7030079150 7.7030079150 Angstrom rprim 0.0000000000 0.5000000000 0.5000000000 diff --git a/abipy/examples/flows/run_conducwork.py b/abipy/examples/flows/run_conducwork.py index de43a2b84..604e5aa26 100755 --- a/abipy/examples/flows/run_conducwork.py +++ b/abipy/examples/flows/run_conducwork.py @@ -15,8 +15,7 @@ from abipy.abio.factories import conduc_kerange_from_inputs -def make_scf_input(structure, pseudos, ngkpt=(2,2,2), shiftk=(0,0,0), - **variables): +def make_scf_input(structure, pseudos, ngkpt=(2,2,2), shiftk=(0,0,0), **variables): """Build and return SCF input given the structure and pseudopotentials""" scf_inp = abilab.AbinitInput(structure, pseudos=pseudos) @@ -31,8 +30,7 @@ def make_scf_input(structure, pseudos, ngkpt=(2,2,2), shiftk=(0,0,0), return scf_inp -def make_nscf_input(structure, pseudos, ngkpt=(2,2,2), shiftk=(0,0,0), - **variables): +def make_nscf_input(structure, pseudos, ngkpt=(2,2,2), shiftk=(0,0,0), **variables): """Build and return NSCF input given the structure and pseudopotentials""" scf_inp = abilab.AbinitInput(structure, pseudos=pseudos) @@ -76,14 +74,14 @@ def build_flow(options): # Kerange Variables nbr_proc = 4 - ngqpt_fine = [16, 16, 16] # The sigma_ngkpt grid must be divisible by the qpt grid + ngqpt_fine = [16, 16, 16] # The sigma_ngkpt grid must be divisible by the qpt grid sigma_ngkpt = [16, 16, 16] - einterp = [1, 5, 0, 0] # Star functions Interpolation + einterp = [1, 5, 0, 0] # Star functions Interpolation sigma_erange = [-0.3, -0.3, "eV"] # Negative value for metals flow = flowtk.Flow(workdir=options.workdir) - # Create inputs Object + # Create inputs scf_input = make_scf_input(structure, pseudos, tolvrs=1e-12, ngkpt=ngkpt, @@ -118,17 +116,18 @@ def build_flow(options): sigma_erange=sigma_erange, einterp=einterp, boxcutmin=boxcutmin, # 1.1 is the default value of the function - mixprec=mixprec # 1 is the default value of the function + mixprec=mixprec # 1 is the default value of the function ) # Here we can change multi to change the variable of a particular dataset - conduc_work = flowtk.ConducWork.from_phwork(phwork=ph_work, # Linking the DDB and DVDB via a |PhononWork| - multi=multi, # The multidataset object + conduc_work = flowtk.ConducWork.from_phwork(phwork=ph_work, # Linking the DDB and DVDB via a PhononWork + multi=multi, # The multidataset object nbr_proc=nbr_proc, # Needed to parallelize the calculation flow=flow, with_kerange=True, # Using Kerange - omp_nbr_thread=1) # 1 is the default value of the function + omp_nbr_thread=1) # 1 is the default value of the function + # If you already have the DDB and DVDB, use from_filepath(ddb_path, dvdb_path, multi, ...) instead of from_phwork flow.register_work(conduc_work) @@ -157,5 +156,3 @@ def main(options): if __name__ == "__main__": sys.exit(main()) - -############################################ diff --git a/abipy/examples/flows/run_qha_vzsisa.py b/abipy/examples/flows/run_qha_vzsisa.py index bd9ff5a73..0dc6b98b0 100755 --- a/abipy/examples/flows/run_qha_vzsisa.py +++ b/abipy/examples/flows/run_qha_vzsisa.py @@ -10,7 +10,7 @@ import abipy.data as abidata from abipy import flowtk -from abipy.flowtk.qha import vZSISAFlow +from abipy.flowtk.qha import VzsisaFlow def build_flow(options): @@ -26,24 +26,25 @@ def build_flow(options): structure = abilab.Structure.from_file(abidata.cif_file("si.cif")) pseudos = abidata.pseudos("14si.pspnc") - # Build input for GS calculation. + # Select k-mesh for electrons and q-mesh for phonons. #ngkpt = [2, 2, 2]; ngqpt = [2, 2, 2] #ngkpt = [4, 4, 4]; ngqpt = [4, 4, 4] ngkpt = [2, 2, 2]; ngqpt = [1, 1, 1] + with_becs = False with_quad = False #with_quad = not structure.has_zero_dynamical_quadrupoles #bo_scales = [0.96, 0.98, 1.0, 1.02, 1.04, 1.06] #ph_scales = [0.98, 1.0, 1.02, 1.04, 1.06] # EinfVib4(D) - bo_scales = [0.96, 0.98, 1, 1.02, 1.04] # EinfVib4(S) - ph_scales = [1, 1.02, 1.04] # EinfVib2(D) + bo_scales = [0.96, 0.98, 1, 1.02, 1.04] # EinfVib4(S) + 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) #, paral_kgb=0) + scf_input.set_vars(ecut=8, nband=4, tolvrs=1e-8, nstep=50) scf_input.set_kmesh(ngkpt=ngkpt, shiftk=[0, 0, 0]) - return vZSISAFlow.from_scf_input(options.workdir, scf_input, bo_scales, ph_scales, ngqpt, + return VzsisaFlow.from_scf_input(options.workdir, scf_input, bo_scales, ph_scales, ngqpt, with_becs, with_quad, edos_ngkpt=None) diff --git a/abipy/flowtk/qha.py b/abipy/flowtk/qha.py index 2fc3dd0bd..0b773e90c 100644 --- a/abipy/flowtk/qha.py +++ b/abipy/flowtk/qha.py @@ -13,25 +13,25 @@ from abipy.flowtk.flows import Flow -class vZSISAFlow(Flow): +class VzsisaFlow(Flow): """ - Flow for QHA calculations with the vZSISA approach + Flow for QHA calculations with the VZSISA approach. .. rubric:: Inheritance Diagram - .. inheritance-diagram:: vZSISAFlow + .. inheritance-diagram:: VzsisaFlow """ @classmethod def from_scf_input(cls, workdir, scf_input, bo_scales, ph_scales, ngqpt, with_becs, with_quad, - edos_ngkpt=None, manager=None) -> vZSISAFlow: + edos_ngkpt=None, manager=None) -> VzsisaFlow: """ 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_scales: - ph_scales: + bo_scales: List of scaling factors for the BO volumes. + ph_scales: List of scaling factors for the BO volumes. 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` @@ -45,8 +45,9 @@ def from_scf_input(cls, workdir, scf_input, bo_scales, ph_scales, ngqpt, with_be """ flow = cls(workdir=workdir, manager=manager) - work = vZSISAWork.from_scf_input(scf_input, bo_scales, ph_scales, ngqpt, with_becs, with_quad, - optcell=2, ionmov=2, edos_ngkpt=edos_ngkpt) + # optcell = 3: constant-volume optimization of cell geometry + work = VzsisaWork.from_scf_input(scf_input, bo_scales, ph_scales, ngqpt, with_becs, with_quad, + optcell=3, ionmov=2, edos_ngkpt=edos_ngkpt) flow.register_work(work) return flow @@ -61,47 +62,45 @@ def finalize(self): # 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] - data["gsr_relax_volumes"] = [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["gsr_edos_path"] = [] - if work.edos_work: - data["gsr_edos_paths"] = [task.gsr_path for task in work.edos_work] - - #entries = [] - #items = zip(self.relax_and_phonon_work.relax_tasks_vol, self.relax_and_phonon_work.ph_works_vol) - #for ivol, (relax_task, ph_work) in enumerate(items): - # ddb_path = ph_work.outdir.has_abiext("DDB") - # with relax_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_edos_path=None, - # ddb_path=ddb_path, - # ) - - #data["entries_for_each_volume"] = entries + 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("vzsisa.json"), indent=4) return super().finalize() -class vZSISAWork(Work): +class VzsisaWork(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:: vZSISAWork + .. inheritance-diagram:: VzsisaWork """ @classmethod def from_scf_input(cls, scf_input, bo_scales, ph_scales, ngqpt, with_becs: bool, with_quad: bool, - optcell: int, ionmov: int, edos_ngkpt=None) -> vZSISAWork: + optcell: int, ionmov: int, edos_ngkpt=None) -> VzsisaWork: """ Build the work from an |AbinitInput| representing a GS-SCF calculation. @@ -172,10 +171,16 @@ def on_all_ok(self): if all(abs(bo_scale - self.ph_scales)) > 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) diff --git a/abipy/flowtk/tests/test_qha.py b/abipy/flowtk/tests/test_qha.py index 609b2eeaf..689463113 100644 --- a/abipy/flowtk/tests/test_qha.py +++ b/abipy/flowtk/tests/test_qha.py @@ -1,19 +1,27 @@ """Tests for qha module.""" from abipy.core.testing import AbipyTest -from abipy.flowtk.qha import QhaFlow +from abipy.flowtk.qha import VzsisaFlow -class TestQha(AbipyTest): +class TestVzsisa(AbipyTest): - def test_qhaflow(self): - """Testing QhaFlow.""" + def test_vzsisa_flow(self): + """Testing VzsisaFlow.""" workdir = self.mkdtemp() - scf_input = self.get_gsinput_alas_ngkpt(ngkpt=[4, 4, 4]) - #v0 = scf_input.structure.volume - #volumes = [0.98 * v0, v0, v0 * 1.02] - flow = QhaFlow.from_scf_input(workdir, scf_input, ngqpt=[2, 2, 2], - with_becs=False, edos_ngkpt=[4, 4, 4]) #, metadata={"mpi_id": 123}) + ngkpt = [4, 4, 4] + ngqpt = [2, 2, 2] + scf_input = self.get_gsinput_alas_ngkpt(ngkpt=ngkpt) + + bo_scales = [0.96, 0.98, 1, 1.02, 1.04] # EinfVib4(S) + ph_scales = [1, 1.02, 1.04] # EinfVib2(D) + + with_becs = False + with_quad = False + #with_quad = not structure.has_zero_dynamical_quadrupoles + + flow = VzsisaFlow.from_scf_input(workdir, scf_input, bo_scales, ph_scales, ngqpt, + with_becs, with_quad, edos_ngkpt=None) flow.allocate() flow.check_status()