4. FCC Si density of states

Let us continue with the VASP tutorials. In particular, let us now extract the density of states for FCC Si as done in the second tutorial. We will from here and in the rest of the tutorials not run VASP in the regular manner and leave those exercise to the reader. Please do so, to appreciate how much simpler AiiDA-VASP makes these operations, in particular when you need to perform many or repeat calculations.

Again, we assume you have completed the previous tutorial.

As you might have noticed we have now populated the database with the silicon structure multiple times. Let us try instead to load some of the structures that already are present to save coding, reuse previous results and save data storage.

In previous tutorial, when developing the call script, we set the label. Please inspect this. E.g. for a lattice constant of 3.9, the label is set to silicon_at_3_9. Remember also that we cannot modify structures already present in the database. If we want to do this, we need to create a new Structure Data using the stored entry as an initializer. In this case we will not modify the structure, so we will simply use the unmodified one.

  1. Let us first create a simple call script that uses the VASP workchain:

    Call script to calculate the total energies for one volume of standard silicon.
    This particular call script set up a standard calculation that execute a calculation for
    the fcc silicon structure.
    # pylint: disable=too-many-arguments, invalid-name, import-outside-toplevel
    from aiida.common.extendeddicts import AttributeDict
    from aiida.orm import Code, Bool, Str
    from aiida.plugins import DataFactory, WorkflowFactory
    from aiida.engine import submit
    from aiida import load_profile
    def get_structure(label):
        from aiida.orm import QueryBuilder
        qb = QueryBuilder()
        qb.append(DataFactory('structure'), filters={'label': {'==': label}}, tag='structure')
        # Pick any structure with this label, here, just the first
        return qb.all()[0][0]
    def main(code_string, incar, kmesh, structure, potential_family, potential_mapping, options):
        """Main method to setup the calculation."""
        # First, we need to fetch the AiiDA datatypes which will
        # house the inputs to our calculation
        dict_data = DataFactory('dict')
        kpoints_data = DataFactory('array.kpoints')
        # Then, we set the workchain you would like to call
        workchain = WorkflowFactory('vasp.verify')
        # And finally, we declare the options, settings and input containers
        settings = AttributeDict()
        inputs = AttributeDict()
        # organize settings
        settings.parser_settings = {'add_dos': True}
        # set inputs for the following WorkChain execution
        # set code
        inputs.code = Code.get_from_string(code_string)
        # set structure
        inputs.structure = structure
        # set k-points grid density
        kpoints = kpoints_data()
        inputs.kpoints = kpoints
        # set parameters
        inputs.parameters = dict_data(dict=incar)
        # set potentials and their mapping
        inputs.potential_family = Str(potential_family)
        inputs.potential_mapping = dict_data(dict=potential_mapping)
        # set options
        inputs.options = dict_data(dict=options)
        # set settings
        inputs.settings = dict_data(dict=settings)
        # set workchain related inputs, in this case, give more explicit output to report
        inputs.verbose = Bool(True)
        # submit the requested workchain with the supplied inputs
        submit(workchain, **inputs)
    if __name__ == '__main__':
        # Code_string is chosen among the list given by 'verdi code list'
        CODE_STRING = 'vasp@mycluster'
        # INCAR equivalent
        # Set input parameters
        INCAR = {'encut': 240, 'ismear': -5, 'lorbit': 11}
        # KPOINTS equivalent
        # Set kpoint mesh
        KMESH = [21, 21, 21]
        # POTCAR equivalent
        # Potential_family is chosen among the list given by
        # 'verdi data vasp-potcar listfamilies'
        POTENTIAL_FAMILY = 'pbe'
        # The potential mapping selects which potential to use, here we use the standard
        # for silicon, this could for instance be {'Si': 'Si_GW'} to use the GW ready
        # potential instead
        POTENTIAL_MAPPING = {'Si': 'Si'}
        # jobfile equivalent
        # In options, we typically set scheduler options.
        # See https://aiida.readthedocs.io/projects/aiida-core/en/latest/scheduler/index.html
        # AttributeDict is just a special dictionary with the extra benefit that
        # you can set and get the key contents with mydict.mykey, instead of mydict['mykey']
        OPTIONS = AttributeDict()
        OPTIONS.account = ''
        OPTIONS.qos = ''
        OPTIONS.resources = {'num_machines': 1, 'num_mpiprocs_per_machine': 16}
        OPTIONS.queue_name = ''
        OPTIONS.max_wallclock_seconds = 3600
        OPTIONS.max_memory_kb = 1024000
        # POSCAR equivalent
        # Set the silicon structure
        STRUCTURE = get_structure('silicon_at_3_9')

    Or you can download it:

    %/$ wget https://github.com/aiida-vasp/aiida-vasp/raw/develop/tutorials/run_fcc_si_dos.py

    Where we have modified the input parameters according to the tutorial density of states for FCC Si. And importantly, we have also told the parser to give us the density of states.

  2. And execute it:

    %/$ python run_fcc_si_dos.py
  3. After a while we check the status:

    %/$ verdi process list -a
    PK  Created    Process label                 Process State      Process status
    ------  ---------  ----------------------------  -----------------  --------------------------------------
    103820  2m ago     VerifyWorkChain               ⏹ Finished [0]
    103821  2m ago     VaspWorkChain                 ⏹ Finished [0]
    103822  2m ago     VaspCalculation               ⏹ Finished [0]
  4. Check the content of the topmost workchain:

    %/$ verdi process show 103820
    (aiida) [efl@efl tutorials]$ verdi process show 103820
    Property       Value
    -------------  ------------------------------------
    type           WorkChainNode
    pk             103820
    uuid           ae0047fe-6351-4585-a60c-ad02dcda8f93
    ctime          2019-10-08 14:30:46.965520+00:00
    mtime          2019-10-08 14:31:53.971080+00:00
    process state  Finished
    exit status    0
    computer       [6] mycluster
    Inputs                     PK  Type
    ---------------------  ------  -------------
    clean_workdir          103818  Bool
    code                   101271  Code
    kpoints                103810  KpointsData
    max_iterations         103817  Int
    options                103814  Dict
    parameters             103811  Dict
    potential_family       103812  Str
    potential_mapping      103813  Dict
    settings               103815  Dict
    structure              103729  StructureData
    verbose                103816  Bool
    verify_max_iterations  103819  Int
    Outputs            PK  Type
    -------------  ------  ----------
    dos            103825  ArrayData
    misc           103826  Dict
    remote_folder  103823  RemoteData
    retrieved      103824  FolderData
    Called        PK  Type
    --------  ------  -------------
    CALL      103821  WorkChainNode
    Log messages
    There are 1 log messages for this calculation
    Run 'verdi process report 103820' to see them

    And as you can see, dos is now listed in the output.

    Now, as you may already know, running with such a dense k-point grid for the initial calculation is usually not a good idea. It is more efficient to pre-converge the electronic states using a more sparse k-point grid and then restart the calculation using a more dense k-point grid when calculating the density of states.