5. Parsing and quantities

As you have realized, the plugin does not parse and store many quantities by default. This is done intentionally in order not to fill up disk space without there being a need for it. However, for most users, you would most likely want to have access to more useful quantities after a VASP calculation are finished and this tutorial touches on this to give you a short introduction on how this is done and where to find more information about this. 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 get a feel for how AiiDA-VASP makes these operations more standardized and reproducible, 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 run_fcc_si_workchain, we set the structure.label in get_structure method. Please inspect this now and notice that 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. This goes for all nodes, structure, calculation, workchain or any other node. This is due to the strict compliance of data provenance. If we want to modify a node, for instance, here the case of the StructureData node, 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.
    """
    from aiida import load_profile
    # pylint: disable=too-many-arguments, invalid-name, import-outside-toplevel
    from aiida.common.extendeddicts import AttributeDict
    from aiida.engine import submit
    from aiida.orm import Bool, Code, Str
    from aiida.plugins import DataFactory, WorkflowFactory
    
    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.vasp')
    
        # 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
        # Code
        inputs.code = Code.get_from_string(code_string)
        # Structure
        inputs.structure = structure
        # k-points grid density
        kpoints = kpoints_data()
        kpoints.set_kpoints_mesh(kmesh)
        inputs.kpoints = kpoints
        # Parameters
        inputs.parameters = dict_data(dict=incar)
        # Potential family and their mapping between element and potential type to use
        inputs.potential_family = Str(potential_family)
        inputs.potential_mapping = dict_data(dict=potential_mapping)
        # Options
        inputs.options = dict_data(dict=options)
        # Settings
        inputs.settings = dict_data(dict=settings)
        # Workchain related inputs, in this case, give more explicit output to report
        inputs.verbose = Bool(True)
        # Submit the workchain with the set 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 = {'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.54'
        # 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': 8}
        OPTIONS.queue_name = ''
        OPTIONS.max_wallclock_seconds = 3600
        OPTIONS.max_memory_kb = 2000000
    
        # POSCAR equivalent
        # Set the silicon structure
        STRUCTURE = get_structure('silicon_at_3_9')
    
        main(CODE_STRING, INCAR, KMESH, STRUCTURE, POTENTIAL_FAMILY, POTENTIAL_MAPPING, OPTIONS)
    

    Or you can download it:

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

    Where we have modified the input parameters according to the tutorial density of states for FCC Si.

    Note

    Importantly, we have told the parser to give us the density of states by adding add_dos: True dictionary entry to the parser_settings key in the settings dictionary. In fact, this is the general way to control what quantities to parse, to enable with add_<quantity>: True or add_<quantity>: False, where <quantity> is a given quantity. Please consult Parsing for the supported quantities and additional details.

  2. Let us now execute run_fcc_si_dos.py by issuing the command:

    $ 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
    ----  ---------  --------------------  ----------------  --------------------------------------------------------------------------------------
    2431  6m ago     VaspWorkChain         ⏹ Finished [0]
    2433  6m ago     VaspCalculation       ⏹ Finished [0]
    
  4. Check the content of the VaspWorkChain:

    $ verdi process show 2431
    Property     Value
    -----------  ------------------------------------
    type         VaspWorkChain
    state        Finished [0]
    pk           2431
    uuid         40ce7bd6-cd38-405e-951e-c56251a0cf1b
    label
    description
    ctime        2022-12-22 11:17:25.967623+01:00
    mtime        2022-12-22 11:19:39.816274+01:00
    
    Inputs               PK  Type
    -----------------  ----  -------------
    clean_workdir      2430  Bool
    code                818  InstalledCode
    kpoints            2422  KpointsData
    max_iterations     2429  Int
    options            2426  Dict
    parameters         2423  Dict
    potential_family   2424  Str
    potential_mapping  2425  Dict
    settings           2427  Dict
    structure          1529  StructureData
    verbose            2428  Bool
    
    Outputs          PK  Type
    -------------  ----  ----------
    dos            2436  ArrayData
    misc           2437  Dict
    remote_folder  2434  RemoteData
    retrieved      2435  FolderData
    
    Called          PK  Type
    ------------  ----  ---------------
    iteration_01  2433  VaspCalculation
    
    Log messages
    ---------------------------------------------
    There are 3 log messages for this calculation
    Run 'verdi process report 2431' to see them
    

    And as you can see, dos is listed in the output, so let us quickly inspect it:

    $ verdi data core.array show 2436
    ...
    

    where we have truncated the output. You can verify both the structure of the stored density of states data and its values.

    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. Of course this can form the foundations of a workflow dedicated to pre-converging results, not just for density of states calculations. We will leave this exercise to the user.

    With this we conclude the initial tutorials that follow the VASP tutorials closely and will now continue with a few tutorials related to how to interact with the plugin and finally conclude with a few rather specific tutorials that you might find interesting or inspiring.