First steps into MBPT calculations

The following tutorial shows how to run all the necessary steps to obtain a the databases needed to start a Yambo calculation. The starting point is a self-consistent calculation of the electronic density, and then a calculation of the electronic wavefunctions through a non-self-consistent DFT calculation. So, the first AiiDA plugin used here is the QuantumEspresso one. Then, it is necessary to convert the results in Yambo-readable format, i.e. NetCDF format: this is done using the p2y executable, included in the yambo package and executed from the same yambo plugin. 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. The example here considers the bulk hexagonal boron nitride hBN.

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. Remember that the pk is the number that identifies the node in the AiiDA database. 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 respective 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>

this is just our choice to build the calculation, it is possible also to do it interactively (jupyter notebooks) or adapting the script to have no input to provide.

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. The parent folder now should be the nscf one.

# -*- 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.yambo 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, 'COPY_DBS': 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 fundamental input that tells the plugin to only run a p2y calculation is the settings key INITALISE: if True, the only code running will be p2y. If False, an automatic procedure is implemented to decide what executables are needed to complete the calculation.