Geometry optimisation#
Note
This notebook can be downloaded as silicon_relax.ipynb and silicon_relax.md
In this example we will perform a geometry optimisation of silicon using VASP. We will use the VASP code and the AiiDA plugin for VASP.
Setting up the environment#
The code block below configures the environment for this example. Please see the single point example for more details.
from aiida_vasp.utils.temp_profile import *
print(load_temp_profile())
# Uncomment the below line to create a localhost Computer if you have not done so
comp = orm.Computer('localhost', 'localhost', transport_type='core.local', scheduler_type='core.direct')
comp.store()
# Some configuration may be needed for first-time user
from pathlib import Path
import os
comp.set_workdir('/tmp/aiida_run/')
comp.configure()
vasp_path = !which mock-vasp
vasp_code = orm.InstalledCode(comp, vasp_path[0], default_calc_job_plugin='vasp.vasp')
vasp_code.label ='mock-vasp'
vasp_code.store()
os.environ['MOCK_VASP_REG_BASE'] = str((Path() / 'mock_registry').absolute())
os.environ['MOCK_VASP_UPLOAD_PREFIX'] = 'relax'
# Upload the POTCAR files
from aiida_vasp.data.potcar import PotcarData, PotcarFileData
from pathlib import Path
print(PotcarData.upload_potcar_family(str(Path('potcars').absolute()), "PBE.EXAMPLE", "PBE.EXAMPLE"))
# Setting up the silicon structure
from ase.build import bulk
si = bulk('Si', 'diamond', 5.4)
si_node = orm.StructureData(ase=si)
Profile<uuid='459c3e84deae41a093ae458d4fa40df3' name='myprofile'>
(1, 1, 1)
Running the relaxation#
Similar to the single point calculation tutorial, we will use a VaspInputGenerator to setup the inputs
for the VaspRelaxWorkChain.
from aiida import orm
from aiida_vasp.protocols.generator import VaspRelaxInputGenerator
upd = VaspRelaxInputGenerator()
# Override the vasp.potential_family input of the builder
upd.get_builder(structure=si_node, code='mock-vasp@localhost', overrides={
'vasp': {'potential_family': 'PBE.EXAMPLE'}}
)
upd.builder
Process class: VaspRelaxWorkChain
Inputs:
relax_settings:
algo: cg
clean_reuse: true
convergence_absolute: false
convergence_max_iterations: 5
convergence_mode: last
convergence_on: true
convergence_positions: 0.1
convergence_shape_angles: 0.1
convergence_shape_lengths: 0.1
convergence_volume: 0.01
double_relax_mode: false
energy_cutoff: null
force_cutoff: 0.03
hybrid_calc_bootstrap: false
hybrid_calc_bootstrap_wallclock: 3600
keep_magnetization: false
keep_sp_workdir: false
perform: true
positions: true
residual_forces_check: true
reuse: false
shape: true
steps: 60
volume: true
static_calc_options: {}
static_calc_parameters: {}
static_calc_settings: {}
structure: Si
vasp:
calc:
metadata:
options:
max_wallclock_seconds: 3600
resources:
num_machines: 1
tot_num_mpiprocs: 1
withmpi: false
clean_workdir: true
code: mock-vasp@localhost
kpoints_spacing: 0.05
magmom_mapping:
Ce: 5
Co: 5
Cr: 5
Dy: 5
Eu: 10
Fe: 5
Gd: 7
Ho: 4
La: 0.6
Lu: 0.6
Mn: 5
Mo: 5
Nd: 3
Ni: 5
Pm: 4
Pr: 2
Sm: 5
Tb: 6
Tm: 2
V: 5
W: 5
Yb: 1
default: 1.0
max_iterations: 5
parameters:
incar:
algo: fast
ediff: 2.0e-06
encut: 250
gga: ps
ismear: 0
ispin: 2
lasph: true
lcharg: false
lorbit: null
lreal: false
lvhar: true
lwave: false
nedos: 100
nelm: 200
nelmin: 4
nwrite: 1
prec: normal
sigma: 0.05
potential_family: PBE.EXAMPLE
potential_mapping:
Si: Si
settings: {}
verbose: false
Here we can see that the input nodes are different from that of the single point calculation (workflow).
Nodes such as the settings and parameters go into a vasp input names rather than at the root level.
This is because the VaspRelaxWorkChain exposes the inputs of the VaspWorkChain in the vasp namespace.
Also note that the structure input is still at the root level.
This is because the atomic structure is an essential input for the relaxation (and another calculations).
Note that there is a also a relax_settings input, which contains control parameters for the relaxation.
You can see the available options by inspecting the field of the builder:
upd.get_input_help('relax_settings')
algo: str
Default: The algorithm to use for relaxation
energy_cutoff: Optional
Default: The cut off energy difference when the relaxation is stopped (e.g. EDIFF)
force_cutoff: float
Default: The maximum force when the relaxation is stopped (e.g. EDIFFG)
steps: int
Default: Number of relaxation steps to perform (eg. NSW)
positions: bool
Default: If True, perform relaxation of the atomic positions
shape: bool
Default: If True, perform relaxation of the cell shape
volume: bool
Default: If True, perform relaxation of the cell volume
convergence_on: bool
Default: If True, perform convergence checks within the workchain
convergence_absolute: bool
Default: If True, use absolute values where possible when performing convergence checks
convergence_max_iterations: int
Default: Maximum iterations for convergence checking
convergence_positions: float
Default: The cutoff value for the convergence check on positions in Angstram. A negative value by pass the check.
convergence_volume: float
Default: The cutoff value for the convergence check on volume between the two structures. A negative value by pass the check.
convergence_shape_lengths: float
Default: The cutoff value for the convergence check on the lengths of the unit cell vectors, between input and the outputs structure. A negative value by pass the check.
convergence_shape_angles: float
Default: The cutoff value for the convergence check on the angles of the unit cell vectors, between input and the outputs structure. A negative value by pass the check.
convergence_mode: str
Default: Mode of the convergence check for positions. 'inout' for checking input/output structure, or 'last' to check only the change of the last step.
reuse: bool
Default: Whether reuse the previous calculation by copying over the remote folder
clean_reuse: bool
Default: Whether to perform a final cleaning of the reused calculations
keep_sp_workdir: bool
Default: Whether to keep the workdir of the final singlepoint calculation
perform: bool
Default: Do not perform any relaxation if set to 'False'
hybrid_calc_bootstrap: bool
Default: Whether to bootstrap hybrid calculation by performing standard DFT first
hybrid_calc_bootstrap_wallclock: int
Default: Wall time limit in second for the bootstrap calculation
keep_magnetization: bool
Default: Whether to keep magnetization from the previous calculation if possible
double_relax_mode: bool
Default: Experimental: Run in double relax mode - launch of the sub workflow is only performed up to two times without checking convergence in the end. This is useful for cases where the convergence is difficult due to change of basis set with variable cell and high-throughput studies.
residual_forces_check: bool
Default: Whether to perform residual force check after relaxation to ensure the forces are below threshold
or simply hit <Shift>+<Tab> after upd.builder.relax_settings in the Notebook code cell.
We can now run the relaxation using the same run_get_node method as in the single point example.
results = upd.run_get_node()
The VaspRelaxWorkChain also has a misc output node, which is in fact the output of the last VaspWorkChain it launched.
results.node.outputs.misc.get_dict()
{'total_energies': {'energy_extrapolated': -10.62045741,
'energy_extrapolated_electronic': -10.62045741},
'stress': [[0.07859658, 0.0, -0.0],
[-0.0, 0.07859658, 0.0],
[0.0, 0.0, 0.07859658]],
'forces': [[0.0, 0.0, -0.0], [-0.0, -0.0, 0.0]],
'fermi_level': 6.03669271,
'band_properties': {'cbm': 6.2549,
'vbm': 5.7938,
'is_direct_gap': False,
'band_gap': 0.4611},
'version': '6.2.0',
'magnetization': [-0.0004373],
'site_magnetization': {'sphere': {'x': {'site_moment': {},
'total_magnetization': {}},
'y': {'site_moment': {}, 'total_magnetization': {}},
'z': {'site_moment': {}, 'total_magnetization': {}}},
'full_cell': [-0.0004373]},
'run_stats': {'mem_usage_base': 30000.0,
'mem_usage_nonl-proj': 478.0,
'mem_usage_fftplans': 1184.0,
'mem_usage_grid': 2809.0,
'mem_usage_one-center': 12.0,
'mem_usage_wavefun': 2194.0,
'total_cpu_time_used': 3.241,
'user_time': 3.221,
'system_time': 0.019,
'elapsed_time': 3.566,
'maximum_memory_used': 100124.0,
'average_memory_used': None},
'run_status': {'nelm': 200,
'nsw': 0,
'last_iteration_index': [1, 13],
'finished': True,
'ionic_converged': None,
'electronic_converged': True,
'consistent_nelm_breach': False,
'contains_nelm_breach': False,
'nbands': 9},
'notifications': [{'name': 'xc_enforced',
'kind': 'ADVICE',
'message': 'You enforced a specific xc type in the INCAR file but a different\ntype was found in the POTCAR file.\nI HOPE YOU KNOW WHAT YOU ARE DOING!',
'regex': 'You enforced a specific xc type in the INCAR file but a different'}]}
The final relaxed structure is stored in the relax.structure output node.
relaxed_node = results.node.outputs.relax.structure
print(f"Volume before relaxations: {si_node.get_cell_volume():3f} A^3")
print(f"Volume after relaxations: {relaxed_node.get_cell_volume():3f} A^3")
Volume before relaxations: 39.366000 A^3
Volume after relaxations: 40.087704 A^3
Why use VaspRelaxWorkChain?#
Surely VASP has the option to run relaxation itself (ISIF=3, IBRION=2, NSW=60), so why do we need a VaspRelaxWorkChain in the first place?
The reason is tha VaspRelaxWorkChain run more strict check and verifies if the structure is really fully relaxed.
Also, since VASP run variable cell relaxation in the fixed-basis mode, the effective plane wave cut off energies changes with the cell volume inside a single calculation.
This means multiple calculations are needed to get fully (and self-consistently) relaxed structures.
In addition, the VaspRelaxWorkChain will run a single point calculation in the end to get the correct energy that is consistent with the plane wave cut-off energy (ENCUT) originally set.
The called and called_descendants property can be used to inspect the sub-processes launched by the VaspRelaxWorkChain
print("Directly called workchains:", results.node.called)
print("All called processes:", results.node.called_descendants)
Directly called workchains: [<WorkChainNode: uuid: c0cafddf-26ff-4d1b-a85f-f9e5f02e6d07 (pk: 23) (aiida.workflows:vasp.v2.vasp)>, <WorkChainNode: uuid: 17140be9-1c65-4676-b73c-dde7f138d668 (pk: 35) (aiida.workflows:vasp.v2.vasp)>, <WorkChainNode: uuid: d613cf15-20c8-474d-985c-91dfa4ecde7f (pk: 46) (aiida.workflows:vasp.v2.vasp)>]
All called processes: [<WorkChainNode: uuid: c0cafddf-26ff-4d1b-a85f-f9e5f02e6d07 (pk: 23) (aiida.workflows:vasp.v2.vasp)>, <CalcJobNode: uuid: 607845b7-09cd-42b3-890a-3ee66c5e65f6 (pk: 28) (aiida.calculations:vasp.vasp)>, <WorkChainNode: uuid: 17140be9-1c65-4676-b73c-dde7f138d668 (pk: 35) (aiida.workflows:vasp.v2.vasp)>, <CalcJobNode: uuid: c6c3639b-711f-46bc-ab76-88db70929226 (pk: 40) (aiida.calculations:vasp.vasp)>, <WorkChainNode: uuid: d613cf15-20c8-474d-985c-91dfa4ecde7f (pk: 46) (aiida.workflows:vasp.v2.vasp)>, <CalcJobNode: uuid: d075c0d4-060c-41ee-97e1-5e12ac2be54a (pk: 51) (aiida.calculations:vasp.vasp)>]