From 6b39602fcf7ca54956e19d99e98db511e51303f3 Mon Sep 17 00:00:00 2001 From: Guido Petretto Date: Thu, 9 Jun 2022 17:32:15 +0200 Subject: [PATCH 1/4] fix nband setting in factories --- abipy/abio/factories.py | 46 ++++++++++++++++++++++++++---- abipy/abio/tests/test_factories.py | 33 +++++++++++++++++++++ 2 files changed, 74 insertions(+), 5 deletions(-) diff --git a/abipy/abio/factories.py b/abipy/abio/factories.py index 8e37cdefb..b06b2e277 100644 --- a/abipy/abio/factories.py +++ b/abipy/abio/factories.py @@ -154,6 +154,34 @@ def _find_scf_nband(structure, pseudos, electrons, spinat=None): return int(nband) +def _find_nscf_nband_from_gsinput(gs_input: AbinitInput) -> int: + """ + Find the value of nband for a NSCF calculation based on the input of a previous gs calculation. + + Args: + gs_input (AbinitInput): the ground state input. + + Returns: + the value of nband. + """ + scf_nband = gs_input.get("nband") + if scf_nband is None: + occopt = gs_input.get("occopt", 1) + tsmear = gs_input.get("tsmear", 0.01) + smearing = aobj.Smearing(occopt, tsmear) + + # only nsppol is used in _find_scf_nband, so just set that with a meaningful value + nsppol = gs_input.get("nsppol", 1) + spin_mode = aobj.SpinMode(mode="test", nsppol=nsppol, nspinor=1, nspden=nsppol) + + charge = gs_input.get("charge", 0) + electrons = aobj.Electrons(spin_mode=spin_mode, smearing=smearing, charge=charge) + spinat = gs_input.get('spinat', None) + scf_nband = _find_scf_nband(gs_input.structure, gs_input.pseudos, electrons, spinat) + + return scf_nband + 10 + + def _get_shifts(shift_mode: str, structure: Structure): """ Gives the shifts based on the selected shift mode and on the symmetry of the structure. @@ -1153,7 +1181,7 @@ def ebands_from_gsinput(gs_input, nband=None, ndivsm=15, accuracy="normal") -> A nscf_ksampling = aobj.KSampling.path_from_structure(ndivsm, gs_input.structure) if nband is None: - nband = gs_input.get("nband", gs_input.structure.num_valence_electrons(gs_input.pseudos)) + 10 + nband = _find_nscf_nband_from_gsinput(gs_input) bands_input.set_vars(nscf_ksampling.to_abivars()) bands_input.set_vars(nband=nband, iscf=-2) @@ -1170,7 +1198,11 @@ def dos_from_gsinput(gs_input, dos_kppa, nband=None, accuracy="normal", pdos=Fal dos_ksampling = aobj.KSampling.automatic_density(dos_input.structure, dos_kppa, chksymbreak=0) dos_input.set_vars(dos_ksampling.to_abivars()) - dos_input.set_vars(iscf=-2, ionmov=0) + + if nband is None: + nband = _find_nscf_nband_from_gsinput(gs_input) + + dos_input.set_vars(nband=nband, iscf=-2) dos_input.set_vars(_stopping_criterion("nscf", accuracy)) if pdos: @@ -1235,13 +1267,17 @@ def scf_for_phonons(structure, pseudos, kppa=None, ecut=None, pawecutdg=None, nb spin_mode="polarized", smearing="fermi_dirac:0.1 eV", charge=0.0, scf_algorithm=None, shift_mode="Symmetric"): + # add the band for nbdbuf, if needed + nbdbuf = 4 + if nband is not None: + nband += nbdbuf + abiinput = scf_input(structure=structure, pseudos=pseudos, kppa=kppa, ecut=ecut, pawecutdg=pawecutdg, nband=nband, accuracy=accuracy, spin_mode=spin_mode, smearing=smearing, charge=charge, scf_algorithm=scf_algorithm, shift_mode=shift_mode) - nbdbuf = 4 - # with no smearing set the minimum number of bands plus some nbdbuf - if smearing is None: + # with no bands set and no smearing the minimum number of bands plus some nbdbuf + if nband is None and smearing is None: nval = structure.num_valence_electrons(pseudos) nval -= abiinput['charge'] nband = int(round(nval / 2) + nbdbuf) diff --git a/abipy/abio/tests/test_factories.py b/abipy/abio/tests/test_factories.py index 113a15ae1..d2ed3299f 100644 --- a/abipy/abio/tests/test_factories.py +++ b/abipy/abio/tests/test_factories.py @@ -7,6 +7,7 @@ from abipy.abio.factories import * from abipy.abio.factories import (BandsFromGsFactory, IoncellRelaxFromGsFactory, HybridOneShotFromGsFactory, ScfForPhononsFactory, PhononsFromGsFactory, PiezoElasticFactory, PiezoElasticFromGsFactory, ShiftMode) +from abipy.abio.factories import _find_nscf_nband_from_gsinput import json write_inputs_to_json = False @@ -23,6 +24,38 @@ def test_shiftmode(self): ShiftMode.from_object({}) +class HelperTest(AbipyTest): + + @classmethod + def setUpClass(cls): + cls.si_structure = abidata.structure_from_cif("si.cif") + cls.si_pseudo = abidata.pseudos("14si.pspnc") + + def test_find_nscf_nband_from_gsinput(self): + gsi = AbinitInput(self.si_structure, self.si_pseudo) + gsi.set_vars(nband=100) + assert _find_nscf_nband_from_gsinput(gsi) == 110 + + gsi = AbinitInput(self.si_structure, self.si_pseudo) + gsi.set_vars(occopt=1, nsppol=1) + assert _find_nscf_nband_from_gsinput(gsi) == 18 + + sc4 = self.si_structure.copy() + sc4.make_supercell([4, 4, 4]) + gsi = AbinitInput(sc4, self.si_pseudo) + gsi.set_vars(occopt=1, nsppol=1) + assert _find_nscf_nband_from_gsinput(gsi) == 292 + + gsi.set_vars(occopt=3) + assert _find_nscf_nband_from_gsinput(gsi) == 318 + + gsi.set_vars(nsppol=2) + assert _find_nscf_nband_from_gsinput(gsi) == 318 + + gsi.set_vars(spinat=[[0, 0, 1]] * len(sc4)) + assert _find_nscf_nband_from_gsinput(gsi) == 382 + + class FactoryTest(AbipyTest): def setUp(self): From 02a6f97eaf237d3951e0a32b7a042fcd1c741226 Mon Sep 17 00:00:00 2001 From: Guido Petretto Date: Sun, 19 Jun 2022 15:58:06 +0200 Subject: [PATCH 2/4] refactor nscf input factories --- abipy/abio/factories.py | 100 +++++++++++++++++++++++---- abipy/abio/tests/test_factories.py | 24 ++++++- abipy/electrons/ebands.py | 11 +++ abipy/electrons/tests/test_ebands.py | 2 + 4 files changed, 119 insertions(+), 18 deletions(-) diff --git a/abipy/abio/factories.py b/abipy/abio/factories.py index b06b2e277..975f73f07 100644 --- a/abipy/abio/factories.py +++ b/abipy/abio/factories.py @@ -1162,15 +1162,18 @@ def scf_input(structure, pseudos, kppa=None, ecut=None, pawecutdg=None, nband=No return abinit_input -def ebands_from_gsinput(gs_input, nband=None, ndivsm=15, accuracy="normal") -> AbinitInput: +def ebands_from_gsinput(gs_input, nband=None, ndivsm=15, accuracy="normal", + projection=None) -> AbinitInput: """ Return an |AbinitInput| object to compute a band structure from a GS SCF input. Args: - gs_input: - nband: - ndivsm: - accuracy: + gs_input: the |AbinitInput| that was used to calculated the charge density. + nband: the number of bands to be used for the calculation. If None it will be + automatically generated. + ndivsm: Number of divisions used to sample the smallest segment of the k-path. + accuracy: Accuracy of the calculation. + projection: which projection should be performed. If None no projection, otherwise "l" or "lm" Return: |AbinitInput| """ @@ -1187,27 +1190,94 @@ def ebands_from_gsinput(gs_input, nband=None, ndivsm=15, accuracy="normal") -> A bands_input.set_vars(nband=nband, iscf=-2) bands_input.set_vars(_stopping_criterion("nscf", accuracy)) + if projection is not None: + if "l" in projection: + bands_input.set_vars(prtdos=3) + if "m" in projection: + bands_input.set_vars(prtdosm=1) + return bands_input -def dos_from_gsinput(gs_input, dos_kppa, nband=None, accuracy="normal", pdos=False): +def nscf_from_gsinput(gs_input, kppa=None, nband=None, accuracy="normal", + shift_mode="Monkhorst-Pack") -> AbinitInput: + """ + Return an |AbinitInput| object to perform a NSCF calculation from a GS SCF input. + + Args: + gs_input: the |AbinitInput| that was used to calculated the charge density. + kppa: defines the kpt sampling used for the NSCF run. If None the kpoint sampling and + shifts will be the same as in the SCF input. + nband: the number of bands to be used for the calculation. If None it will be + automatically generated. + accuracy: accuracy of the calculation. + shift_mode: the mode to be used for the shifts. Options are "Gamma", "Monkhorst-Pack", + "Symmetric", "OneSymmetric". See ShiftMode object for more details. Only used if kppa + is not None. + Return: |AbinitInput| + """ # create a copy to avoid messing with the previous input - dos_input = gs_input.deepcopy() - dos_input.pop_irdvars() + nscf_input = gs_input.deepcopy() + nscf_input.pop_irdvars() - dos_ksampling = aobj.KSampling.automatic_density(dos_input.structure, dos_kppa, chksymbreak=0) - dos_input.set_vars(dos_ksampling.to_abivars()) + if kppa is not None: + shift_mode = ShiftMode.from_object(shift_mode) + shifts = _get_shifts(shift_mode, gs_input.structure) + dos_ksampling = aobj.KSampling.automatic_density(nscf_input.structure, kppa, chksymbreak=0, shifts=shifts) + nscf_input.set_vars(dos_ksampling.to_abivars()) if nband is None: nband = _find_nscf_nband_from_gsinput(gs_input) - dos_input.set_vars(nband=nband, iscf=-2) - dos_input.set_vars(_stopping_criterion("nscf", accuracy)) + nscf_input.set_vars(nband=nband, iscf=-2) + nscf_input.set_vars(_stopping_criterion("nscf", accuracy)) + + return nscf_input + + +def dos_from_gsinput(gs_input, kppa=None, nband=None, accuracy="normal", dos_method="tetra", + projection="l", shift_mode="Monkhorst-Pack") -> AbinitInput: + """ + Return an |AbinitInput| object to perform a DOS calculation from a GS SCF input. + + Args: + gs_input: the |AbinitInput| that was used to calculated the charge density. + kppa: defines the kpt sampling used for the NSCF run. If None the kpoint sampling and + shifts will be the same as in the SCF input. + nband: the number of bands to be used for the calculation. If None it will be + automatically generated. + accuracy: accuracy of the calculation. + dos_method: method to calculate the DOS in abinit (NB: not the one used from postprocessing + in abipy). Set to "tetra" for the tetrahedron method (prtdos 2 or 3). If "smearing", + occopt and tsmear will be taken from gs_input else a "smearing-type: smearing value" + (prtdos 1 or 4). + projection: which projection should be performed. If None no projection, otherwise "l" or "lm" + shift_mode: the mode to be used for the shifts. Options are "Gamma", "Monkhorst-Pack", + "Symmetric", "OneSymmetric". See ShiftMode object for more details. Only used if kppa + is not None. + + Return: |AbinitInput| + """ + dos_input = nscf_from_gsinput(gs_input, kppa=kppa, nband=nband, accuracy=accuracy, shift_mode=shift_mode) + + if dos_method == "tetra": + if projection is not None and "l" in projection: + dos_input.set_vars(prtdos=3) + else: + dos_input.set_vars(prtdos=2) + else: + if dos_method != "smearing": + smear_obj = aobj.Smearing.as_smearing(dos_method) + dos_input.set_vars(smear_obj.to_abivars()) + + if projection is not None and "l" in projection.lower(): + dos_input.set_vars(prtdos=4) + else: + dos_input.set_vars(prtdos=1) - if pdos: - # FIXME - raise NotImplementedError() + if projection is not None and "m" in projection.lower(): + dos_input.set_vars(prtdosm=1) return dos_input diff --git a/abipy/abio/tests/test_factories.py b/abipy/abio/tests/test_factories.py index d2ed3299f..98d05c7c7 100644 --- a/abipy/abio/tests/test_factories.py +++ b/abipy/abio/tests/test_factories.py @@ -372,17 +372,35 @@ def test_scf_input(self): inp["ecut"] = 2 self.abivalidate_input(inp) - def test_ebands_dos_from_gsinput(self): + def test_nscf_ebands_dos_from_gsinput(self): """Testing ebands_from_gsinput and dos_from_gsinput""" - from abipy.abio.factories import ebands_from_gsinput, dos_from_gsinput + from abipy.abio.factories import ebands_from_gsinput, dos_from_gsinput, nscf_from_gsinput gs_inp = gs_input(self.si_structure, self.si_pseudo, kppa=None, ecut=2, spin_mode="unpolarized") + + nscf_inp = nscf_from_gsinput(gs_inp, kppa=None, nband=120) + self.assertArrayEqual(gs_inp["ngkpt"], nscf_inp["ngkpt"]) + self.assertEqual(nscf_inp["nband"], 120) + ebands_inp = ebands_from_gsinput(gs_inp, nband=None, ndivsm=15, accuracy="normal") self.abivalidate_input(ebands_inp) + ebands_inp = ebands_from_gsinput(gs_inp, nband=None, ndivsm=15, accuracy="normal", projection="lm") + self.assertEqual(ebands_inp["prtdos"], 3) + self.assertEqual(ebands_inp["prtdosm"], 1) + dos_kppa = 3000 - edos_inp = dos_from_gsinput(gs_inp, dos_kppa, nband=None, accuracy="normal", pdos=False) + edos_inp = dos_from_gsinput(gs_inp, dos_kppa, nband=None, accuracy="normal") self.abivalidate_input(edos_inp) + edos_inp = dos_from_gsinput(gs_inp, dos_method="smearing", projection="lm") + self.assertEqual(gs_inp["occopt"], edos_inp["occopt"]) + self.assertEqual(edos_inp["prtdos"], 4) + self.assertEqual(edos_inp["prtdosm"], 1) + + edos_inp = dos_from_gsinput(gs_inp, dos_method="marzari5: 0.01 eV", projection="l") + self.assertEqual(edos_inp["occopt"], 5) + self.assertNotIn("prtdosm", edos_inp) + factory_obj = BandsFromGsFactory(nband=None, ndivsm=15, accuracy="normal") self.assertMSONable(factory_obj) ebands_input_obj = factory_obj.build_input(gs_inp) diff --git a/abipy/electrons/ebands.py b/abipy/electrons/ebands.py index 26a1d6ccd..60c0fbf68 100644 --- a/abipy/electrons/ebands.py +++ b/abipy/electrons/ebands.py @@ -4741,6 +4741,17 @@ def plotly_up_minus_down(self, e0="fermie", fig=None, xlims=None, **kwargs): return fig + def to_pymatgen(self): + """ + Return a pymatgen DOS object from an Abipy |ElectronDos| object. + """ + + from pymatgen.electronic_structure.dos import Dos + den = {s: d.values for d, s in zip(self.spin_dos, [PmgSpin.up, PmgSpin.down])} + pmg_dos = Dos(energies=self.spin_dos[0].mesh, densities=den, efermi=self.fermie) + + return pmg_dos + class ElectronDosPlotter(NotebookWriter): """ diff --git a/abipy/electrons/tests/test_ebands.py b/abipy/electrons/tests/test_ebands.py index 9ed07a9ed..ddc9212ed 100644 --- a/abipy/electrons/tests/test_ebands.py +++ b/abipy/electrons/tests/test_ebands.py @@ -293,6 +293,8 @@ def test_silicon_ebands(self): with self.assertRaises(TypeError): ElectronDos.as_edos({}, {}) + assert si_ebands_kmesh.to_pymatgen() + mu = si_edos.find_mu(8) imu = si_edos.tot_idos.find_mesh_index(mu) self.assert_almost_equal(si_edos.tot_idos[imu][1], 8, decimal=2) From 8da7ca011217b0439cddd152a6a1ddebef70a68e Mon Sep 17 00:00:00 2001 From: Guido Petretto Date: Mon, 20 Jun 2022 13:23:48 +0200 Subject: [PATCH 3/4] improve handling of dos projection --- abipy/abio/factories.py | 29 ++++++++++++++++++++--------- abipy/abio/tests/test_factories.py | 6 ++++-- 2 files changed, 24 insertions(+), 11 deletions(-) diff --git a/abipy/abio/factories.py b/abipy/abio/factories.py index 975f73f07..99e31cc1c 100644 --- a/abipy/abio/factories.py +++ b/abipy/abio/factories.py @@ -1190,11 +1190,14 @@ def ebands_from_gsinput(gs_input, nband=None, ndivsm=15, accuracy="normal", bands_input.set_vars(nband=nband, iscf=-2) bands_input.set_vars(_stopping_criterion("nscf", accuracy)) - if projection is not None: - if "l" in projection: - bands_input.set_vars(prtdos=3) - if "m" in projection: - bands_input.set_vars(prtdosm=1) + if projection is None: + pass + elif projection == "l": + bands_input.set_vars(prtdos=3) + elif projection == "lm": + bands_input.set_vars(prtdos=3, prtdosm=1) + else: + raise ValueError(f"Unrecognized value for projection: {projection}") return bands_input @@ -1262,19 +1265,27 @@ def dos_from_gsinput(gs_input, kppa=None, nband=None, accuracy="normal", dos_met dos_input = nscf_from_gsinput(gs_input, kppa=kppa, nband=nband, accuracy=accuracy, shift_mode=shift_mode) if dos_method == "tetra": - if projection is not None and "l" in projection: + if projection is None: + dos_input.set_vars(prtdos=2) + elif projection == "l": dos_input.set_vars(prtdos=3) + elif projection == "lm": + dos_input.set_vars(prtdos=3, prtdosm=1) else: - dos_input.set_vars(prtdos=2) + ValueError(f"Unrecognized value for projection: {projection}") else: if dos_method != "smearing": smear_obj = aobj.Smearing.as_smearing(dos_method) dos_input.set_vars(smear_obj.to_abivars()) - if projection is not None and "l" in projection.lower(): + if projection is None: + dos_input.set_vars(prtdos=1) + elif projection == "l": dos_input.set_vars(prtdos=4) + elif projection == "lm": + raise ValueError("lm projection is only allowed for dos_method 'tetra'") else: - dos_input.set_vars(prtdos=1) + ValueError(f"Unrecognized value for projection: {projection}") if projection is not None and "m" in projection.lower(): dos_input.set_vars(prtdosm=1) diff --git a/abipy/abio/tests/test_factories.py b/abipy/abio/tests/test_factories.py index 98d05c7c7..560ec49ce 100644 --- a/abipy/abio/tests/test_factories.py +++ b/abipy/abio/tests/test_factories.py @@ -392,10 +392,12 @@ def test_nscf_ebands_dos_from_gsinput(self): edos_inp = dos_from_gsinput(gs_inp, dos_kppa, nband=None, accuracy="normal") self.abivalidate_input(edos_inp) - edos_inp = dos_from_gsinput(gs_inp, dos_method="smearing", projection="lm") + edos_inp = dos_from_gsinput(gs_inp, dos_method="smearing", projection="l") self.assertEqual(gs_inp["occopt"], edos_inp["occopt"]) self.assertEqual(edos_inp["prtdos"], 4) - self.assertEqual(edos_inp["prtdosm"], 1) + + with self.assertRaises(ValueError): + edos_inp = dos_from_gsinput(gs_inp, dos_method="smearing", projection="lm") edos_inp = dos_from_gsinput(gs_inp, dos_method="marzari5: 0.01 eV", projection="l") self.assertEqual(edos_inp["occopt"], 5) From cc342460a5096cfdeca5ad9f276bd2b3f41a8599 Mon Sep 17 00:00:00 2001 From: Guido Petretto Date: Mon, 11 Sep 2023 11:26:08 +0200 Subject: [PATCH 4/4] fix qha to fit new phonopy API --- abipy/dfpt/qha.py | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/abipy/dfpt/qha.py b/abipy/dfpt/qha.py index 892a80448..d15a62702 100644 --- a/abipy/dfpt/qha.py +++ b/abipy/dfpt/qha.py @@ -403,7 +403,8 @@ def write_phonopy_qha_inputs(self, tstart=0, tstop=2100, num=211, path=None) -> with open(os.path.join(path, "thermal_properties-{}.yaml".format(j)), 'wt') as f: f.write("\n".join(lines)) - def get_phonopy_qha(self, tstart=0, tstop=2100, num=211, eos='vinet', t_max=None, energy_plot_factor=None): + def get_phonopy_qha(self, tstart=0, tstop=2100, num=211, eos='vinet', t_max=None, + energy_plot_factor=None, pressure=None): """ Creates an instance of phonopy.qha.core.QHA that can be used generate further plots and output data. The object is returned right after the construction. The "run()" method should be executed @@ -419,6 +420,7 @@ def get_phonopy_qha(self, tstart=0, tstop=2100, num=211, eos='vinet', t_max=None "murnaghan" and "birch_murnaghan". Passed to phonopy's QHA. t_max: maximum temperature. Passed to phonopy's QHA. energy_plot_factor: factor multiplying the energies. Passed to phonopy's QHA. + pressure: pressure value, passed to phonopy. Returns: An instance of phonopy.qha.core.QHA """ @@ -437,7 +439,18 @@ def get_phonopy_qha(self, tstart=0, tstop=2100, num=211, eos='vinet', t_max=None en = self.energies + self.volumes * self.pressure / abu.eVA3_GPa - qha_p = QHA_phonopy(self.volumes, en, temperatures, cv, entropy, fe, eos, t_max, energy_plot_factor) + qha_p = QHA_phonopy( + volumes=self.volumes, + electronic_energies=en, + temperatures=temperatures, + cv=cv, + entropy=entropy, + fe_phonon=fe, + pressure=pressure, + eos=eos, + t_max=t_max, + energy_plot_factor=energy_plot_factor + ) return qha_p