Merge pull request #1 from molssi-seamm/dev

Dev

.github/workflows/BranchCI.yaml

@@ -10,4 +10,4 @@ jobs:
     name: Branch CI
     uses: molssi-seamm/devops/.github/workflows/BranchCI.yaml@main
     with:
-      src : seamm_ase
+      src : seamm_geometric

HISTORY.rst

@@ -1,6 +1,11 @@
 =======
 History
 =======
+2024.10.15 -- Added various convergence metrics to results
+    * Added maximum force, gradient, etc. to results.
+    * Improved control over when to calculate the Hessian matrix.
+    * Fixed an issue converting the units of the gradients.
+
 2024.8.2.1 -- Calculator required to have 'implemented_properties' and 'nolabel' defined
 
 2024.8.2 -- Initial version

seamm_geometric/seamm_geometric.py

@@ -140,9 +140,9 @@ def calculate_gradients(self, coordinates):
         n_atoms = self._working_configuration.n_atoms
 
         if self.logger.isEnabledFor(logging.DEBUG):
-            print("\nnew coordinates")
+            self.logger.debug("\nnew coordinates")
             for i in range(n_atoms):
-                print(
+                self.logger.debug(
                     f"   {coordinates[i][0]:8.3f} {coordinates[i][1]:8.3f} "
                     f"{coordinates[i][2]:8.3f}"
                 )
@@ -213,6 +213,7 @@ def calculate_gradients(self, coordinates):
             traceback.print_exc(file=sys.stdout)
         except Exception as e:
             printer.job(f"Caught exception in step {self.step}: {str(e)}")
+            traceback.print_exc(file=sys.stdout)
             with open(step_dir / "stderr.out", "a") as fd:
                 traceback.print_exc(file=fd)
             raise
@@ -270,12 +271,29 @@ def calculate_gradients(self, coordinates):
                 )
 
         if "gradients,units" in data:
-            units = data["gradients,units"]
+            funits = data["gradients,units"]
+        else:
+            funits = "kJ/mol/Å"
+
+        # Get the measures of convergence
+        max_force = np.max(np.linalg.norm(gradients, axis=1))
+        self._data["max_force"].append(max_force)
+        self._results["maximum_gradient"] = Q_(max_force, funits).m_as("kJ/mol/Å")
+        rms_force = np.sqrt(np.mean(np.linalg.norm(gradients, axis=1) ** 2))
+        self._data["rms_force"].append(rms_force)
+        self._results["rms_gradient"] = Q_(rms_force, funits).m_as("kJ/mol/Å")
+
+        if self._step > 1:
+            step = coordinates - self._last_coordinates
+            max_step = np.max(np.linalg.norm(step, axis=1))
         else:
-            units = "kJ/mol/Å"
+            max_step = 0.0
+        self._data["max_step"].append(max_step)
+        self._results["maximum_step"] = max_step
+        self._last_coordinates = np.array(coordinates)
 
         # Units!
-        gradients = np.array(gradients) * Q_(1.0, units).to("E_h/a_0").magnitude
+        gradients = np.array(gradients) * Q_(1.0, funits).to("E_h/a_0").magnitude
 
         return energy, gradients
 
@@ -309,7 +327,7 @@ def describe_geomeTRIC_optimizer(self, P=None, short=False, natoms=None):
 
         if target == "minimum":
             text = "The structure will be optimized"
-        elif target == "transition_state":
+        elif target == "transition state":
             text = "The transition state will be optimized"
         elif self.is_expr(target):
             text = (
@@ -320,7 +338,18 @@ def describe_geomeTRIC_optimizer(self, P=None, short=False, natoms=None):
             raise ValueError(f"Unknown target {target}")
 
         text += " using the geomeTRIC optimizer with {coordinate system}. "
-        text += "Convergence will be reached when "
+
+        if self.is_expr(P["calculate hessian"]):
+            text += (
+                "'{calculate hessian}' will determine whether and how often to "
+                "calculate the Hessian matrix."
+            )
+        elif P["calculate hessian"] != "never":
+            text += (
+                "The Hessian matrix will be calculated for {calculate hessian} step."
+            )
+
+        text += " Convergence will be reached when "
         tmp_text = ["\n"]
         criteria = (
             "energy change criterion",
@@ -350,7 +379,7 @@ def describe_geomeTRIC_optimizer(self, P=None, short=False, natoms=None):
                 tmp_text.append(line.center(70))
             criteria = formula["criteria"]
 
-        result = str(__(text, **P))
+        result = str(__(text, dedent=False, indent=3 * " ", **P))
         result += "\n".join(tmp_text)
 
         table = {
@@ -365,23 +394,10 @@ def describe_geomeTRIC_optimizer(self, P=None, short=False, natoms=None):
 
         tmp = tabulate(table, headers="keys", tablefmt="rounded_outline")
 
-        text = "\n"
-        text += "\n"
-        text += tmp
-        text += "\n"
-        text += "\n"
-        result += str(__(text, indent=8 * " ", wrap=False, dedent=False))
-
-        if self.is_expr(P["calculate hessian"]):
-            result += (
-                "\n\n'{calculate hession}'will determine whether and how often to "
-                "calculate the Hessian matrix."
-            )
-        elif P["calculate hessian"] != "never":
-            result += (
-                "\n\nThe Hessian matrix will be calculated for "
-                "{calculate hessian} step."
-            )
+        result += "\n"
+        result += str(__(tmp, indent=7 * " ", wrap=False))
+        result += "\n"
+        result += "\n"
 
         return result
 
@@ -395,6 +411,16 @@ def run_geomeTRIC_optimizer(self, P, PP):
         PP : dict
             The current values of the parameters, formatted for printing
         """
+        self._data = {
+            "step": [],
+            "energy": [],
+            "max_force": [],
+            "rms_force": [],
+            "max_step": [],
+        }
+        self._last_coordinates = None
+        self._step = 0
+
         # Create the directory
         directory = Path(self.directory)
         self._working_directory = directory / "geomeTRIC"
@@ -487,7 +513,7 @@ def run_geomeTRIC_optimizer(self, P, PP):
             "maxiter": max_steps,
             "hessian": P["calculate hessian"],
             "frequency": P["calculate hessian"] != "never",
-            "transition": P["target"] == "Transition state",
+            "transition": "transition" in P["target"].lower(),
             "coordsys": coordsys,
         }