diff --git a/aiida_cp2k/calculations/__init__.py b/aiida_cp2k/calculations/__init__.py index 801f4b89..2531addb 100644 --- a/aiida_cp2k/calculations/__init__.py +++ b/aiida_cp2k/calculations/__init__.py @@ -205,6 +205,11 @@ def define(cls, spec): "ERROR_OUT_OF_WALLTIME", message="The calculation stopped prematurely because it ran out of walltime.", ) + spec.exit_code( + 450, + "ERROR_SCF_NOT_CONVERGED", + message="SCF cycle did not converge for the given threshold.", + ) spec.exit_code( 500, "ERROR_GEOMETRY_CONVERGENCE_NOT_REACHED", diff --git a/aiida_cp2k/parsers/__init__.py b/aiida_cp2k/parsers/__init__.py index 6045ef79..5ea71c98 100644 --- a/aiida_cp2k/parsers/__init__.py +++ b/aiida_cp2k/parsers/__init__.py @@ -24,6 +24,10 @@ class Cp2kBaseParser(parsers.Parser): def parse(self, **kwargs): """Receives in input a dictionary of retrieved nodes. Does all the logic here.""" + self.SEVERE_ERRORS = [ + self.exit_codes.ERROR_OUTPUT_CONTAINS_ABORT, + ] + try: _ = self.retrieved except common.NotExistent: @@ -68,13 +72,15 @@ def _parse_stdout(self): # Check the standard output for errors. exit_code = self._check_stdout_for_errors(output_string) - if exit_code: + + # Return the error code if an error was severe enough to stop the parsing. + if exit_code in self.SEVERE_ERRORS: return exit_code # Parse the standard output. result_dict = utils.parse_cp2k_output(output_string) self.out("output_parameters", orm.Dict(dict=result_dict)) - return None + return exit_code def _parse_final_structure(self): """CP2K trajectory parser.""" @@ -100,6 +106,11 @@ def _check_stdout_for_errors(self, output_string): """This function checks the CP2K output file for some basic errors.""" if "ABORT" in output_string: + if ( + "SCF run NOT converged. To continue the calculation regardless" + in output_string + ): + return self.exit_codes.ERROR_SCF_NOT_CONVERGED return self.exit_codes.ERROR_OUTPUT_CONTAINS_ABORT if "exceeded requested execution time" in output_string: @@ -219,7 +230,9 @@ def _parse_stdout(self): # Check the standard output for errors. exit_code = self._check_stdout_for_errors(output_string) - if exit_code: + + # Return the error code if an error was severe enough to stop the parsing. + if exit_code in self.SEVERE_ERRORS: return exit_code # Parse the standard output. @@ -257,7 +270,7 @@ def _parse_stdout(self): self.out("output_bands", bnds) self.out("output_parameters", orm.Dict(dict=result_dict)) - return None + return exit_code class Cp2kToolsParser(Cp2kBaseParser): @@ -275,7 +288,9 @@ def _parse_stdout(self): # Check the standard output for errors. exit_code = self._check_stdout_for_errors(output_string) - if exit_code: + + # Return the error code if an error was severe enough to stop the parsing. + if exit_code in self.SEVERE_ERRORS: return exit_code # Parse the standard output. @@ -296,4 +311,4 @@ def _parse_stdout(self): pass self.out("output_parameters", orm.Dict(dict=result_dict)) - return None + return exit_code diff --git a/aiida_cp2k/utils/__init__.py b/aiida_cp2k/utils/__init__.py index 5ddccbe0..1102222b 100644 --- a/aiida_cp2k/utils/__init__.py +++ b/aiida_cp2k/utils/__init__.py @@ -14,6 +14,7 @@ Cp2kInput, add_ext_restart_section, add_first_snapshot_in_reftraj_section, + add_ignore_convergence_failure, add_wfn_restart_section, increase_geo_opt_max_iter_by_factor, ) @@ -24,6 +25,7 @@ check_resize_unit_cell, get_input_multiplicity, get_kinds_section, + get_last_convergence_value, merge_dict, merge_Dict, ot_has_small_bandgap, @@ -33,11 +35,13 @@ __all__ = [ "Cp2kInput", "add_ext_restart_section", + "add_ignore_convergence_failure", "add_first_snapshot_in_reftraj_section", "add_wfn_restart_section", "check_resize_unit_cell", "get_input_multiplicity", "get_kinds_section", + "get_last_convergence_value", "HARTREE2EV", "HARTREE2KJMOL", "increase_geo_opt_max_iter_by_factor", diff --git a/aiida_cp2k/utils/input_generator.py b/aiida_cp2k/utils/input_generator.py index 86b34f59..99237421 100644 --- a/aiida_cp2k/utils/input_generator.py +++ b/aiida_cp2k/utils/input_generator.py @@ -240,3 +240,13 @@ def increase_geo_opt_max_iter_by_factor(input_dict, factor): factor * params.get("MOTION", {}).get("GEO_OPT", {}).get("MAX_ITER", 100) ) return Dict(params) + + +@calcfunction +def add_ignore_convergence_failure(input_dict): + """Add IGNORE_CONVERGENCE_FAILURE for non converged SCF runs.""" + params = input_dict.get_dict() + params.setdefault("FORCE_EVAL", {}).setdefault("DFT", {}).setdefault("SCF", {})[ + "IGNORE_CONVERGENCE_FAILURE" + ] = ".TRUE." + return Dict(params) diff --git a/aiida_cp2k/utils/workchains.py b/aiida_cp2k/utils/workchains.py index 095def81..0e1a9474 100644 --- a/aiida_cp2k/utils/workchains.py +++ b/aiida_cp2k/utils/workchains.py @@ -6,6 +6,8 @@ ############################################################################### """AiiDA-CP2K utilities for workchains""" +import re + from aiida.engine import calcfunction from aiida.orm import Dict from aiida.plugins import DataFactory @@ -100,6 +102,34 @@ def ot_has_small_bandgap(cp2k_input, cp2k_output, bandgap_thr_ev): return using_ot and is_bandgap_small +def get_last_convergence_value(input_string): + """Search for last "OT CG" and returns the SCF gradient. + + If no "OT CG", searches for last "DIIS/Diag" and returns the gradient. + + Args: + input_string (str): the cp2k output string. + + Returns: + float or None: the SCF gradient or None if not found. + """ + # Search all "OT CG" lines and gets the 6th column. + ot_cg_pattern = r"OT CG\s+\S+\s+\S+\s+\S+\s+\S+\s+([\d.E+-]+)" + ot_cg_matches = re.findall(ot_cg_pattern, input_string) + + if ot_cg_matches: + return float(ot_cg_matches[-1]) # Last value found for "OT CG". + + # Search for "DIIS/Diag" lines and returns the 5th column. + diis_diag_pattern = r"DIIS/Diag\.\s+\S+\s+\S+\s+\S+\s+([\d.E+-]+)" + diis_diag_matches = re.findall(diis_diag_pattern, input_string) + + if diis_diag_matches: + return float(diis_diag_matches[-1]) # Last value found for "DIIS/Diag". + + return None # No value found. + + @calcfunction def check_resize_unit_cell(struct, threshold): """Returns the multiplication factors for the cell vectors to respect, in every direction: diff --git a/aiida_cp2k/workchains/base.py b/aiida_cp2k/workchains/base.py index 76cc312a..cf5225d8 100644 --- a/aiida_cp2k/workchains/base.py +++ b/aiida_cp2k/workchains/base.py @@ -80,17 +80,33 @@ def overwrite_input_structure(self): @engine.process_handler(priority=401, exit_codes=[ Cp2kCalculation.exit_codes.ERROR_OUT_OF_WALLTIME, Cp2kCalculation.exit_codes.ERROR_OUTPUT_INCOMPLETE, + Cp2kCalculation.exit_codes.ERROR_SCF_NOT_CONVERGED, Cp2kCalculation.exit_codes.ERROR_MAXIMUM_NUMBER_OPTIMIZATION_STEPS_REACHED, ], enabled=False) def restart_incomplete_calculation(self, calc): """This handler restarts incomplete calculations.""" content_string = calc.outputs.retrieved.base.repository.get_object_content(calc.base.attributes.get('output_filename')) - # CP2K was updating geometry - continue with that. - restart_geometry_transformation = re.search(r"Max. gradient\s+=", content_string) or re.search(r"OPT\| Maximum gradient\s*[-+]?\d*\.?\d+", content_string) or "MD| Step number" in content_string - end_inner_scf_loop = "Total energy: " in content_string - # The message is written in the log file when the CP2K input parameter `LOG_PRINT_KEY` is set to True. - if not (restart_geometry_transformation or end_inner_scf_loop or "Writing RESTART" in content_string): + # CP2K was updating geometry. + possible_geometry_restart = re.search(r"Max. gradient\s+=", content_string) or re.search(r"OPT\| Maximum gradient\s*[-+]?\d*\.?\d+", content_string) or "MD| Step number" in content_string + + # CP2K wrote a wavefunction restart file. + possible_scf_restart = "Total energy: " in content_string + + # External restart file was written. + possible_ext_restart = "Writing RESTART" in content_string + + # Check if calculation aborted due to SCF convergence failure. + scf_didnt_converge_and_aborted = "SCF run NOT converged. To continue the calculation regardless" in content_string + good_scf_gradient = None + if scf_didnt_converge_and_aborted: + scf_gradient = utils.get_last_convergence_value(content_string) + scf_restart_thr = 1e-5 # if ABORT for not SCF convergence, but SCF gradient is small, continue + good_scf_gradient = (scf_gradient is not None) and (scf_gradient < scf_restart_thr) + + # Condition for allowing restart. + restart_possible = any([possible_geometry_restart, possible_scf_restart, possible_ext_restart]) and good_scf_gradient is not False + if not restart_possible: # The message is written in the log file when the CP2K input parameter `LOG_PRINT_KEY` is set to True. self.report("It seems that the restart of CP2K calculation wouldn't be able to fix the problem as the " "previous calculation didn't produce any output to restart from. " "Sending a signal to stop the Base work chain.") @@ -103,7 +119,7 @@ def restart_incomplete_calculation(self, calc): params = utils.add_wfn_restart_section(params, orm.Bool('kpoints' in self.ctx.inputs)) - if restart_geometry_transformation: + if possible_geometry_restart: # Check if we need to fix restart snapshot in REFTRAJ MD first_snapshot = None try: @@ -114,9 +130,13 @@ def restart_incomplete_calculation(self, calc): pass params = utils.add_ext_restart_section(params) + is_geo_opt = params.get_dict().get("GLOBAL", {}).get("RUN_TYPE") in ["GEO_OPT", "CELL_OPT"] + if is_geo_opt and good_scf_gradient: + self.report("The SCF was not converged, but the SCF gradient is small and we are optimising geometry. Enabling IGNORE_CONVERGENCE_FAILURE.") + params = utils.add_ignore_convergence_failure(params) + if calc.exit_code == Cp2kCalculation.exit_codes.ERROR_MAXIMUM_NUMBER_OPTIMIZATION_STEPS_REACHED: # If the maximum number of optimization steps is reached, we increase the number of steps by 40%. - print(type(params)) params = utils.increase_geo_opt_max_iter_by_factor(params, 1.4) self.ctx.inputs.parameters = params # params (new or old ones) that include the necessary restart information. diff --git a/examples/workchains/example_base_geo_opt_ignore_converge_restart.py b/examples/workchains/example_base_geo_opt_ignore_converge_restart.py new file mode 100644 index 00000000..1d07d3d7 --- /dev/null +++ b/examples/workchains/example_base_geo_opt_ignore_converge_restart.py @@ -0,0 +1,146 @@ +############################################################################### +# Copyright (c), The AiiDA-CP2K authors. # +# SPDX-License-Identifier: MIT # +# AiiDA-CP2K is hosted on GitHub at https://github.com/aiidateam/aiida-cp2k # +# For further information on the license, see the LICENSE.txt file. # +############################################################################### +"""An example testing the restart calculation handler for ENERGY run in CP2K.""" + +import os +import sys + +import ase.io +import click +from aiida.common import NotExistent +from aiida.engine import run_get_node +from aiida.orm import Dict, SinglefileData, load_code +from aiida.plugins import DataFactory, WorkflowFactory + +Cp2kBaseWorkChain = WorkflowFactory("cp2k.base") +StructureData = DataFactory("core.structure") + + +def example_base(cp2k_code): + """Run simple DFT calculation through a workchain.""" + + thisdir = os.path.dirname(os.path.realpath(__file__)) + + print("Testing CP2K ENERGY on H2O (DFT) through a workchain...") + + # Basis set. + basis_file = SinglefileData( + file=os.path.join(thisdir, "..", "files", "BASIS_MOLOPT") + ) + + # Pseudopotentials. + pseudo_file = SinglefileData( + file=os.path.join(thisdir, "..", "files", "GTH_POTENTIALS") + ) + + # Structure. + structure = StructureData( + ase=ase.io.read(os.path.join(thisdir, "..", "files", "h2o.xyz")) + ) + + # Parameters. + parameters = Dict( + { + "GLOBAL": { + "RUN_TYPE": "GEO_OPT", + }, + "FORCE_EVAL": { + "METHOD": "Quickstep", + "DFT": { + "BASIS_SET_FILE_NAME": "BASIS_MOLOPT", + "POTENTIAL_FILE_NAME": "GTH_POTENTIALS", + "QS": { + "EPS_DEFAULT": 1.0e-16, + "WF_INTERPOLATION": "ps", + "EXTRAPOLATION_ORDER": 3, + }, + "MGRID": { + "NGRIDS": 4, + "CUTOFF": 450, + "REL_CUTOFF": 70, + }, + "XC": { + "XC_FUNCTIONAL": { + "_": "LDA", + }, + }, + "POISSON": { + "PERIODIC": "none", + "PSOLVER": "MT", + }, + "SCF": { + "MAX_SCF": 10, # not enough to converge + "EPS_SCF": "1.e-6", + "PRINT": {"RESTART": {"_": "ON"}}, + }, + }, + "SUBSYS": { + "KIND": [ + { + "_": "O", + "BASIS_SET": "DZVP-MOLOPT-SR-GTH", + "POTENTIAL": "GTH-LDA-q6", + }, + { + "_": "H", + "BASIS_SET": "DZVP-MOLOPT-SR-GTH", + "POTENTIAL": "GTH-LDA-q1", + }, + ], + }, + }, + } + ) + + # Construct process builder. + builder = Cp2kBaseWorkChain.get_builder() + + # Switch on resubmit_unconverged_geometry disabled by default. + builder.handler_overrides = Dict( + {"restart_incomplete_calculation": {"enabled": True}} + ) + + # Input structure. + builder.cp2k.structure = structure + builder.cp2k.parameters = parameters + builder.cp2k.code = cp2k_code + builder.cp2k.file = { + "basis": basis_file, + "pseudo": pseudo_file, + } + builder.cp2k.metadata.options = { + "resources": { + "num_machines": 1, + "num_mpiprocs_per_machine": 1, + }, + "max_wallclock_seconds": 1 * 10 * 60, + } + + print("Submitted calculation...") + _, process_node = run_get_node(builder) + + if process_node.exit_status == 0: + print("Work chain is finished correctly.") + else: + print("ERROR! Work chain failed.") + sys.exit(1) + + +@click.command("cli") +@click.argument("codelabel") +def cli(codelabel): + """Click interface.""" + try: + code = load_code(codelabel) + except NotExistent: + print(f"The code '{codelabel}' does not exist") + sys.exit(1) + example_base(code) + + +if __name__ == "__main__": + cli()