Merge pull request #221 from gbrunin/eph

Added workflow for e-ph computations of the mobility in semiconductors
This commit is contained in:
Matteo Giantomassi 2021-05-22 23:07:17 +02:00 committed by GitHub
commit 52d7abe24d
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
6 changed files with 290 additions and 50 deletions

View File

@ -12,6 +12,7 @@ from monty.json import jsanitize, MontyDecoder, MSONable
from pymatgen.util.serialization import pmg_serialize
from abipy.core.structure import Structure
from abipy.abio.inputs import AbinitInput, MultiDataset
import abipy.core.abinit_units as abu
__all__ = [
@ -1459,7 +1460,8 @@ def conduc_from_inputs(scf_input, nscf_input, tmesh, ddb_ngqpt, eph_ngqpt_fine,
def conduc_kerange_from_inputs(scf_input, nscf_input, tmesh, ddb_ngqpt, eph_ngqpt_fine,
sigma_ngkpt, sigma_erange, einterp=(1, 5, 0, 0), boxcutmin=1.1, mixprec=1):
sigma_ngkpt, sigma_erange, sigma_kerange=None, epad=0.25*abu.eV_Ha,
einterp=(1, 5, 0, 0), boxcutmin=1.1, mixprec=1):
"""
Returns a list of inputs in the form of a MultiDataset to perform a set of calculations to determine the conductivity.
This part require a ground state |AbinitInput| and a non self-consistent |AbinitInput|. You will also need
@ -1473,6 +1475,9 @@ def conduc_kerange_from_inputs(scf_input, nscf_input, tmesh, ddb_ngqpt, eph_ngqp
eph_ngqpt_fine: the fine qpoints grid that will be interpolated.
sigma_ngkpt: The fine grid of kpt inside the sigma interval
sigma_erange: The energy range for Sigma_nk
sigma_kerange: The energy window for the WFK generation (should be larger than sigma_erange). Can be specified directly
or determined from epad, if sigma_kerange is None.
epad: Additional energy range to sigma_erange to determine sigma_kerange automatically.
einterp: The interpolation used. By default it is a star-function interpolation.
boxcutmin: For the last task only, 1.1 is often used to decrease memory and is faster over the Abinit default of 2.
mixprec: For the last task only, 1 is often used to make the EPH calculation faster. Note that Abinit default is 0.
@ -1484,9 +1489,24 @@ def conduc_kerange_from_inputs(scf_input, nscf_input, tmesh, ddb_ngqpt, eph_ngqp
extension = MultiDataset.replicate_input(nscf_input, 4)
multi.extend(extension)
# If sigma_kerange is None, we determine it automatically based on epad
# We have to consider semiconductors and metals, electrons and holes
if sigma_kerange is None:
h_range = sigma_erange[0]
e_range = sigma_erange[1]
sigma_kerange = [h_range, e_range]
if h_range < 0:
sigma_kerange[0] -= epad
if e_range < 0:
sigma_kerange[1] -= epad
if h_range > 0:
sigma_kerange[0] += epad
if e_range > 0:
sigma_kerange[1] += epad
# Modify the second nscf input to get a task that calculate the kpt in the sigma interval (Kerange.nc file)
multi[2].set_vars(optdriver=8, wfk_task='"wfk_kpts_erange"', kptopt=1,
sigma_ngkpt=sigma_ngkpt, einterp=einterp, sigma_erange=sigma_erange)
sigma_ngkpt=sigma_ngkpt, einterp=einterp, sigma_erange=sigma_kerange)
# Modify the third nscf input to get a task that add the kpt of Kerange.nc to the WFK file
multi[3].set_vars(optdriver=0, iscf=-2, kptopt=0, ddb_ngqpt=ddb_ngqpt)

View File

@ -1857,6 +1857,73 @@ with the Abinit version you are using. Please contact the AbiPy developers.""" %
return multi
def make_wfk_kerange_input(self, sigma_kerange, sigma_ngkpt, einterp=(1, 5, 0, 0)):
"""
Return |MultiDataset| inputs for a WFK calculation using the kerange trick.
This method should be called with the input associated to the NSCF run that produces
the WFK file used for the interpolation
Args:
sigma_kerange: The energy window for the WFK generation.
sigma_ngkpt: The fine grid of kpt inside the sigma interval.
einterp: The interpolation used. By default it is a star-function interpolation.
"""
# Create a MultiDataset from the nscf input
multi = MultiDataset.from_inputs([self, self])
# Modify the first nscf input to get a task that calculate the kpt in the sigma interval (Kerange.nc file)
multi[0].set_vars(optdriver=8, wfk_task='"wfk_kpts_erange"', kptopt=1,
sigma_ngkpt=sigma_ngkpt, sigma_nshiftk=1, sigma_shiftk=[0,0,0],
einterp=einterp, sigma_erange=sigma_kerange, prtwf=0)
multi[0].pop_vars(["iscf","tolwfr","prtden"])
# Modify the third nscf input to get a task that add the kpt of Kerange.nc to the WFK file
multi[1].set_vars(optdriver=0, iscf=-2, kptopt=0, ngkpt=sigma_ngkpt)
return multi
def make_eph_transport_input(self, ddb_ngqpt, sigma_erange, tmesh, eph_ngqpt_fine=None,
mixprec=1, boxcutmin=1.1, ibte_prep=0, ibte_niter=200, ibte_abs_tol=1e-3):
"""
Return an |AbinitInput| to perform phonon-limited transport calculations.
This method is usually called with with the input associated to the NSCF run that produces
the WFK file used to start the EPH run so that we can directly inherit the k-mesh
Args:
ddb_ngqpt: the coarse qpt grid used to compute the DDB and DVDB files in the phonon_work.
eph_ngqpt_fine: the fine qpt grid used for the Fourier interpolation.
sigma_erange: Energy window for k-states (see Abinit variable)
tmesh: The mesh of temperatures (in Kelvin)
boxcutmin: For the last task only, 1.1 is often used to decrease memory and is faster over the Abinit default of 2.
mixprec: For the last task only, 1 is often used to make the EPH calculation faster. Note that Abinit default is 0.
ibte_prep: Set it to 1 to activate the iterative Boltzmann equation. Default is RTA.
ibte_niter: Number of iterations to solve the Boltzmann equation.
ibte_abs_tol: Stopping criterion for the IBTE.
"""
eph_ngqpt_fine = self.get("ngkpt") if eph_ngqpt_fine is None else eph_ngqpt_fine
nbdbuf = 0 if self.get("nbdbuf") is None else self.get("nbdbuf")
nband = self.get("nband")-nbdbuf # If nbdbuf is used in the NSCF computation,
# it cannot be used in the EPH driver and should
# be removed
new = self.new_with_vars(
optdriver=7, # Enter EPH driver.
eph_task=-4, # Compute imag part of sigma_eph.
ddb_ngqpt=ddb_ngqpt, # Ab-initio coarse q-mesh used to produce the DDB/DVDB files.
eph_ngqpt_fine=eph_ngqpt_fine, # Interpolate DFPT potentials on this denser q-mesh.
sigma_erange=sigma_erange,
tmesh=tmesh,
mixprec=mixprec,
boxcutmin=boxcutmin,
ibte_prep=ibte_prep,
ibte_niter=ibte_niter,
ibte_abs_tol=ibte_abs_tol,
nband=nband,
)
new.pop_vars(["iscf","prtwf","tolwfr","prtden","nbdbuf","kptopt"])
#new.add_phbbands_vars()
return new
#def make_eph_zpr_input(self, ddb_ngqpt, tmesh, eph_ngqpt_fine=None,,
# mixprec=1, boxcutmin=1.1):
# """
@ -1885,37 +1952,6 @@ with the Abinit version you are using. Please contact the AbiPy developers.""" %
# #new.add_phbbands_vars()
# return ne
#def make_eph_transport_input(self, ddb_ngqpt, sigma_erange, tmesh, eph_ngqpt_fine=None,,
# mixprec=1, boxcutmin=1.1, ibte_prep=0):
# """
# Return an |AbinitInput| to perform phonon-limited transport calculations.
# This method is usually called with with the input associated to the NSCF run that produces
# the WFK file used to start the EPH run so that we can directly inherit the k-mesh
# Args:
# ddb_ngqpt: the coarse qpt grid used to compute the DDB and DVDB files in the phonon_work.
# eph_ngqpt_fine: the fine qpt grid used for the Fourier interpolation.
# sigma_erange: Energy window for k-states (see Abinit variable)
# tmesh: The mesh of temperatures (in Kelvin)
# boxcutmin: For the last task only, 1.1 is often used to decrease memory and is faster over the Abinit default of 2.
# mixprec: For the last task only, 1 is often used to make the EPH calculation faster. Note that Abinit default is 0.
# ibte_prep: Set it to 1 to activate the iterative Boltzmann equation. Default is RTA.
# """
# eph_ngqpt_fine = self.get("ngkpt") if eph_ngqpt_fine is None else eph_ngqpt_fine
# new = self.new_with_vars(
# optdriver=7, # Enter EPH driver.
# eph_task=-4, # Compute imag part of sigma_eph.
# ddb_ngqpt=ddb_ngqpt, # Ab-initio coarse q-mesh used to produce the DDB/DVDB files.
# eph_ngqpt_fine=eph_ngqpt_fine, # Interpolate DFPT potentials on this denser q-mesh.
# sigma_erange=sigma_erange,
# tmesh=tmesh,
# mixprec=mixprec,
# boxcutmin=boxcutmin,
# ibte_prep=ibte_prep,
# )
# #new.add_phbbands_vars()
# return new
#def make_eph_isotc_input(self, ddb_ngqpt, eph_fsewin, eph_ngqpt_fine=None,,
# mixprec=1, boxcutmin=1.1):
# """

View File

@ -625,12 +625,14 @@ class RtaRobot(Robot, RobotWithEbands):
#def get_mobility_mu_dataframe(self, eh=0, component='xx', itemp=0, spin=0, **kwargs):
@add_fig_kwargs
def plot_mobility_kconv(self, eh=0, component='xx', itemp=0, spin=0, fontsize=14, ax=None, **kwargs):
def plot_mobility_kconv(self, eh=0, bte=['serta','mrta','ibte'], component='xx', itemp=0, spin=0, fontsize=14, ax=None, **kwargs):
"""
Plot the convergence of the mobility as a function of the number of k-points.
Plot the convergence of the mobility as a function of the number of k-points,
for different transport formalisms included in the computation.
Args:
eh: 0 for electrons, 1 for holes.
bte: list of transport formalism to plot (serta, mrta, ibte)
component: Cartesian component to plot ('xx', 'xy', ...)
itemp: temperature index.
spin: Spin index.
@ -639,36 +641,42 @@ class RtaRobot(Robot, RobotWithEbands):
Returns: |matplotlib-Figure|
"""
if 'ibte' in bte and not self.all_have_ibte:
print("At least some IBTE results are missing ! Will remove ibte from the bte list")
bte.remove('ibte')
ax, fig, plt = get_ax_fig_plt(ax=ax)
ax.grid(True)
i, j = abu.s2itup(component)
irta = 0
res, temps = []
res, temps = [], []
for ncfile in self.abifiles:
#kptrlattx, kptrlatty, kptrlattz = ncfile.ngkpt
kptrlatt = ncfile.reader.read_value("kptrlatt")
kptrlattx = kptrlatt[0, 0]
kptrlatty = kptrlatt[1, 1]
kptrlattz = kptrlatt[2, 2]
# nctkarr_t('mobility_mu',"dp", "three, three, two, ntemp, nsppol, nrta")]
mobility = ncfile.reader.read_variable("mobility_mu")[irta, spin, itemp, eh, j, i]
#print(mobility)
res.append([kptrlattx, mobility])
mob_serta, mob_mrta, mob_ibte = -1, -1, -1
if 'serta' in bte:
mob_serta = ncfile.reader.read_variable("mobility_mu")[0, spin, itemp, eh, j, i]
if 'mrta' in bte:
mob_mrta = ncfile.reader.read_variable("mobility_mu")[1, spin, itemp, eh, j, i]
if 'ibte' in bte:
mob_ibte = ncfile.reader.read_variable("ibte_mob")[spin, itemp, eh, j, i]
res.append([kptrlattx, mob_serta, mob_mrta, mob_ibte])
temps.append(ncfile.tmesh[itemp])
res.sort(key=lambda t: t[0])
res = np.array(res)
res = np.array(res, dtype=object)
#print(res)
size = 14
ylabel = r"%s mobility (cm$^2$/(V$\cdot$s))" % {0: "Electron", 1: "Hole"}[eh]
ax.set_ylabel(ylabel, size=size)
#if "title" not in kwargs:
# title = r"$\frac{1}{N_k} \sum_{nk} \delta(\epsilon - \epsilon_{nk})$"
# ax.set_title(title, fontsize=fontsize)
from fractions import Fraction
ratio1 = Fraction(kptrlatty, kptrlattx)
ratio2 = Fraction(kptrlattz, kptrlattx)
@ -680,7 +688,15 @@ class RtaRobot(Robot, RobotWithEbands):
ax.set_xlabel(r'Homogeneous $N_k \times$ ' + text1 + r'$N_k \times$ ' + text2 + r'$N_k$ $\mathbf{k}$-point grid',
size=size)
ax.plot(res[:,0], res[:,1], **kwargs)
ax.set_title(component+" component, T = {0} K".format(ncfile.tmesh[itemp]),size=size)
if 'serta' in bte:
ax.plot(res[:,0], res[:,1], '-ob', label='SERTA')
if 'mrta' in bte:
ax.plot(res[:,0], res[:,2], '-or', label='MRTA')
if 'ibte' in bte:
ax.plot(res[:,0], res[:,3], '-og', label='IBTE')
ax.set_xticks(res[:,0].astype(float))
ax.legend(loc="best", shadow=True, fontsize=fontsize)
return fig
@ -773,7 +789,25 @@ class RtaRobot(Robot, RobotWithEbands):
#if self.all_have_ibte:
yield self.plot_ibte_mrta_serta_conv(show=False)
yield self.plot_ibte_vs_rta_rho(show=False)
#self.plot_mobility_kconv(eh=0, component='xx', itemp=0, spin=0, fontsize=14, ax=None, **kwargs):
# Determine the independent component. For the time being,
# only consider the cubic case separately
abifile = self.abifiles[0]
if 'cubic' in abifile.structure.spget_summary():
components = ['xx']
else:
components = ['xx','yy','zz']
# Determine the type of carriers for which the mobility is computed
eh_list = []
if abifile.sigma_erange[0] > 0:
eh_list.append(1)
if abifile.sigma_erange[1] > 0:
eh_list.append(0)
for eh in eh_list:
for comp in components:
yield self.plot_mobility_kconv(eh=eh, component=comp, itemp=0, spin=0, show=False)
#def get_panel(self):
# """
@ -813,7 +847,7 @@ if __name__ == "__main__":
#plt.tick_params(labelsize=14)
#ax = plt.gca()
robot.plot_mobility_kconv(ax=None, color='k', label=r'$N_{{q_{{x,y,z}}}}$ = $N_{{k_{{x,y,z}}}}$')
robot.plot_mobility_kconv(ax=None, color='k')
#fileslist = ['conv_fine/k27x27x27/q27x27x27/Sio_DS1_TRANSPORT.nc',
# 'conv_fine/k30x30x30/q30x30x30/Sio_DS1_TRANSPORT.nc',

View File

@ -0,0 +1,129 @@
#!/usr/bin/env python
r"""
Flow for E-PH calculations
==========================
This flow computes the phonon-limited mobility for different dense k/q meshes in AlAs.
"""
import sys
import os
import abipy.data as abidata
import abipy.abilab as abilab
import abipy.flowtk as flowtk
from abipy.abilab import abiopen
import abipy.core.abinit_units as abu
def build_flow(options):
# Working directory (default is the name of the script with '.py' removed and "run_" replaced by "flow_")
if not options.workdir:
options.workdir = os.path.basename(sys.argv[0]).replace(".py", "").replace("run_", "flow_")
# Initialize the flow
flow = flowtk.Flow(workdir=options.workdir, manager=options.manager)
# Initialize the tmesh, sigma_kerange, sigma_erange and dense meshes for the eph integrations
tmesh = [300,300,1]
sigma_kerange = [0, 0.5*abu.eV_Ha] # 0.5 eV range for the WFK of electrons
sigma_erange = [0, 0.25*abu.eV_Ha] # 0.25 eV range for the EPH computation for electrons
dense_meshes = [[30,30,30],
[40,40,40]]
# Initialize the structure and pseudos
structure = abidata.structure_from_ucell("AlAs")
pseudos = abidata.pseudos("13al.981214.fhi", "33as.pspnc")
# Ground-state computation for 1) the phonons and 2) the WFK generation
scf_input = abilab.AbinitInput(structure, pseudos=pseudos)
scf_input.set_vars(
nband=8,
ecut=2.0,
ngkpt=[4, 4, 4],
nshiftk=1,
shiftk=[0, 0, 0],
tolvrs=1.0e-10,
diemac=9.0,
prtden=1,
#iomode=3,
)
# Initialize the first work and add the ground-state task
work_init = flowtk.Work()
work_init.register_scf_task(scf_input)
# Band structure calculation to make sure everything is ok
# Also allows to compare the results obtained with abitk to
# check the SKW interpolation works as needed
bs_input = scf_input.make_ebands_input(tolwfr=1e-12, ndivsm=10, nb_extra=5)
bs_input.set_vars(nstep=100,nbdbuf=1)
work_init.register_nscf_task(bs_input, deps={work_init[0]: "DEN"})
# NSCF input for the WFK needed to interpolate with kerange
nscf_input = abilab.AbinitInput(structure, pseudos)
nscf_input.set_vars(
ecut=2,
nband=8,
iscf=-2,
tolwfr=1e-20,
prtwf=1,
ngkpt=[16, 16, 16], # Should be dense enough so that the kerange interpolation works
shiftk=[0.0, 0.0, 0.0],
)
work_init.register_nscf_task(nscf_input, deps={work_init[0]: "DEN"})
flow.register_work(work_init)
# Add the phonon work to the flow
ddb_ngqpt = [4, 4, 4]
ph_work = flowtk.PhononWork.from_scf_task(work_init[0], qpoints=ddb_ngqpt, is_ngqpt=True, with_becs=True)
flow.register_work(ph_work)
# We loop over the dense meshes
for i, sigma_ngkpt in enumerate(dense_meshes):
# Use the kerange trick to generate a WFK file
kerange_input, wfk_input = nscf_input.make_wfk_kerange_input(sigma_kerange=sigma_kerange, sigma_ngkpt=sigma_ngkpt).split_datasets()
work_eph = flowtk.Work()
work_eph.register_kerange_task(kerange_input, deps={work_init[2]: "WFK"})
work_eph.register_nscf_task(wfk_input, deps={work_init[0]: "DEN", work_eph[0]: "KERANGE.nc"})
# Generate the input file for the transport calculation
eph_input = wfk_input.make_eph_transport_input(ddb_ngqpt=ddb_ngqpt, sigma_erange=sigma_erange,
tmesh=tmesh, eph_ngqpt_fine=sigma_ngkpt, ibte_prep=1)
# We compute the phonon dispersion to be able to check they are ok
if i==0:
eph_input.set_qpath(20)
work_eph.register_eph_task(eph_input, deps={work_eph[1]: "WFK", ph_work: ["DDB", "DVDB"]})
flow.register_work(work_eph)
flow.allocate(use_smartio=True)
return flow
# 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())

View File

@ -52,6 +52,7 @@ __all__ = [
"ElasticTask",
"SigmaTask",
"EphTask",
"KerangeTask",
"OpticTask",
"AnaddbTask",
"set_user_config_taskmanager",
@ -3940,6 +3941,12 @@ class EphTask(AbinitTask):
color_rgb = np.array((255, 128, 0)) / 255
class KerangeTask(AbinitTask):
"""
Class for kerange calculations.
"""
color_rgb = np.array((255, 128, 0)) / 255
class ManyBodyTask(AbinitTask):
"""
Base class for Many-body tasks (Screening, Sigma, Bethe-Salpeter)

View File

@ -20,7 +20,7 @@ from . import wrappers
from .nodes import Dependency, Node, NodeError, NodeResults, FileNode #, check_spectator
from .tasks import (Task, AbinitTask, ScfTask, NscfTask, DfptTask, PhononTask, ElasticTask, DdkTask, EffMassTask,
BseTask, RelaxTask, DdeTask, BecTask, ScrTask, SigmaTask, TaskManager,
DteTask, EphTask, CollinearThenNonCollinearScfTask)
DteTask, EphTask, KerangeTask, CollinearThenNonCollinearScfTask)
from .utils import Directory
from .netcdf import ETSF_Reader, NetcdfReader
@ -467,6 +467,20 @@ class NodeContainer(metaclass=abc.ABCMeta):
seq_manager = TaskManager.from_user_config().new_with_fixed_mpi_omp(1, 1)
kwargs.update({"manager": seq_manager})
if eph_inp.get("eph_task",0) == -4:
max_cores = TaskManager.from_user_config().qadapter.max_cores
natom3 = 3 * len(eph_inp.structure)
nprocs = max(max_cores - max_cores % natom3, 1)
new_manager = TaskManager.from_user_config().new_with_fixed_mpi_omp(nprocs, 1)
kwargs.update({"manager": new_manager})
return self.register_task(*args, **kwargs)
def register_kerange_task(self, *args, **kwargs):
""" Register a kerange task."""
kwargs["task_class"] = KerangeTask
seq_manager = TaskManager.from_user_config().new_with_fixed_mpi_omp(1, 1)
kwargs.update({"manager": seq_manager})
return self.register_task(*args, **kwargs)
def walknset_vars(self, task_class=None, *args, **kwargs):