5. Running VASP calculations explicitly¶
We will now try to build a simple call script that executes a workchain, basically a of one-shot VASP calculation for wurtzite-AlN. When that is done, we will modify this example is to also use the relax workchain, which enables relaxations of the structure.
Before running VASP calculation¶
It is recommended to well symmetrize the crystal strucutre of interest if the space group type is known. This can be done by using spglib. Spglib with python interface can be installed either pip or conda, e.g.:
% pip install spglib
or:
% conda install -c conda-forge spglib
The usage of spglib is found at https://atztogo.github.io/spglib/. Let’s take an example of POSCAR of wurtzite-type SiC structure, which can be obtained from the Materials project database, https://materialsproject.org/materials/mp-7140/
Si2 C2
1.0
3.092007 0.000000 0.000000
-1.546004 2.677757 0.000000
0.000000 0.000000 5.073347
Si C
2 2
direct
0.333333 0.666667 0.499589 Si
0.666667 0.333333 0.999589 Si
0.333333 0.666667 0.875411 C
0.666667 0.333333 0.375411 C
This is already symmetrized but we can symmetrize more if we want more number of digits for \(\sqrt{3}\) and \(1/3\) that often appear in hexagonal crystals. An ad-hoc script to symmetrize this may be like
import numpy as np
import spglib
poscar_lines = """Si2 C2
1.0
3.092007 0.000000 0.000000
-1.546004 2.677757 0.000000
0.000000 0.000000 5.073347
Si C
2 2
direct
0.333333 0.666667 0.499589 Si
0.666667 0.333333 0.999589 Si
0.333333 0.666667 0.875411 C
0.666667 0.333333 0.375411 C""".splitlines()
lattice = np.genfromtxt(poscar_lines[2:5]).reshape(3, 3)
points = np.genfromtxt(poscar_lines[8:12]).reshape(4, -1)[:, :3]
numbers = [14, 14, 6, 6]
cell = (lattice, points, numbers)
sym_cell = spglib.refine_cell(cell)
The space group type is found by
print(spglib.get_spacegroup(cell)), which should give P6_3mc
(186) in this example, and we have to be sure at least that this is
the correct one. Then, with:
np.set_printoptions(precision=15)
[print(sym_cell[i]) for i in range(3)]
we see:
[[ 3.092007293580808 0. 0. ]
[-1.546003646790404 2.677756864927749 0. ]
[ 0. 0. 5.073347 ]]
[[0.333333333333333 0.666666666666667 0.499589 ]
[0.666666666666667 0.333333333333333 0.999589 ]
[0.333333333333333 0.666666666666667 0.875411 ]
[0.666666666666667 0.333333333333333 0.375411 ]]
[14 14 6 6]
The third value of the atomic position, e.g., 0.499589, may be
expected to be zero or 0.5 by feeling, but this is not possible to be
done by spglib, because it has freedome against rigid shift along c
that doesn’t alter the symmetry.
Use of AiiDA-VASP¶
A typical script to launch a VASP caluculation (wurtzite-type SiC) is something like:
import numpy as np
from aiida.common.extendeddicts import AttributeDict
from aiida.manage.configuration import load_profile
from aiida.orm import Bool, Str, Code, Int, Float
from aiida.plugins import DataFactory, WorkflowFactory
from aiida.engine import submit
load_profile()
def launch_aiida(structure, code_string, resources,
label="SiC VASP calculation"):
Dict = DataFactory('dict')
KpointsData = DataFactory("array.kpoints")
incar_dict = {
'PREC': 'Accurate',
'IBRION': -1,
'EDIFF': 1e-8,
'NELMIN': 5,
'NELM': 100,
'ENCUT': 500,
'IALGO': 38,
'ISMEAR': 0,
'SIGMA': 0.01,
'GGA': 'PS',
'LREAL': False,
'LCHARG': False,
'LWAVE': False,
}
kpoints = KpointsData()
kpoints.set_kpoints_mesh([6, 6, 4], offset=[0, 0, 0.5])
options = {'resources': resources,
'account': '',
'max_memory_kb': 1024000,
'max_wallclock_seconds': 3600 * 10}
potential_family = 'pbe'
potential_mapping = {'Si': 'Si', 'C': 'C'}
parser_settings = {'add_energies': True,
'add_forces': True,
'add_stress': True}
code = Code.get_from_string(code_string)
Workflow = WorkflowFactory('vasp.vasp')
builder = Workflow.get_builder()
builder.code = code
builder.parameters = Dict(dict=incar_dict)
builder.structure = structure
builder.settings = Dict(dict={'parser_settings': parser_settings})
builder.potential_family = Str(potential_family)
builder.potential_mapping = Dict(dict=potential_mapping)
builder.kpoints = kpoints
builder.options = Dict(dict=options)
builder.metadata.label = label
builder.metadata.description = label
builder.clean_workdir = Bool(False)
node = submit(builder)
return node
def get_structure_SiC():
"""Set up SiC cell
Si C
1.0
3.0920072935808083 0.0000000000000000 0.0000000000000000
-1.5460036467904041 2.6777568649277486 0.0000000000000000
0.0000000000000000 0.0000000000000000 5.0733470000000001
Si C
2 2
Direct
0.3333333333333333 0.6666666666666665 0.4995889999999998
0.6666666666666667 0.3333333333333333 0.9995889999999998
0.3333333333333333 0.6666666666666665 0.8754109999999998
0.6666666666666667 0.3333333333333333 0.3754109999999997
"""
StructureData = DataFactory('structure')
a = 3.092
c = 5.073
lattice = [[a, 0, 0],
[-a / 2, a / 2 * np.sqrt(3), 0],
[0, 0, c]]
structure = StructureData(cell=lattice)
for pos_direct, symbol in zip(
([1. / 3, 2. / 3, 0],
[2. / 3, 1. / 3, 0.5],
[1. / 3, 2. / 3, 0.375822],
[2. / 3, 1. / 3, 0.875822]), ('Si', 'Si', 'C', 'C')):
pos_cartesian = np.dot(pos_direct, lattice)
structure.append_atom(position=pos_cartesian, symbols=symbol)
return structure
def main(code_string, resources):
structure = get_structure_SiC()
launch_aiida(structure, code_string, resources)
if __name__ == '__main__':
code_string = 'vasp@saga'
resources = {'num_machines': 1, 'num_mpiprocs_per_machine': 20}
main(code_string, resources)
Once the calculation is done, we locate the <pk> by using:
verdi process list
Pick the most recent VaspWorkChain process and then we can watch the results using
for instance the verdi shell:
verdi shell
And then we load the node:
In [1]: n = load_node(<pk>)
In [2]: n.outputs.energies.get_array('energy_extrapolated')
Out[2]: array([-31.80518222])
In [3]: n.outputs.stress.get_array('final')
Out[3]:
array([[-29.89502712, 0. , 0. ],
[ 0. , -29.89502712, 0. ],
[ 0. , 0. , -29.47075517]])
When we want to fully relax a crystal structure, the above script is modified as follows:
Replace
WorkflowFactory('vasp.vasp')byWorkflowFactory('vasp.relax')Remove
IBRIONfromincar_dictAdd the following setting:
relax = AttributeDict() relax.perform = Bool(True) # Turn on relaxation of the structure relax.force_cutoff = Float(1e-5) # Relax force cutoff relax.steps = Int(10) # Relax number of ionic steps relax.positions = Bool(True) # Relax atomic positions relax.shape = Bool(True) # Relax cell shape (alpha, beta, gamma) relax.volume = Bool(True) # Relax volume builder.relax = relax builder.verbose = Bool(True)
The lattice parameters of the relax crystal structure is found by
In [1]: n = load_node(<PK>)
In [2]: n.outputs.relax__structure.cell
Out[2]:
[[3.07798535, 0.0, 0.0],
[-1.53899268, 2.66561351, 0.0],
[0.0, 0.0, 5.04931673]]
In [3]: n.outputs.stress.get_array('final')
Out[3]:
array([[-0.01708304, 0. , 0. ],
[ 0. , -0.01708304, 0. ],
[ 0. , 0. , -0.00809151]])
There are more options for the relax workchain, e.g., running VASP several time iteratively until convergence, which is used in the bulk modulus example in the next section.
After the relaxation, sometimes the crystal symmetry can be slightly broken by the VASP calculation, especially for hexagonal crystals. So it is recommended to symmetrize the final structure if this is the case.