# 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

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

# 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()
kpoints.set_kpoints_mesh(kmesh)
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_MAPPING = {'Si': 'Si'}

# jobfile equivalent
# In options, we typically set scheduler options.
# 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')

main(CODE_STRING, INCAR, KMESH, STRUCTURE, POTENTIAL_FAMILY, POTENTIAL_MAPPING, OPTIONS)


%/$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
label
description
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.