Yambo-AiiDA Tutorial¶
The following tutorial shows how to run all the necessary steps to obtain a G0W0 result with Yambo for bulk hBN. In order to keep the tutorial light in terms of computational resources and time of execution, calculations are not fully converged with respect to parameters such as k-points, empty bands or G-vectors.
SCF step (Quantum ESPRESSO)¶
Using the AiiDA quantumespresso.pw plugin, we begin with submitting an SCF calculation.
We are going to use the pk
of the SCF calculation in the next steps.
We use the PwBaseWorkChain to submit a pw calculation, in such a way to have automatic
error handling and restarting from failed runs.
For details on how to use the quantumespresso.pw plugin, please refer to the plugins documentation page. Remember to replace the codename and pseudo-family with those configured in your AiiDA installation. NB: Yambo can be used only with norm-conserving pseudopotentials!
# -*- coding: utf-8 -*-
from __future__ import absolute_import
from __future__ import print_function
import sys
import os
from aiida.plugins import DataFactory
from aiida.engine import submit
from aiida_quantumespresso.utils.pseudopotential import validate_and_prepare_pseudos_inputs
from aiida_quantumespresso.workflows.pw.base import PwBaseWorkChain
from ase import Atoms
import argparse
def get_options():
parser = argparse.ArgumentParser(description='SCF calculation.')
parser.add_argument(
'--code',
type=int,
dest='code_pk',
required=True,
help='The pw codename to use')
parser.add_argument(
'--pseudo',
type=str,
dest='pseudo_family',
required=True,
help='The pseudo family to use')
parser.add_argument(
'--time',
type=int,
dest='max_wallclock_seconds',
required=False,
default=30*60,
help='max wallclock in seconds')
parser.add_argument(
'--nodes',
type=int,
dest='num_machines',
required=False,
default=1,
help='number of machines')
parser.add_argument(
'--mpi',
type=int,
dest='num_mpiprocs_per_machine',
required=False,
default=1,
help='number of mpi processes per machine')
parser.add_argument(
'--threads',
type=int,
dest='num_cores_per_mpiproc',
required=False,
default=1,
help='number of threads per mpi process')
parser.add_argument(
'--queue_name',
type=str,
dest='queue_name',
required=False,
default=None,
help='queue(PBS) or partition(SLURM) name')
parser.add_argument(
'--qos',
type=str,
dest='qos',
required=False,
default=None,
help='qos name')
parser.add_argument(
'--account',
type=str,
dest='account',
required=False,
default=None,
help='account name')
args = parser.parse_args()
###### setting the machine options ######
options = {
'code_pk': args.code_pk,
'pseudo_family': args.pseudo_family,
'max_wallclock_seconds': args.max_wallclock_seconds,
'resources': {
"num_machines": args.num_machines,
"num_mpiprocs_per_machine": args.num_mpiprocs_per_machine,
"num_cores_per_mpiproc": args.num_cores_per_mpiproc,
},
'custom_scheduler_commands': u"export OMP_NUM_THREADS="+str(args.num_cores_per_mpiproc),
}
if args.queue_name:
options['queue_name']=args.queue_name
if args.qos:
options['qos']=args.qos
if args.account:
options['account']=args.account
return options
def main(options):
###### setting the lattice structure ######
alat = 2.4955987320 # Angstrom
the_cell = [[1.000000*alat, 0.000000, 0.000000],
[-0.500000*alat, 0.866025*alat, 0.000000],
[0.000000, 0.000000, 6.4436359260]]
atoms = Atoms('BNNB', [(1.2477994910, 0.7204172280, 0.0000000000),
(-0.0000001250, 1.4408346720, 0.0000000000),
(1.2477994910, 0.7204172280, 3.2218179630),
(-0.0000001250,1.4408346720, 3.2218179630)],
cell = [1,1,1])
atoms.set_cell(the_cell, scale_atoms=False)
atoms.set_pbc([True,True,True])
StructureData = DataFactory('structure')
structure = StructureData(ase=atoms)
###### setting the kpoints mesh ######
KpointsData = DataFactory('array.kpoints')
kpoints = KpointsData()
kpoints.set_kpoints_mesh([6,6,2])
###### setting the scf parameters ######
Dict = DataFactory('dict')
params_scf = {
'CONTROL': {
'calculation': 'scf',
'verbosity': 'high',
'wf_collect': True
},
'SYSTEM': {
'ecutwfc': 130.,
'force_symmorphic': True,
'nbnd': 20
},
'ELECTRONS': {
'mixing_mode': 'plain',
'mixing_beta': 0.7,
'conv_thr': 1.e-8,
'diago_thr_init': 5.0e-6,
'diago_full_acc': True
},
}
parameter_scf = Dict(dict=params_scf)
###### creation of the workchain ######
builder = PwBaseWorkChain.get_builder()
builder.pw.structure = structure
builder.pw.parameters = parameter_scf
builder.kpoints = kpoints
builder.pw.metadata.options.max_wallclock_seconds = \
options['max_wallclock_seconds']
builder.pw.metadata.options.resources = \
dict = options['resources']
if 'queue_name' in options:
builder.pw.metadata.options.queue_name = options['queue_name']
if 'qos' in options:
builder.pw.metadata.options.qos = options['qos']
if 'account' in options:
builder.metadata.options.account = options['account']
builder.pw.metadata.options.custom_scheduler_commands = options['custom_scheduler_commands']
builder.pw.code = load_node(options['code_pk'])
builder.pw.pseudos = validate_and_prepare_pseudos_inputs(
builder.pw.structure, pseudo_family = Str(options['pseudo_family']))
return builder
if __name__ == "__main__":
options = get_options()
builder = main(options)
running = submit(builder)
print("Submitted PwBaseWorkchain; with pk=<{}>".format(running.pk))
As you can notice, we use the “argparse” module to provide some inputs from the shell, like the code and the pseudos to be used in the simulation. In practice, the command to be run would be like:
verdi run name_of_the_script.py --code <pk of the pw code> --pseudo <name of the pseudofamily>
NSCF step (Quantum ESPRESSO) for G0W0¶
Using the pk
of the SCF calculation, we now run a NSCF calculation as the starting point for the GW calculation.
# -*- coding: utf-8 -*-
from __future__ import absolute_import
from __future__ import print_function
import sys
import os
from aiida.plugins import DataFactory
from aiida.engine import submit
from aiida_quantumespresso.utils.pseudopotential import validate_and_prepare_pseudos_inputs
from aiida_quantumespresso.workflows.pw.base import PwBaseWorkChain
from ase import Atoms
import argparse
def get_options():
parser = argparse.ArgumentParser(description='SCF calculation.')
parser.add_argument(
'--code',
type=int,
dest='code_pk',
required=True,
help='The pw codename to use')
parser.add_argument(
'--pseudo',
type=str,
dest='pseudo_family',
required=True,
help='The pseudo family to use')
parser.add_argument(
'--time',
type=int,
dest='max_wallclock_seconds',
required=False,
default=30*60,
help='max wallclock in seconds')
parser.add_argument(
'--nodes',
type=int,
dest='num_machines',
required=False,
default=1,
help='number of machines')
parser.add_argument(
'--mpi',
type=int,
dest='num_mpiprocs_per_machine',
required=False,
default=1,
help='number of mpi processes per machine')
parser.add_argument(
'--threads',
type=int,
dest='num_cores_per_mpiproc',
required=False,
default=1,
help='number of threads per mpi process')
parser.add_argument(
'--queue_name',
type=str,
dest='queue_name',
required=False,
default=None,
help='queue(PBS) or partition(SLURM) name')
parser.add_argument(
'--qos',
type=str,
dest='qos',
required=False,
default=None,
help='qos name')
parser.add_argument(
'--account',
type=str,
dest='account',
required=False,
default=None,
help='account name')
args = parser.parse_args()
###### setting the machine options ######
options = {
'code_pk': args.code_pk,
'pseudo_family': args.pseudo_family,
'max_wallclock_seconds': args.max_wallclock_seconds,
'resources': {
"num_machines": args.num_machines,
"num_mpiprocs_per_machine": args.num_mpiprocs_per_machine,
"num_cores_per_mpiproc": args.num_cores_per_mpiproc,
},
'custom_scheduler_commands': u"export OMP_NUM_THREADS="+str(args.num_cores_per_mpiproc),
}
if args.queue_name:
options['queue_name']=args.queue_name
if args.qos:
options['qos']=args.qos
if args.account:
options['account']=args.account
return options
def main(options):
###### setting the lattice structure ######
alat = 2.4955987320 # Angstrom
the_cell = [[1.000000*alat, 0.000000, 0.000000],
[-0.500000*alat, 0.866025*alat, 0.000000],
[0.000000, 0.000000, 6.4436359260]]
atoms = Atoms('BNNB', [(1.2477994910, 0.7204172280, 0.0000000000),
(-0.0000001250, 1.4408346720, 0.0000000000),
(1.2477994910, 0.7204172280, 3.2218179630),
(-0.0000001250,1.4408346720, 3.2218179630)],
cell = [1,1,1])
atoms.set_cell(the_cell, scale_atoms=False)
atoms.set_pbc([True,True,True])
StructureData = DataFactory('structure')
structure = StructureData(ase=atoms)
###### setting the kpoints mesh ######
KpointsData = DataFactory('array.kpoints')
kpoints = KpointsData()
kpoints.set_kpoints_mesh([6,6,2])
###### setting the scf parameters ######
Dict = DataFactory('dict')
params_nscf = {
'CONTROL': {
'calculation': 'nscf',
'verbosity': 'high',
'wf_collect': True
},
'SYSTEM': {
'ecutwfc': 130.,
'force_symmorphic': True,
'nbnd': 150
},
'ELECTRONS': {
'mixing_mode': 'plain',
'mixing_beta': 0.7,
'conv_thr': 1.e-8,
'diago_thr_init': 5.0e-6,
'diago_full_acc': True
},
}
parameter_nscf = Dict(dict=params_nscf)
###### creation of the workchain ######
builder = PwBaseWorkChain.get_builder()
builder.pw.structure = structure
builder.pw.parameters = parameter_nscf
builder.kpoints = kpoints
builder.pw.metadata.options.max_wallclock_seconds = \
options['max_wallclock_seconds']
builder.pw.metadata.options.resources = \
dict = options['resources']
if 'queue_name' in options:
builder.pw.metadata.options.queue_name = options['queue_name']
if 'qos' in options:
builder.pw.metadata.options.qos = options['qos']
if 'account' in options:
builder.metadata.options.account = options['account']
builder.pw.metadata.options.custom_scheduler_commands = options['custom_scheduler_commands']
builder.pw.code = load_node(options['code_pk'])
builder.pw.pseudos = validate_and_prepare_pseudos_inputs(
builder.pw.structure, pseudo_family = Str(options['pseudo_family']))
return builder
if __name__ == "__main__":
options = get_options()
builder = main(options)
running = submit(builder)
print("Submitted PwBaseWorkchain; with pk=<{}>".format(running.pk))
P2Y step (Yambo)¶
Now we use the Yambo plugin to run the p2y code, converting the Quantum ESPRESSO files into a NetCDF Yambo database.
# -*- coding: utf-8 -*-
from __future__ import absolute_import
from __future__ import print_function
import sys
import os
from aiida.plugins import DataFactory, CalculationFactory
from aiida.engine import submit
from aiida_yambo.calculations.gw import YamboCalculation
import argparse
def get_options():
parser = argparse.ArgumentParser(description='SCF calculation.')
parser.add_argument(
'--code',
type=int,
dest='code_pk',
required=True,
help='The yambo codename to use')
parser.add_argument(
'--precode',
type=int,
dest='precode_pk',
required=True,
help='The p2y codename to use')
parser.add_argument(
'--parent',
type=int,
dest='parent_pk',
required=True,
help='The parent to use')
parser.add_argument(
'--time',
type=int,
dest='max_wallclock_seconds',
required=False,
default=30*60,
help='max wallclock in seconds')
parser.add_argument(
'--nodes',
type=int,
dest='num_machines',
required=False,
default=1,
help='number of machines')
parser.add_argument(
'--mpi',
type=int,
dest='num_mpiprocs_per_machine',
required=False,
default=1,
help='number of mpi processes per machine')
parser.add_argument(
'--threads',
type=int,
dest='num_cores_per_mpiproc',
required=False,
default=1,
help='number of threads per mpi process')
parser.add_argument(
'--queue_name',
type=str,
dest='queue_name',
required=False,
default=None,
help='queue(PBS) or partition(SLURM) name')
parser.add_argument(
'--qos',
type=str,
dest='qos',
required=False,
default=None,
help='qos name')
parser.add_argument(
'--account',
type=str,
dest='account',
required=False,
default=None,
help='account name')
args = parser.parse_args()
###### setting the machine options ######
options = {
'code_pk': args.code_pk,
'precode_pk': args.precode_pk,
'parent_pk': args.parent_pk,
'max_wallclock_seconds': args.max_wallclock_seconds,
'resources': {
"num_machines": args.num_machines,
"num_mpiprocs_per_machine": args.num_mpiprocs_per_machine,
"num_cores_per_mpiproc": args.num_cores_per_mpiproc,
},
'custom_scheduler_commands': u"export OMP_NUM_THREADS="+str(args.num_cores_per_mpiproc),
}
if args.queue_name:
options['queue_name']=args.queue_name
if args.qos:
options['qos']=args.qos
if args.account:
options['account']=args.account
return options
def main(options):
Dict = DataFactory('dict')
###### creation of the YamboCalculation ######
builder = YamboCalculation.get_builder()
builder.metadata.options.max_wallclock_seconds = \
options['max_wallclock_seconds']
builder.metadata.options.resources = \
dict = options['resources']
if 'queue_name' in options:
builder.metadata.options.queue_name = options['queue_name']
if 'qos' in options:
builder.metadata.options.qos = options['qos']
if 'account' in options:
builder.metadata.options.account = options['account']
builder.metadata.options.custom_scheduler_commands = options['custom_scheduler_commands']
builder.parameters = Dict(dict={})
builder.precode_parameters = Dict(dict={})
builder.settings = Dict(dict={'INITIALISE': True, 'PARENT_DB': False})
builder.code = load_node(options['code_pk'])
builder.preprocessing_code = load_node(options['precode_pk'])
builder.parent_folder = load_node(options['parent_pk']).outputs.remote_folder
return builder
if __name__ == "__main__":
options = get_options()
builder = main(options)
running = submit(builder)
print("Submitted YamboCalculation; with pk=<{}>".format(running.pk))
G0W0 (Yambo)¶
Now we are ready to run a G0W0 calculation in the plasmon-pole approximation (PPA), in particular we compute the direct band gap at Gamma of hBN.
# -*- coding: utf-8 -*-
from __future__ import absolute_import
from __future__ import print_function
import sys
import os
from aiida.plugins import DataFactory, CalculationFactory
from aiida.engine import submit
from aiida_yambo.calculations.gw import YamboCalculation
import argparse
def get_options():
parser = argparse.ArgumentParser(description='SCF calculation.')
parser.add_argument(
'--code',
type=int,
dest='code_pk',
required=True,
help='The yambo codename to use')
parser.add_argument(
'--precode',
type=int,
dest='precode_pk',
required=True,
help='The p2y codename to use')
parser.add_argument(
'--parent',
type=int,
dest='parent_pk',
required=True,
help='The parent to use')
parser.add_argument(
'--time',
type=int,
dest='max_wallclock_seconds',
required=False,
default=30*60,
help='max wallclock in seconds')
parser.add_argument(
'--nodes',
type=int,
dest='num_machines',
required=False,
default=1,
help='number of machines')
parser.add_argument(
'--mpi',
type=int,
dest='num_mpiprocs_per_machine',
required=False,
default=1,
help='number of mpi processes per machine')
parser.add_argument(
'--threads',
type=int,
dest='num_cores_per_mpiproc',
required=False,
default=1,
help='number of threads per mpi process')
parser.add_argument(
'--queue_name',
type=str,
dest='queue_name',
required=False,
default=None,
help='queue(PBS) or partition(SLURM) name')
parser.add_argument(
'--qos',
type=str,
dest='qos',
required=False,
default=None,
help='qos name')
parser.add_argument(
'--account',
type=str,
dest='account',
required=False,
default=None,
help='account name')
args = parser.parse_args()
###### setting the machine options ######
options = {
'code_pk': args.code_pk,
'precode_pk': args.precode_pk,
'parent_pk': args.parent_pk,
'max_wallclock_seconds': args.max_wallclock_seconds,
'resources': {
"num_machines": args.num_machines,
"num_mpiprocs_per_machine": args.num_mpiprocs_per_machine,
"num_cores_per_mpiproc": args.num_cores_per_mpiproc,
},
'custom_scheduler_commands': u"export OMP_NUM_THREADS="+str(args.num_cores_per_mpiproc),
}
if args.queue_name:
options['queue_name']=args.queue_name
if args.qos:
options['qos']=args.qos
if args.account:
options['account']=args.account
return options
def main(options):
###### setting the gw parameters ######
Dict = DataFactory('dict')
params_gw = {
'ppa': True,
'gw0': True,
'HF_and_locXC': True,
'em1d': True,
'Chimod': 'hartree',
#'EXXRLvcs': 40,
#'EXXRLvcs_units': 'Ry',
'BndsRnXp': [1, 10],
'NGsBlkXp': 2,
'NGsBlkXp_units': 'Ry',
'GbndRnge': [1, 10],
'DysSolver': "n",
'QPkrange': [[1, 1, 8, 9]],
'X_all_q_CPU': "1 1 1 1",
'X_all_q_ROLEs': "q k c v",
'SE_CPU': "1 1 1",
'SE_ROLEs': "q qp b",
}
params_gw = Dict(dict=params_gw)
###### creation of the YamboCalculation ######
builder = YamboCalculation.get_builder()
builder.metadata.options.max_wallclock_seconds = \
options['max_wallclock_seconds']
builder.metadata.options.resources = \
dict = options['resources']
if 'queue_name' in options:
builder.metadata.options.queue_name = options['queue_name']
if 'qos' in options:
builder.metadata.options.qos = options['qos']
if 'account' in options:
builder.metadata.options.account = options['account']
builder.metadata.options.custom_scheduler_commands = options['custom_scheduler_commands']
builder.parameters = params_gw
builder.precode_parameters = Dict(dict={})
builder.settings = Dict(dict={'INITIALISE': False, 'PARENT_DB': False})
builder.code = load_node(options['code_pk'])
builder.preprocessing_code = load_node(options['precode_pk'])
builder.parent_folder = load_node(options['parent_pk']).outputs.remote_folder
return builder
if __name__ == "__main__":
options = get_options()
builder = main(options)
running = submit(builder)
print("Submitted YamboCalculation; with pk=<{}>".format(running.pk))
The quasiparticle corrections and the renormalization factors can be accessed from the Yambo calculation (yambo_calc) using the output bands and array data:
yambo_calc = load_node(pk)
energies_DFT = yambo_calc.outputs.array_ndb.get_array('E_0')
QP_corrections = yambo_calc.outputs.array_ndb.get_array('E_minus_Eo')
Z_factors = yambo_calc.outputs.array_ndb.get_array('Z')
kpoint_band_array = yambo_calc.outputs.array_ndb.get_array('qp_table')
kpoints = y.outputs.bands_quasiparticle.get_kpoints()
To retrieve additional files:
settings = ParameterData(dict={"ADDITIONAL_RETRIEVE_LIST":['r-*','o-*','LOG/l-*01',
'aiida.out/ndb.QP','aiida.out/ndb.HF_and_locXC']})
builder.use_settings(settings)
This selects the additional files that will be retrieved and parsed after a calculation. Supported
files include the report files r-*
, text outputs o-*
, logs, the quasiparticle
database for GW calculations aiida.out/ndb.QP
, and the Hartree-Fock and local exchange
db aiida.out/ndb.HF_and_locXC
.