CASTEP protocol revision (#269)

Add the verification-PBE-v1 protocol and verification-PBE-v1-a0 protocol.

*  Use a fixed cut off energy of 800 eV. This is to make sure the consistency between the oxides and the unaries set. Previously, the cut off energy is determined from the pseudopotentials on the code side automatically. However, to construct the convex hull, the same cut off energy should be used for all calculations. Previously oxides are calculated at about 750 eV cut off.
 * Fixed a bug with the soft_elements.yml file. This file provides a list of "soft" elements where the default cut off energy is very low - which can give some large pulay stress. If all of the elements are "soft" a default 326 eV cut of energy is set regardlessly. The yml files expands No and Y into False and True, which need to be properly quoted. In addition, this override should not act if cut_off_energy is set explicitly.
* The verification-PBE-v1 uses new potentials ("C19V2"), for lanthanides and actinides. The verification-PBE-v1-a0 uses instead the standard pseudo set.
* Added extractor of TS term for CASTEP

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', 'verification-PBE-v1-a0')),
+            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))
@@ -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,83 @@ fast:
                     geom_force_tol: 0.05
                     geom_energy_tol: 2.0e-5
                     geom_stress_tol: 0.1
+
+
+verification-PBE-v1-a0:
+    name: 'verification-PBE-v1-a0'
+    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"]
+
+
+verification-PBE-v1:
+    name: 'verification-PBE-v1'
+    description: 'Protocol for the oxides validation study.'
+    relax:
+        relax_options:
+            max_meta_iterations: 5
+
+        base:
+            pseudos_family: 'C19MK2'
+            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

tests/workflows/relax/test_castep.py

@@ -240,3 +240,28 @@ def test_input_generator(castep_code, nacl, si):  # pylint: disable=invalid-name
     )
     assert builder.calc.settings is None
     assert builder.base.kpoints_spacing == pytest.approx(0.023873, abs=1e-6)
+
+
+def test_otfg_upload(with_otfg):
+    """
+    Test uploading customized OTFG family
+    """
+
+    # Initial upload
+    ensure_otfg_family('C19V2')
+    assert OTFGGroup.objects.get(label='C19V2')
+
+    # Second call should not error
+    ensure_otfg_family('C19V2')
+    assert OTFGGroup.objects.get(label='C19V2')
+
+    # Second call with forced update
+    ensure_otfg_family('C19V2', force_update=True)
+
+    group = OTFGGroup.objects.get(label='C19V2')
+    found = False
+    for node in group.nodes:
+        if node.element == 'La':
+            assert node.entry == 'La 2|2.3|5|6|7|50U:60:51:52:43{4f0.1}(qc=4.5)'
+            found = True
+    assert found