Skip to content

Commit

Permalink
Gw factory functions (#249)
Browse files Browse the repository at this point in the history 
* fix nband setting in factories

* refactor nscf input factories

* improve handling of dos projection

* Added ion_ioncell reference input json files.
Added comparison of the reference input json files.

* Added scr_* and self_* input factories.

* Updated g0w0_with_ppmodel_inputs comments.

Co-authored-by: Guido Petretto <[email protected]>
  • Loading branch information
@davidwaroquiers @gpetretto
davidwaroquiers and gpetretto authored Jun 27, 2022
1 parent 2ca0177 commit f8cae90
Show file tree
Hide file tree
Showing 4 changed files with 267 additions and 69 deletions.
264 changes: 198 additions & 66 deletions abipy/abio/factories.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -379,6 +407,42 @@ def ion_ioncell_relax_and_ebands_input(structure, pseudos,
return relax_multi + ebands_multi


def scr_from_nscfinput(nscf_input, nband=None, ecuteps=3.0, ecutwfn=None, inclvkb=2, w_type="RPA", sc_mode="one_shot", hilbert=None, accuracy="normal"):
"""Return a screening input."""
scr_input = nscf_input.deepcopy()
scr_input.pop_irdvars()
if nband is None:
nband = nscf_input.get("nband")
screening = aobj.Screening(ecuteps, nband, w_type=w_type, sc_mode=sc_mode,
hilbert=hilbert, ecutwfn=ecutwfn, inclvkb=inclvkb)

scr_input.set_vars(screening.to_abivars())
scr_input.set_vars(_stopping_criterion("screening", accuracy)) # Dummy

return scr_input


def sigma_from_inputs(nscf_input, scr_input, nband=None, ecutwfn=None, ecuteps=None, ecutsigx=None, ppmodel="godby", gw_qprange=1, accuracy="normal"):
"""Return a sigma input."""
self_input = nscf_input.deepcopy()
self_input.pop_irdvars()
if nband is None:
nband = nscf_input.get("nband")
screening = aobj.Screening(ecuteps=scr_input["ecuteps"], nband=scr_input["nband"],
w_type="RPA",
sc_mode="one_shot",
hilbert=None,
ecutwfn=scr_input["ecutwfn"],)
ecuteps = ecuteps if ecuteps is not None else screening.ecuteps
self_energy = aobj.SelfEnergy(se_type="gw", sc_mode="one_shot", nband=nband, ecutsigx=ecutsigx, screening=screening,
gw_qprange=gw_qprange, ppmodel=ppmodel, ecuteps=ecuteps, ecutwfn=ecutwfn, gwpara=2)

self_input.set_vars(self_energy.to_abivars())
self_input.set_vars(_stopping_criterion("sigma", accuracy)) # Dummy

return self_input


def g0w0_with_ppmodel_inputs(structure, pseudos,
kppa, nscf_nband, ecuteps, ecutsigx,
ecut=None, pawecutdg=None, shifts=(0.0, 0.0, 0.0),
Expand Down Expand Up @@ -418,55 +482,34 @@ def g0w0_with_ppmodel_inputs(structure, pseudos,
"""

structure = Structure.as_structure(structure)
multi = MultiDataset(structure, pseudos, ndtset=4)

# Set the cutoff energies.
multi.set_vars(_find_ecut_pawecutdg(ecut, pawecutdg, multi.pseudos, accuracy))

scf_ksampling = aobj.KSampling.automatic_density(structure, kppa, chksymbreak=0, shifts=shifts)
# Scf input
# Note that kppa and shift_mode here are dummy as, they will be overwritten just after.
scf_inp = scf_input(structure=structure, pseudos=pseudos, kppa=kppa,
ecut=ecut, pawecutdg=pawecutdg, nband=None, accuracy=accuracy,
spin_mode=spin_mode, smearing=smearing, charge=charge, scf_algorithm=scf_algorithm,
shift_mode="Monkhorst-Pack")
ksampling = aobj.KSampling.automatic_density(structure, kppa, chksymbreak=0, shifts=shifts)
scf_inp.set_vars(ksampling.to_abivars())
scf_inp.set_vars(istwfk="*1")
# The following is to keep the previous behavior
# - The spinat is not set
# - The number of bands is not adapted for spinat
# TODO: Should we consider changing that and update the reference files accordingly ?
scf_inp.pop_vars('spinat')
scf_electrons = aobj.Electrons(spin_mode=spin_mode, smearing=smearing, algorithm=scf_algorithm,
charge=charge, nband=None, fband=None)

if scf_electrons.nband is None:
scf_electrons.nband = _find_scf_nband(structure, multi.pseudos, scf_electrons)

multi[0].set_vars(scf_ksampling.to_abivars())
multi[0].set_vars(scf_electrons.to_abivars())
multi[0].set_vars(_stopping_criterion("scf", accuracy))

nscf_ksampling = aobj.KSampling.automatic_density(structure, kppa, chksymbreak=0, shifts=shifts)
nscf_electrons = aobj.Electrons(spin_mode=spin_mode, smearing=smearing, algorithm={"iscf": -2},
charge=charge, nband=nscf_nband, fband=None)

multi[1].set_vars(nscf_ksampling.to_abivars())
multi[1].set_vars(nscf_electrons.to_abivars())
multi[1].set_vars(_stopping_criterion("nscf", accuracy))
# nbdbuf

# Screening.
if scr_nband is None: scr_nband = nscf_nband
screening = aobj.Screening(ecuteps, scr_nband, w_type="RPA", sc_mode="one_shot",
hilbert=None, ecutwfn=None, inclvkb=inclvkb)

multi[2].set_vars(nscf_ksampling.to_abivars())
multi[2].set_vars(nscf_electrons.to_abivars())
multi[2].set_vars(screening.to_abivars())
multi[2].set_vars(_stopping_criterion("screening", accuracy)) # Dummy

# Sigma.
if sigma_nband is None: sigma_nband = nscf_nband
self_energy = aobj.SelfEnergy("gw", "one_shot", sigma_nband, ecutsigx, screening,
gw_qprange=gw_qprange, ppmodel=ppmodel)

multi[3].set_vars(nscf_ksampling.to_abivars())
multi[3].set_vars(nscf_electrons.to_abivars())
multi[3].set_vars(self_energy.to_abivars())
multi[3].set_vars(_stopping_criterion("sigma", accuracy)) # Dummy

# TODO: Cannot use istwfk != 1.
multi.set_vars(istwfk="*1")

return multi
nband = _find_scf_nband(structure, scf_inp.pseudos, scf_electrons)
scf_inp.set_vars(nband=nband)
# Non-Scf input
nscf_inp = nscf_from_gsinput(gs_input=scf_inp, nband=nscf_nband, accuracy=accuracy)
# Scr input
scr_inp = scr_from_nscfinput(nscf_input=nscf_inp, nband=scr_nband, ecuteps=ecuteps, ecutwfn=None,
inclvkb=inclvkb, w_type="RPA", sc_mode="one_shot", hilbert=None, accuracy="normal")
# Sigma input
sigma_inp = sigma_from_inputs(nscf_input=nscf_inp, scr_input=scr_inp, nband=sigma_nband,
ecutwfn=None, ecuteps=None, ecutsigx=ecutsigx,
ppmodel=ppmodel, gw_qprange=gw_qprange)
return MultiDataset.from_inputs([scf_inp, nscf_inp, scr_inp, sigma_inp])


def g0w0_convergence_inputs(structure, pseudos, kppa, nscf_nband, ecuteps, ecutsigx, scf_nband, ecut,
Expand Down Expand Up @@ -1118,15 +1161,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|
"""
Expand All @@ -1137,29 +1183,111 @@ 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)
bands_input.set_vars(_stopping_criterion("nscf", accuracy))

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


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()

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)

nscf_input.set_vars(nband=nband, iscf=-2)
nscf_input.set_vars(_stopping_criterion("nscf", accuracy))

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)
dos_input.set_vars(_stopping_criterion("nscf", accuracy))
return nscf_input

if pdos:
# FIXME
raise NotImplementedError()

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 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:
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 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:
ValueError(f"Unrecognized value for projection: {projection}")

if projection is not None and "m" in projection.lower():
dos_input.set_vars(prtdosm=1)

return dos_input

Expand Down Expand Up @@ -1219,13 +1347,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)
Expand Down
Loading

0 comments on commit f8cae90

Please sign in to comment.