"""
VASP NEB workchain.
Contains the VaspNEBWorkChain class definition which uses the BaseRestartWorkChain.
"""
import numpy as np
from aiida import __version__ as aiida_version
from aiida import orm
from aiida.common.exceptions import InputValidationError, NotExistent
from aiida.common.extendeddicts import AttributeDict
from aiida.common.lang import override
from aiida.engine import while_
from aiida.engine.processes.workchains.restart import BaseRestartWorkChain, ProcessHandlerReport, process_handler
from aiida.plugins import CalculationFactory
# pylint: disable=too-many-branches, too-many-statements
from packaging import version
from aiida_vasp.assistant.parameters import ParametersMassage
from aiida_vasp.calcs.neb import VaspNEBCalculation
from aiida_vasp.data.potcar import PotcarData
from aiida_vasp.parsers.content_parsers.potcar import MultiPotcarIo
from aiida_vasp.utils.workchains import compose_exit_code
# Additional tags for VTST calculations - these are not the tags used by standard VASP
VTST_ADDITIONAL_TAGS = {
'iopt': 'TAG for VTST',
'maxmove': 'Maximum ionic movement',
'ilbfgsmem': 'Number of steps saved when building the inverse Hessian matrix',
'lglobal': 'Optimize the NEB globally instead of image-by-image',
'lautoscale': 'Automatically determines INVCURV',
'invcurv': 'Initial inverse curvature, used to construct the inverse Hessian matrix',
'llineopt': 'Use a force based line minimizer for translation',
'fdstep': 'Finite difference step size for line optimizer',
'timestep': 'Dynamical time step',
'sdalpha': 'Ratio between force and step size',
'ftimemax': 'Maximum dynamical time step allowed',
'ftimedec': 'Factor to decrease dt',
'ftimeinc': 'Factor to increase dt',
'falpha': 'Parameter that controls velocity damping',
'fnmin': 'Minium number of iterations before adjust alpha and dt',
'lclimb': 'Use climbing image mode',
'ichain': 'Indicates which method to run. NEB (ICHAIN=0) is the default',
'ltangentold': 'Flag to turn on the old central difference tangent',
'ldneb': 'Flag to turn on modified double nudging',
'lnebcell': 'Flag to turn on SS-NEB. Used with ISIF=3 and IOPT=3.',
'jacobian': 'Controls weight of lattice to atomic motion. Ω is volume and N is the number of atoms.',
}
[docs]
class VaspNEBWorkChain(BaseRestartWorkChain):
"""
A NEB workchain
Error handling enriched wrapper around VaspNEBCalculation.
Deliberately conserves most of the interface (required inputs) of the VaspNEBCalculation class, but
makes it possible for a user to interact with a workchain and not a calculation.
In addition, implement restarts of calculation when the calculation is net full converged for error handling.
"""
_verbose = False
_process_class = CalculationFactory('vasp.neb')
_norm_disp_threshold = 4.0
[docs]
@classmethod
def define(cls, spec):
super(VaspNEBWorkChain, cls).define(spec)
spec.expose_inputs(
cls._process_class,
exclude=('potential', 'kpoints', 'dynamics', 'metadata'),
)
spec.input(
'kpoints',
valid_type=orm.KpointsData,
required=False,
)
spec.input(
'kpoints_spacing',
valid_type=orm.Float,
required=False,
)
spec.input(
'potential_family',
valid_type=orm.Str,
required=True,
)
spec.input(
'potential_mapping',
valid_type=orm.Dict,
required=True,
)
spec.input(
'options',
valid_type=orm.Dict,
required=True,
)
spec.input(
'max_iterations',
valid_type=orm.Int,
required=False,
default=lambda: orm.Int(5),
help="""
The maximum number of iterations to perform.
""",
)
spec.input(
'clean_workdir',
valid_type=orm.Bool,
required=False,
default=lambda: orm.Bool(False),
help="""
If True, clean the work dir upon the completion of a successful calculation.
""",
)
spec.input(
'verbose',
valid_type=orm.Bool,
required=False,
default=lambda: orm.Bool(True),
help="""
If True, enable more detailed output during workchain execution.
""",
)
spec.input(
'dynamics.positions_dof',
valid_type=orm.List,
required=False,
help="""
Site dependent flag for selective dynamics when performing relaxation
""",
)
spec.input(
'ldau_mapping',
valid_type=orm.Dict,
required=False,
help="Mappings, see the doc string of 'get_ldau_keys'",
)
spec.input(
'kpoints_spacing',
valid_type=orm.Float,
required=False,
help='Spacing for the kpoints in units A^-1 * 2pi (CASTEP style `kpoints_mp_spacing`)',
)
spec.input(
'kpoints_spacing_vasp',
valid_type=orm.Float,
required=False,
help='Spacing for the kpoints in units A^-1 (VASP style)',
)
spec.outline(
cls.setup,
while_(cls.should_run_process)(
cls.run_process,
cls.inspect_process,
),
cls.results,
)
spec.expose_outputs(cls._process_class)
spec.exit_code(
0,
'NO_ERROR',
message='the sun is shining',
)
spec.exit_code(
700,
'ERROR_NO_POTENTIAL_FAMILY_NAME',
message='the user did not supply a potential family name',
)
spec.exit_code(
701,
'ERROR_POTENTIAL_VALUE_ERROR',
message='ValueError was returned from get_potcars_from_structure',
)
spec.exit_code(
702,
'ERROR_POTENTIAL_DO_NOT_EXIST',
message='the potential does not exist',
)
spec.exit_code(
703,
'ERROR_IN_PARAMETER_MASSAGER',
message='the exception: {exception} was thrown while massaging the parameters',
)
spec.exit_code(
501,
'SUB_NEB_CALCULATION_ERROR',
message='Unrecoverable error in launched NEB calculations.',
)
[docs]
def setup(self):
super().setup()
# Setup the initial inputs
self.ctx.inputs = self.exposed_inputs(self._process_class)
# Stage the neb images
self.ctx.neb_images = self.inputs.neb_images
# Handle and convert additional inputs and store them in self.ctx.inputs
exit_code = self._setup_vasp_inputs()
if exit_code is not None:
return exit_code
# Sanity checks
self._check_neb_inputs()
return None
[docs]
@process_handler(priority=500, exit_codes=[VaspNEBCalculation.exit_codes.ERROR_IONIC_NOT_CONVERGED]) # pylint: disable=no-member
def handle_unconverged(self, node):
"""
Handle the problem where the NEB optimization is not converged.
Note that VASP could reach NSW before the actual convergence.
Hence this check is necessary even for finished runs.
"""
if 'misc' not in node.outputs:
self.report('Cannot found the `misc` output containing the NEB run data')
return None
misc_dict = node.outputs.misc.get_dict()
neb_data = misc_dict.get('neb_data')
if neb_data is None:
self.report('Cannot found the `neb_data` dictionary containing the NEB run data')
return None
converged = [tmp.get('neb_converged', False) for tmp in neb_data.values()]
if not all(converged):
self.report('At least one image is not converged in the run. Restart required.')
# Attach images
out = self._attach_output_structure(node)
if out is not None:
return out
self.report(f'Successfully handled unconverged calculation {node}.')
return ProcessHandlerReport()
self.report(f'Cannot handle ionically unconverged calculation {node}.')
return None
[docs]
@process_handler(priority=900, exit_codes=[VaspNEBCalculation.exit_codes.ERROR_DID_NOT_FINISH]) # pylint: disable=no-member
def handle_unfinished(self, node):
"""
Handle the case where the calculations is not fully finished.
This checks the existing of the run_stats field in the parsed per-image misc output
"""
finished = []
# Since 1.6.3 the nested namespaces are handled properly.
if version.parse(aiida_version) >= version.parse('1.6.3'):
if 'misc' not in node.outputs:
self.report('Cannot found the `misc` output containing the parsed per-image data')
return None
misc_dict = node.outputs.misc.get_dict()
if 'run_status' in misc_dict:
finished = {key: value.get('finished', False) for key, value in misc_dict['run_status'].items()}
if not all(finished.values()):
self.report('At least one image did not reach the end of VASP execution - calculation not finished!')
out = self._attach_output_structure(node)
if out is not None:
return out
# No further process handling is needed
self.report(f'Successfully handled unfinished calculation {node}.')
return ProcessHandlerReport(do_break=True)
self.report(f'Cannot handle unfinished calculation {node}.')
return None
[docs]
def _attach_output_structure(self, node):
"""
Attached the output structure of a children node as the inputs for the
next workchain launch.
"""
output_images = AttributeDict() # A dictionary holding the structures with keys like 'image_xx'
output_images = node.outputs['structure']
nout = len(output_images)
nexists = len(self.inputs.neb_images)
if nout != nexists:
self.report(f'Number of parsed images: {nout} does not equal to the images need to restart: {nexists}.')
return ProcessHandlerReport(do_break=True, exit_code=self.exit_codes.SUB_NEB_CALCULATION_ERROR) # pylint: disable=no-member
self.report(f'Attached output structures from the previous calculation {node} as new inputs.')
self.ctx.inputs.neb_images = output_images
return None
[docs]
@override
def results(self):
"""Attach the outputs specified in the output specification from the last completed process."""
node = self.ctx.children[self.ctx.iteration - 1]
# We check the `is_finished` attribute of the work chain and not the successfulness of the last process
# because the error handlers in the last iteration can have qualified a "failed" process as satisfactory
# for the outcome of the work chain and so have marked it as `is_finished=True`.
max_iterations = self.inputs.max_iterations.value # type: ignore[union-attr]
if not self.ctx.is_finished and self.ctx.iteration >= max_iterations:
self.report(
f'reached the maximum number of iterations {max_iterations}: '
f'last ran {self.ctx.process_name}<{node.pk}>'
)
return self.exit_codes.ERROR_MAXIMUM_ITERATIONS_EXCEEDED # pylint: disable=no-member
self.report(f'work chain completed after {self.ctx.iteration} iterations')
# Simply attach the output of the last children
self.out_many({key: node.outputs[key] for key in node.outputs})
return None
# The code below should be moved for utility module, but I keep them here for now
FELEMS = [
'La',
'Ce',
'Pr',
'Nd',
'Pm',
'Sm',
'Eu',
'Gd',
'Tb',
'Dy',
'Ho',
'Er',
'Tm',
'Yb',
'Lu',
'Ac',
'Th',
'Pa',
'U',
'Np',
'Pu',
'Am',
'Cm',
'Bk',
'Cf',
'Es',
'Fm',
'Md',
'No',
'Lr',
]
[docs]
def get_ldau_keys(structure, mapping, utype=2, jmapping=None, felec=False):
"""
Setup LDAU mapping.
In VASP, the U for each species has to be defined in the order that they appear in POSCAR.
This is a helper to make sure the values of U are associated to each species.
:param structure: The structure, either StructureData or ase.Atoms is fine.
:param mapping: A dictionary in the format of ``{"Mn": [d, 4]...}`` for U.
:param utype: The type of LDA+U, default to 2, which is the one with only one parameter.
:param jmapping: A dictionary in the format of ``{"Mn": [d, 4]...}`` but for J.
:param felec: Whether we are dealing with f electrons; will increase lmaxmix if we are.
:returns: A dictionary to be used to update the raw input parameters for VASP.
"""
if isinstance(structure, orm.StructureData):
species = MultiPotcarIo.potentials_order(structure)
else:
# For ASE atoms, we keep the order of species occurrence no sorting is done
species = []
for symbol in structure.get_chemical_symbols():
if symbol not in species:
species.append(symbol)
lsymbols = {'d': 2, 'f': 3, 'p': 1}
if jmapping is None:
jmapping = {}
# Setup U array
ldauu = []
ldauj = []
ldaul = []
count = 0
for specie in species:
if specie in mapping:
uvalue = mapping[specie][1]
j = jmapping.get(specie, 0.0)
ldaul.append(lsymbols[mapping[specie][0]])
ldauu.append(mapping[specie][1])
j = jmapping.get(specie, 0.0)
ldauj.append(j)
if specie in FELEMS:
felec = True
# Count the number of valid mappings
if uvalue != 0.0 or j != 0.0:
count += 1
else:
ldauu.append(0.0)
ldauj.append(0.0)
ldaul.append(-1)
if count > 0:
# Only enable U is there is any non-zero value
output = {
'ldauu': ldauu,
'ldauj': ldauj,
'ldautype': utype,
'lmaxmix': 6 if felec else 4,
'ldaul': ldaul,
'ldau': True,
}
else:
output = {}
return output