handling of SCF not converged if walltime not exceeded

Author: cpignedoli

aiida_cp2k/calculations/__init__.py

@@ -205,6 +205,7 @@ 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 thegiven threshold.")
         spec.exit_code(
             500,
             "ERROR_GEOMETRY_CONVERGENCE_NOT_REACHED",

aiida_cp2k/parsers/__init__.py

@@ -98,6 +98,8 @@ def _parse_final_structure(self):
 
     def _check_stdout_for_errors(self, output_string):
         """This function checks the CP2K output file for some basic errors."""
+        if "ABORT" in output_string and "SCF run NOT converged. To continue the calculation regardless" in output_string:
+            return self.exit_codes.ERROR_SCF_NOT_CONVERGED
 
         if "ABORT" in output_string:
             return self.exit_codes.ERROR_OUTPUT_CONTAINS_ABORT
@@ -107,6 +109,9 @@ def _check_stdout_for_errors(self, output_string):
 
         if "PROGRAM STOPPED IN" not in output_string:
             return self.exit_codes.ERROR_OUTPUT_INCOMPLETE
+        
+        if "SCF run NOT converged ***" in output_string:
+            return self.exit_codes.ERROR_SCF_NOT_CONVERGED
 
         if "MAXIMUM NUMBER OF OPTIMIZATION STEPS REACHED" in output_string:
             return self.exit_codes.ERROR_MAXIMUM_NUMBER_OPTIMIZATION_STEPS_REACHED

aiida_cp2k/utils/__init__.py

@@ -15,6 +15,7 @@
     add_ext_restart_section,
     add_first_snapshot_in_reftraj_section,
     add_wfn_restart_section,
+    add_ignore_convergence_failure,
 )
 from .parser import parse_cp2k_output, parse_cp2k_output_advanced, parse_cp2k_trajectory
 from .workchains import (
@@ -27,11 +28,13 @@
     merge_Dict,
     ot_has_small_bandgap,
     resize_unit_cell,
+    get_last_convergence_value,
 )
 
 __all__ = [
     "Cp2kInput",
     "add_ext_restart_section",
+    "add_ignore_convergence_failure",
     "add_first_snapshot_in_reftraj_section",
     "add_wfn_restart_section",
     "parse_cp2k_output",
@@ -42,6 +45,7 @@
     "check_resize_unit_cell",
     "get_input_multiplicity",
     "get_kinds_section",
+    "get_last_convergence_value",
     "merge_dict",
     "merge_Dict",
     "ot_has_small_bandgap",

aiida_cp2k/utils/input_generator.py

@@ -230,3 +230,9 @@ def add_first_snapshot_in_reftraj_section(input_dict, first_snapshot):
     params = input_dict.get_dict()
     params["MOTION"]["MD"]["REFTRAJ"]["FIRST_SNAPSHOT"] = first_snapshot
     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["FORCE_EVAL"]["DFT"]["SCF"]["IGNORE_CONVERGENCE_FAILURE"] = ".TRUE."
+    return Dict(params)

aiida_cp2k/utils/workchains.py

@@ -6,6 +6,7 @@
 ###############################################################################
 """AiiDA-CP2K utilities for workchains"""
 
+import re
 from aiida.engine import calcfunction
 from aiida.orm import Dict
 from aiida.plugins import DataFactory
@@ -99,6 +100,32 @@ def ot_has_small_bandgap(cp2k_input, cp2k_output, bandgap_thr_ev):
     is_bandgap_small = min_bandgap_ev < 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])  # RLast value found for  "DIIS/Diag"
+    
+    return None  # No value found
 
 @calcfunction
 def check_resize_unit_cell(struct, threshold):

aiida_cp2k/workchains/base.py

@@ -78,14 +78,19 @@ 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,
     ], 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'))
+        
+        # Check if time wase exceeded
+        walltime_exceeded = calc.exit_status == Cp2kCalculation.exit_codes.ERROR_OUT_OF_WALLTIME.status
 
         # CP2K was updating geometry - continue with that.
         restart_geometry_transformation = "Max. gradient              =" in 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):
             self.report("It seems that the restart of CP2K calculation wouldn't be able to fix the problem as the "
@@ -99,7 +104,12 @@ def restart_incomplete_calculation(self, calc):
         params = self.ctx.inputs.parameters
 
         params = utils.add_wfn_restart_section(params, orm.Bool('kpoints' in self.ctx.inputs))
-
+        
+        # After version 9.1 Check if calculation aborted due to SCF convergence failure and in case ignore_convergence_failure is set
+        scf_gradient = utils.get_last_convergence_value(content_string)
+        scf_restart_thr = 1e-5 # if ABORT for not SCF convergence, but RMS gradient is small, continue
+        ignore_convergence_failure = "SCF run NOT converged. To continue the calculation regardless" in content_string
+        is_geo_opt = params.get_dict().get("GLOBAL", {}).get("RUN_TYPE") == "GEO_OPT"
         if restart_geometry_transformation:
             # Check if we need to fix restart snapshot in REFTRAJ MD
             first_snapshot = None
@@ -111,8 +121,14 @@ def restart_incomplete_calculation(self, calc):
                 pass
             params = utils.add_ext_restart_section(params)
 
+        if ignore_convergence_failure and scf_gradient and scf_gradient < scf_restart_thr and not walltime_exceeded and is_geo_opt:
+            # add ignore_convergence_failure to the input parameters
+            self.report("The SCF was not converged, but the SCF gradient is small and we are doing GEO_OPT. Adding IGNORE_CONVERGENCE_FAILURE.")
+            params = utils.add_ignore_convergence_failure(params)
+
         self.ctx.inputs.parameters = params  # params (new or old ones) that include the necessary restart information.
         self.report(
             "The CP2K calculation wasn't completed. The restart of the calculation might be able to "
             "fix the problem.")
         return engine.ProcessHandlerReport(False)
+

examples/workchains/example_base_geo_opt_ignore_converge_restart.py

@@ -0,0 +1,147 @@
+###############################################################################
+# 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",
+                "WALLTIME": "00:02:30",  
+            },
+            "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 * 1 * 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(3)
+
+
+@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()