Implement CommonBandsWorkChain for CASTEP

aiida_common_workflows/workflows/bands/castep/generator.py

@@ -0,0 +1,81 @@
+# -*- coding: utf-8 -*-
+"""Implementation of `aiida_common_workflows.common.bands.generator.CommonBandsInputGenerator` for CASTEP."""
+from aiida import engine, orm
+
+from aiida_common_workflows.generators import CodeType
+
+from ..generator import CommonBandsInputGenerator
+
+__all__ = ('CastepCommonBandsInputGenerator',)
+
+
+class CastepCommonBandsInputGenerator(CommonBandsInputGenerator):
+    """Generator of inputs for the CastepCommonBandsWorkChain"""
+
+    @classmethod
+    def define(cls, spec):
+        """Define the specification of the input generator.
+
+        The ports defined on the specification are the inputs that will be accepted by the ``get_builder`` method.
+        """
+        super().define(spec)
+        spec.inputs['engines']['bands']['code'].valid_type = CodeType('castep.castep')
+
+    def _construct_builder(self, **kwargs) -> engine.ProcessBuilder:
+        """Construct a process builder based on the provided keyword arguments.
+
+        The keyword arguments will have been validated against the input generator specification.
+        """
+        # pylint: disable=too-many-branches,too-many-statements,too-many-locals
+
+        required_keys = ('engines', 'parent_folder', 'bands_kpoints')
+        for key in required_keys:
+            if key not in kwargs:
+                raise ValueError(f'Required key `{key}` is missing in the function argument.')
+        engines = kwargs['engines']
+        parent_folder = kwargs['parent_folder']
+        bands_kpoints = kwargs['bands_kpoints']
+
+        # From the parent folder, we retrieve the calculation that created it. Note
+        # that we are sure it exists (it wouldn't be the same for WorkChains). We then check
+        # that it is a CastepCalculation and create the builder.
+        parent_castep_calc = parent_folder.creator
+        if parent_castep_calc.process_type != 'aiida.calculations:castep.castep':
+            raise ValueError('The `parent_folder` has not been created by a CastepCalculation')
+        builder_castep_calc = parent_castep_calc.get_builder_restart()
+
+        # Construct the builder of the `common_bands_wc` from the builder of a CastepCalculation.
+        builder_common_bands_wc = self.process_class.get_builder()
+        builder_common_bands_wc.scf.calc_options = orm.Dict(dict=dict(builder_castep_calc.metadata.options))
+        # Ensure we use castep_bin for restart, instead of the check file
+        #builder_common_bands_wc.scf.options = orm.Dict(dict={'use_castep_bin': True})
+
+        builder_castep_calc.metadata = {}
+
+        # Attach inputs of the calculation
+        for key, value in builder_castep_calc.items():
+            if value and key not in ['metadata', 'structure']:
+                builder_common_bands_wc.scf.calc[key] = value
+
+        # Updated the structure (in case we have one in output)
+        if 'output_structure' in parent_castep_calc.outputs:
+            builder_common_bands_wc.structure = parent_castep_calc.outputs.output_structure
+        else:
+            builder_common_bands_wc.structure = parent_castep_calc.inputs.structure
+
+        try:
+            engb = engines['bands']
+        except KeyError:
+            raise ValueError('The engines dictionary passed must contains a key named `bands`.')
+
+        builder_common_bands_wc.scf.calc.code = engb['code']
+
+        if 'options' in engb:
+            builder_common_bands_wc.scf.calc.metadata.options = engb['options']
+
+        # Set the `bandskpoints` and the `parent_calc_folder` for restart
+        builder_common_bands_wc.bands_kpoints = bands_kpoints
+        builder_common_bands_wc.scf.continuation_folder = parent_folder
+        builder_common_bands_wc.run_separate_scf = orm.Bool(False)
+
+        return builder_common_bands_wc

aiida_common_workflows/workflows/bands/castep/workchain.py

@@ -0,0 +1,41 @@
+# -*- coding: utf-8 -*-
+"""Implementation of `aiida_common_workflows.common.relax.workchain.CommonRelaxWorkChain` for CASTEP."""
+from logging import getLogger
+
+from aiida.engine import calcfunction
+from aiida.orm import Float
+from aiida.plugins import WorkflowFactory
+
+from ..workchain import CommonBandsWorkChain
+from .generator import CastepCommonBandsInputGenerator
+
+__all__ = ('CastepCommonBandsWorkChain',)
+
+
+@calcfunction
+def get_fermi_energy(bands):
+    """Extract the Fermi energy from the BandsData output"""
+    efermi = bands.get_attribute('efermi')
+    if isinstance(efermi, list):
+        efermi = efermi[0]
+        logger = getLogger(__name__)
+        logger.warning('Spin polarised calculation - using the efermi energy of the first spin channel.')
+    return Float(efermi)
+
+
+class CastepCommonBandsWorkChain(CommonBandsWorkChain):
+    """Implementation of `aiida_common_workflows.common.bands.workchain.CommonBandsWorkChain` for CASTEP."""
+
+    _process_class = WorkflowFactory('castep.bands')
+    _generator_class = CastepCommonBandsInputGenerator
+
+    def convert_outputs(self):
+        """Convert the outputs of the sub workchain to the common output specification."""
+        if 'band_structure' not in self.ctx.workchain.outputs:
+            self.report('CastepBandsWorkChain concluded without returning bands!')
+            return self.exit_codes.ERROR_SUB_PROCESS_FAILED
+
+        self.out('fermi_energy', get_fermi_energy(self.ctx.workchain.outputs['band_structure']))
+
+        self.out('bands', self.ctx.workchain.outputs['band_structure'])
+        return None

aiida_common_workflows/workflows/relax/castep/additional_otfg_families.yml

@@ -0,0 +1,32 @@
+# Family based on C19 with updated potentials for lanthanide and actinides
+C19V2:
+  - "C19"
+  - "La 2|2.3|5|6|7|50U:60:51:52:43{4f0.1}(qc=4.5)"
+  - "Ce 2|2.2|8|9|10|50U:60:51:52:43{5d0.1}(qc=4.5)"
+  - "Pr 2|2.1|10|12|13|50U:60:51:52:43{5d0.1}(qc=5)"
+  - "Nd 2|2.1|10|12|13|50U:60:51:52:43{5d0.1}(qc=5)"
+  - "Pm 2|2.1|8|9|11|50U:60:51:52:43{5d0.1,4f4}(qc=5.5)"
+  - "Sm 2|2.1|9|10|12|50U:60:51:52:43{5d0.1,4f5}(qc=5.5)"
+  - "Eu 2|2.1|9|10|12|50U:60:51:52:43{5d0.1,4f6}(qc=5.5)"
+  - "Gd 3|2.1|9|10|12|50U:60:51:52:43(qc=5.5)"
+  - "Tb 2|2.2|12|13|15|50U:60:51:52:43{5d0.1}(qc=5)"
+  - "Dy 2|2.0|12|13|15|50U:60:51:52:43{5d0.1}(qc=6.5)"
+  - "Ho 2|2.0|12|13|15|50U:60:51:52:43{5d0.1}(qc=6.5)"
+  - "Er 2|2.1|10|12|13|50U:60:51:52:43{6s0.1,5d0.1}(qc=6)"
+  - "Tm 2|2.1|10|12|13|50U:60:51:52:43{5d0.1,4f12}(qc=6)"
+  - "Yb 2|2.1|10|12|13|50U:60:51:52:43{5d0.1,4f13}(qc=6)"
+  - "Ac 2|2.4|7|7|9|60U:70U2U2:61:62:53{6d0.1,5f0.1}(qc=5)"
+  - "Th 2|2.2|7|7|9|60U:70U2U2:61:62:53{5f0.1}(qc=5)"
+  - "Pa 2|2.2|8|9|10|60U:70U2U2:61:62:53(qc=5)"
+  - "U  2|2.2|8|9|10|60U:70U2U2:61:62:53(qc=5)"
+  - "Np 2|2.2|9|10|12|60U:70U2U2:61:62:53(qc=5)"
+  - "Pu 2|2.2|9|10|12|60U:70U2U2:61:62:53{6d0.1}(qc=5.5)"
+  - "Am 2|2.2|9|10|12|60U:70U2U2:61:62:53{6d0.1}(qc=5.5)"
+  - "Cm 2|2.2|9|10|12|60U:70U2U2:61:62:53(qc=5.5)"
+  # Elements below in this family are not tested again all-electron codes - use with CAUTION
+  - "Bk 2|2.2|9|10|12|60U:70U2U2:61:62:53{6d0.1}(qc=5.5)"
+  - "Cf 2|2.2|9|10|12|60U:70U2U2:61:62:53{6d0.1}(qc=5.5)"
+  - "Es 2|2.1|10|12|13|60U:70U2U2:61:62:53{6d0.1}(qc=6)"
+  - "Fm 2|2.1|10|12|13|60U:70U2U2:61:62:53{6d0.1}(qc=6)"
+  - "Md 2|2.0|10|12|13|60U:70U2U2:61:62:53{6d0.1,5f12}(qc=6)"
+  - "No 2|2.0|10|12|13|60U:70U2U2:61:62:53{6d0.1,5f13}(qc=6)"

aiida_common_workflows/workflows/relax/castep/extractors.py

@@ -0,0 +1,37 @@
+# -*- coding: utf-8 -*-
+"""
+Collects some functions to postprocess a `CastepCommonRelaxWorkChain`.
+"""
+from aiida.common import LinkType
+from aiida.orm import WorkChainNode
+from aiida.plugins import WorkflowFactory
+
+CastepCommonRelaxWorkChain = WorkflowFactory('common_workflows.relax.castep')
+
+
+def get_ts_energy(common_relax_workchain):
+    """
+    Return the TS value of a concluded CastepCommonRelaxWorkChain.
+
+    CASTEP reports three quantities related to the energy:
+
+      - "free energy": ``Final free energy`` -  the total energy minus the TS term.
+        This is the energy that gets minimised and is consistent with the forces calculated.
+      - "total energy": ``Final energy`` - the total Khon-Sham energy.
+      - "extrapolated 0K energy": ``NB est. 0K energy`` - the result of E-0.5TS to give better convergence.
+      - "enthalpy": Is the free energy minus to PV term under finite temperature (for geometry optimisation).
+
+    The TS term can extrapolated by subtracting the free energy from the total energy.
+    """
+    if not isinstance(common_relax_workchain, WorkChainNode):
+        return ValueError('The input is not a workchain (instance of `WorkChainNode`)')
+    if common_relax_workchain.process_class != CastepCommonRelaxWorkChain:
+        return ValueError('The input workchain is not a `CastepCommonRelaxWorkChain`')
+
+    castep_base_wc = common_relax_workchain.get_outgoing(link_type=LinkType.CALL_WORK).one().node
+    e_ks = castep_base_wc.outputs.output_parameters['total energy']
+    free_e = castep_base_wc.outputs.output_parameters['free energy']
+
+    ts = e_ks - free_e  #pylint: disable=invalid-name
+
+    return ts

aiida_common_workflows/workflows/relax/castep/generator.py

@@ -18,6 +18,7 @@
 from ..generator import CommonRelaxInputGenerator
 
 # pylint: disable=import-outside-toplevel, too-many-branches, too-many-statements
+KNOWN_BUILTIN_FAMILIES = ('C19', 'NCP19', 'QC5', 'C17', 'C9')
 
 __all__ = ('CastepCommonRelaxInputGenerator',)
 
@@ -46,6 +47,13 @@ def define(cls, spec):
         The ports defined on the specification are the inputs that will be accepted by the ``get_builder`` method.
         """
         super().define(spec)
+        spec.input(
+            'protocol',
+            valid_type=ChoiceType(('fast', 'moderate', 'precise', 'verification-PBE-v1')),
+            default='moderate',
+            help='The protocol to use for the automated input generation. This value indicates the level of precision '
+            'of the results and computational cost that the input parameters will be selected for.',
+        )
         spec.inputs['spin_type'].valid_type = ChoiceType((SpinType.NONE, SpinType.COLLINEAR, SpinType.NON_COLLINEAR))
         spec.inputs['relax_type'].valid_type = ChoiceType(tuple(RelaxType))
         spec.inputs['electronic_type'].valid_type = ChoiceType((ElectronicType.METAL, ElectronicType.INSULATOR))
@@ -73,7 +81,7 @@ def _construct_builder(self, **kwargs) -> engine.ProcessBuilder:
 
         # Because the subsequent generators may modify this dictionary and convert things
         # to AiiDA types, here we make a full copy of the original protocol
-        protocol = copy.deepcopy(self.get_protocol(protocol))
+        protocol = copy.deepcopy(self._protocols[protocol])
         code = engines['relax']['code']
 
         override = {'base': {'calc': {'metadata': {'options': engines['relax']['options']}}}}
@@ -161,11 +169,13 @@ def _construct_builder(self, **kwargs) -> engine.ProcessBuilder:
         # Raise the cut off energy for very soft pseudopotentials
         # this is because the small basis set will give rise to errors in EOS / variable volume
         # relaxation even with the "fine" option
-        with open(str(pathlib.Path(__file__).parent / 'soft_elements.yml')) as fhandle:
-            soft_elements = yaml.safe_load(fhandle)
-        symbols = [kind.symbol for kind in structure.kinds]
-        if all(sym in soft_elements for sym in symbols):
-            param['cut_off_energy'] = 326  # eV, approximately 12 Ha
+        if 'cut_off_energy' not in protocol['relax']['base']['calc']['parameters']:
+            with open(str(pathlib.Path(__file__).parent / 'soft_elements.yml')) as fhandle:
+                soft_elements = yaml.safe_load(fhandle)
+            symbols = [kind.symbol for kind in structure.kinds]
+            if all(sym in soft_elements for sym in symbols):
+                param['cut_off_energy'] = 326  # eV, approximately 12 Ha
+                param.pop('basis_precision', None)
 
         # Apply the overrides
         if param:
@@ -342,7 +352,8 @@ def generate_inputs_base(
         'kpoints_spacing': orm.Float(merged['kpoints_spacing'] / 2 / pi),
         'max_iterations': orm.Int(merged['max_iterations']),
         'pseudos_family': orm.Str(otfg_family.label),
-        'calc': calc_dictionary
+        'calc': calc_dictionary,
+        'ensure_gamma_centering': orm.Bool(merged.get('ensure_gamma_centering', False)),
     }
 
     return dictionary
@@ -401,18 +412,42 @@ def generate_inputs_calculation(
     return dictionary
 
 
-def ensure_otfg_family(family_name):
-    """Add common OTFG families if they do not exist"""
+def ensure_otfg_family(family_name, force_update=False):
+    """
+    Add common OTFG families if they do not exist
+    NOTE: CASTEP also supports UPF families, but it is not enabled here, since no UPS based protocol
+    has been implemented.
+    """
 
     from aiida.common import NotExistent
     from aiida_castep.data.otfg import upload_otfg_family
 
     # Ensure family name is a str
     if isinstance(family_name, orm.Str):
         family_name = family_name.value
-
     try:
         OTFGGroup.objects.get(label=family_name)
     except NotExistent:
-        description = f"CASTEP built-in on-the-fly generated pseudos libraray '{family_name}'"
-        upload_otfg_family([family_name], family_name, description, stop_if_existing=True)
+        has_family = False
+    else:
+        has_family = True
+
+    # Check if it is builtin family
+    if family_name in KNOWN_BUILTIN_FAMILIES:
+        if not has_family:
+            description = f"CASTEP built-in on-the-fly generated pseudos libraray '{family_name}'"
+            upload_otfg_family([family_name], family_name, description, stop_if_existing=True)
+        return
+
+    # Not an known family - check if it in the additional settings list
+    # Load configuration from the settings
+    with open(str(pathlib.Path(__file__).parent / 'additional_otfg_families.yml')) as handle:
+        additional = yaml.safe_load(handle)
+
+    if family_name in additional:
+        if not has_family or force_update:
+            description = f"Modified CASTEP built-in on-the-fly generated pseudos libraray '{family_name}'"
+            upload_otfg_family(additional[family_name], family_name, description, stop_if_existing=False)
+    elif not has_family:
+        # No family found - and it is not recognized
+        raise RuntimeError(f"Family name '{family_name}' is not recognized!")

aiida_common_workflows/workflows/relax/castep/protocol.yml

@@ -83,3 +83,43 @@ fast:
                     geom_force_tol: 0.05
                     geom_energy_tol: 2.0e-5
                     geom_stress_tol: 0.1
+
+
+verification-PBE-v1:
+    name: 'verification-PBE-v1'
+    description: 'Protocol for the oxides validation study.'
+    relax:
+        relax_options:
+            max_meta_iterations: 5
+
+        base:
+            pseudos_family: 'C19V2'
+            max_iterations: 5
+            kpoints_spacing: 0.06    # K point spacing in A^-1, not the A^-1 2Pi that CASTEP defaults to
+            ensure_gamma_centering: True   # Ensure that the kpoint grid is Gamma-centered
+            calc:
+                parameters:
+                    task: geometryoptimisation
+                    xc_functional: pbe
+                    symmetry_generate: true
+                    snap_to_symmetry: true
+                    calculate_stress: true
+                    cut_off_energy: 800              # Fixed cut off energy for handling both unaries and oxides
+                    write_otfg: false
+                    opt_strategy: speed
+                    write_bib: false
+                    write_checkpoint: minimal
+                    finite_basis_corr: none   # No finite basis corr as we are NOT moving any atoms....
+                    max_scf_cycles: 100
+                    elec_energy_tol: 1e-8
+                    geom_force_tol: 0.03
+                    geom_energy_tol: 1.0e-5
+                    geom_stress_tol: 0.05
+                    grid_scale: 2
+                    fine_grid_scale : 3
+                    # Revised smearing scheme for the oxide validation project
+                    smearing_width: 0.06122561905370023   # Equivalent to 0.0045 Ry
+                    smearing_scheme: FERMIDIRAC
+                    perc_extra_bands: 80
+                settings:
+                    ADDITIONAL_RETRIEVE_TEMPORARY_LIST: ["aiida.castep_bin"]

aiida_common_workflows/workflows/relax/castep/soft_elements.yml

@@ -16,7 +16,7 @@
 - Kr
 - Rb
 - Sr
-- Y
+- 'Y'
 - Zr
 - In
 - Sn
@@ -49,7 +49,7 @@
 - Ac
 - Th
 - Pa
-- No
+- 'No'
 - Lr
 - Rf
 - Db

aiida_common_workflows/workflows/relax/generator.py

@@ -6,12 +6,11 @@
 
 from aiida_common_workflows.common import ElectronicType, RelaxType, SpinType
 from aiida_common_workflows.generators import ChoiceType, InputGenerator
-from aiida_common_workflows.protocol import ProtocolRegistry
 
 __all__ = ('CommonRelaxInputGenerator',)
 
 
-class CommonRelaxInputGenerator(InputGenerator, ProtocolRegistry, metaclass=abc.ABCMeta):
+class CommonRelaxInputGenerator(InputGenerator, metaclass=abc.ABCMeta):
     """Input generator for the common relax workflow.
 
     This class should be subclassed by implementations for specific quantum engines. After calling the super, they can