diff --git a/aiida_common_workflows/workflows/bands/castep/__init__.py b/aiida_common_workflows/workflows/bands/castep/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/aiida_common_workflows/workflows/bands/castep/generator.py b/aiida_common_workflows/workflows/bands/castep/generator.py new file mode 100644 index 00000000..88a57f22 --- /dev/null +++ b/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 diff --git a/aiida_common_workflows/workflows/bands/castep/workchain.py b/aiida_common_workflows/workflows/bands/castep/workchain.py new file mode 100644 index 00000000..3f4c4f16 --- /dev/null +++ b/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 diff --git a/aiida_common_workflows/workflows/relax/castep/additional_otfg_families.yml b/aiida_common_workflows/workflows/relax/castep/additional_otfg_families.yml new file mode 100644 index 00000000..d7326774 --- /dev/null +++ b/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)" diff --git a/aiida_common_workflows/workflows/relax/castep/extractors.py b/aiida_common_workflows/workflows/relax/castep/extractors.py new file mode 100644 index 00000000..211b1f05 --- /dev/null +++ b/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 diff --git a/aiida_common_workflows/workflows/relax/castep/generator.py b/aiida_common_workflows/workflows/relax/castep/generator.py index 36c2c8e8..b7949cb4 100644 --- a/aiida_common_workflows/workflows/relax/castep/generator.py +++ b/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,8 +412,12 @@ 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 @@ -410,9 +425,29 @@ def ensure_otfg_family(family_name): # 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!") diff --git a/aiida_common_workflows/workflows/relax/castep/protocol.yml b/aiida_common_workflows/workflows/relax/castep/protocol.yml index b08a6b47..8a70749e 100644 --- a/aiida_common_workflows/workflows/relax/castep/protocol.yml +++ b/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"] diff --git a/aiida_common_workflows/workflows/relax/castep/soft_elements.yml b/aiida_common_workflows/workflows/relax/castep/soft_elements.yml index 3fcd2c1b..a97da9b5 100644 --- a/aiida_common_workflows/workflows/relax/castep/soft_elements.yml +++ b/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 diff --git a/aiida_common_workflows/workflows/relax/generator.py b/aiida_common_workflows/workflows/relax/generator.py index e84061d0..d50f7777 100644 --- a/aiida_common_workflows/workflows/relax/generator.py +++ b/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 diff --git a/tests/workflows/relax/test_castep.py b/tests/workflows/relax/test_castep.py index e27d6a1f..f3035de3 100644 --- a/tests/workflows/relax/test_castep.py +++ b/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