1. Structure relaxation¶
In this tutorial, we will first perform a calculation of the electronic
structure of wurtzite-AlN, following by a relaxation of the structure.
The idea here is to understand how we can use the builder to set up
the inputs for the call script and how to interface with the bundled
workchain for relaxing structures.
Before continuing it is worthwhile to consider symmetry and work with that in a more general manner. It is good practice to familiar ourself about these topics from time to time as much of it is easily forgotten and taken for granted, giving us problems downstream that could be avoided with a bit better oversight initially.
Symmetrizing structures¶
It is recommended to symmetrize the crystal structure of interest
if the space group type is known. This can be done by using for instance
spglib, which should already be installed when installing AiiDA-VASP.
Let’s take an example POSCAR of wurtzite-type SiC structure,
which can be obtained from the Materials project database:
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 it even more if we want more number of digits for \(\sqrt{3}\) and \(1/3\) that often appear in hexagonal crystals. An ad-hoc Python script to symmetrize this may look 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)
print('Spacegroup:', spglib.get_spacegroup(cell))
np.set_printoptions(precision=15)
[print(sym_cell[i]) for i in range(3)]
Upon saving and executing this Python script, the space group should be printed as P6_3mc (186).
And the symmetrized cell should be:
[[ 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 component of the atomic position, i.e. 0.499589, may be
expected to be zero or 0.5, but this can never be achieved exactly by spglib,
or other codes working with numerics. The take home message from this is that
you should know what kind of structure you are looking into and from time to time make checks
to make sure you are still where you think you are. Let us continue and relax this structure
using VASP and the plugin.
Relaxation of the structure¶
Here we will relax the structure, but first, we will perform an initial run without relaxing the structure so that you get familiar with the system and also are able to make the quick changes to enable relaxation from the calling script.
Initial static run¶
First assemble a script to launch a VASP calculation using the wurtzite-type SiC structure . The scripts to launch certain calculations can be designed in many different way. Let us fetch an example AiiDA-VASP run file:
$ wget https://github.com/aiida-vasp/aiida-vasp/raw/master/tutorials/run_sic.py
Inspect the file, which has the following content:
import numpy as np from aiida.common.extendeddicts import AttributeDict from aiida.engine import submit from aiida.manage.configuration import load_profile from aiida.orm import Bool, Code, Float, Int, Str from aiida.plugins import DataFactory, WorkflowFactory load_profile() def launch_aiida(structure, code_string, options, potential_family, potential_mapping, 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, 'LREAL': False, 'LCHARG': False, 'LWAVE': False, } kpoints = KpointsData() kpoints.set_kpoints_mesh([6, 6, 4], offset=[0, 0, 0.5]) 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': 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, options, potential_family, potential_mapping): structure = get_structure_SiC() launch_aiida(structure, code_string, options, potential_family, potential_mapping) if __name__ == '__main__': code_string = 'vasp@mycluster' options = { 'resources': { 'num_machines': 1, 'num_mpiprocs_per_machine': 8 }, 'account': '', 'max_memory_kb': 2000000, 'max_wallclock_seconds': 3600 } potential_family = 'PBE.54' potential_mapping = {'Si': 'Si', 'C': 'C'} main(code_string, options, potential_family, potential_mapping)
Change the relevant bits, e.g. the
code_string, theoptionsand thepotential_familyfor your stored setup. Please consult previous tutorials for details on this.Run the modified script to launch the calculation:
$ python run_sic.py
Wait a bit and once the
VaspWorkChainis finalized, get itsPKusingverdi process list -a, here2663.Check the output nodes of
2663:$ verdi process show 2663 Property Value ----------- ------------------------------------ type VaspWorkChain state Finished [0] pk 2663 uuid cc1f4f29-9c96-4e17-ba2f-3a140f789d49 label SiC VASP calculation description SiC VASP calculation ctime 2023-03-22 10:47:17.841728+01:00 mtime 2023-03-22 10:49:45.096366+01:00 Inputs PK Type ----------------- ---- ------------- clean_workdir 2660 Bool code 3 InstalledCode kpoints 2658 KpointsData max_iterations 2661 Int options 2659 Dict parameters 2653 Dict potential_family 2656 Str potential_mapping 2657 Dict settings 2655 Dict structure 2654 StructureData verbose 2662 Bool Outputs PK Type ------------- ---- ---------- energies 2668 ArrayData forces 2670 ArrayData misc 2669 Dict remote_folder 2666 RemoteData retrieved 2667 FolderData stress 2671 ArrayData Called PK Type ------------ ---- --------------- iteration_01 2665 VaspCalculation Log messages --------------------------------------------- There are 3 log messages for this calculation Run 'verdi process report 2663' to see them
And inspect for instance the
energiesoutput node:$ verdi data array show 2668 { "electronic_steps": [ 1 ], "energy_extrapolated": [ -30.09913368 ], "energy_extrapolated_electronic": [ -30.09913368 ] }Do not take these values for granted and compare them to yours. They depend on the system you executed, potential used etc.
We can also inspect it from Python. In order to make all the AiiDA machinery available in Python, we can use
verdi shell. Startverdi shelland then load the node:$ verdi shell noPython 3.10.10 (main, Mar 5 2023, 22:26:53) [GCC 12.2.1 20230201] Type 'copyright', 'credits' or 'license' for more information IPython 7.34.0 -- An enhanced Interactive Python. Type '?' for help. In [1]: n = load_node(2663) In [2]: n.outputs.energies.get_array('energy_extrapolated') Out[2]: array([-30.09913368]) In [3]: n.outputs.stress.get_array('final') Out[3]: array([[-0.01000677, 0. , 0. ], [-0. , -0.01000677, 0. ], [ 0. , 0. , 0.34873446]])Exit
verdi shellby typingexit.
Now that we have verified that the script works, let us extend it to enable relaxation of the structure.
Relaxation run¶
Let us now modify the script so that we perform a structure relaxation. If we want to fully relax the crystal structure, we need to modify the script accordingly.
Open the
run_sic.pylaunch script again.Replace
WorkflowFactory('vasp.vasp')withWorkflowFactory('vasp.relax')Remove the
IBRIONentry fromincar_dictAdd add after
builder.clean_workdir = Bool(False)the following: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)
In the end, the necessary changes to
run_sic.pycan be summarize as:--- /home/docs/checkouts/readthedocs.org/user_builds/aiida-vasp/checkouts/master/tutorials/run_sic.py +++ /home/docs/checkouts/readthedocs.org/user_builds/aiida-vasp/checkouts/master/tutorials/run_sic_relax.py @@ -15,7 +15,6 @@ incar_dict = { 'PREC': 'Accurate', - 'IBRION': -1, 'EDIFF': 1e-8, 'NELMIN': 5, 'NELM': 100, @@ -34,7 +33,7 @@ parser_settings = {'add_energies': True, 'add_forces': True, 'add_stress': True} code = Code.get_from_string(code_string) - Workflow = WorkflowFactory('vasp.vasp') + Workflow = WorkflowFactory('vasp.relax') builder = Workflow.get_builder() builder.code = code builder.parameters = Dict(dict={'incar': incar_dict}) @@ -46,6 +45,15 @@ builder.options = Dict(dict=options) builder.metadata.label = label builder.metadata.description = label + 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 cutoff + 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) builder.clean_workdir = Bool(False) node = submit(builder) return node
The plugin will then set the correct
ISIFandEDIFFGetc. The point of using dedicated settings like this is two folded: (i) it makes more sense to the user, and (ii) it makes this workflow independent on VASP and can in principle be executed with any other backend, say Quantum Espresso as long as the conversion in the backend is done properly. This is the role of theParametersMassage. The long term goal of the development of these plugins is that we will eventually have a more unified interface for the workflows that in principle can be code independent.Save the modified script and relaunch it.
Locate the
PKof the finalizedRelaxWorkChain, in this case2698and launchverdi shellagain. Then locate the relaxed structure and the stress:$ verdi shell Python 3.10.10 (main, Mar 5 2023, 22:26:53) [GCC 12.2.1 20230201] Type 'copyright', 'credits' or 'license' for more information IPython 7.34.0 -- An enhanced Interactive Python. Type '?' for help. In [1]: n = load_node(2698) In [2]: n.outputs.relax.structure.cell Out[2]: [[3.09208384, 0.0, 0.0], [-1.54604192, 2.67782316, 0.0], [0.0, 0.0, 5.07299923]] In [3]: n.outputs.stress.get_array('final') Out[3]: array([[-0.00109338, 0. , 0. ], [ 0. , -0.00109338, 0. ], [ 0. , 0. , -0.00031531]])
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, depending on what you want to use it for.