Merge pull request #1299 from pyiron/qn_separate_update

[minor] Separate update functions in Quasi Newton

pyiron_atomistics/interactive/quasi_newton.py

@@ -7,6 +7,52 @@
 )
 
 
+def get_SR(dx, dg, H, threshold=1e-4):
+    """
+    Args:
+        dx ((n,)-numpy.ndarray/list): Change in positions
+        dg ((n,)-numpy.ndarray/list): Change in derivatives
+        H ((n,n)-numpy.ndarray/list): Current Hessian
+        threshold (float): Minimum value that the denominator can have (the
+            resulting Hessian matrix might explode if too small)
+
+    Returns:
+        ((n,n)-numpy.ndarray): Updated Hessian matrix
+    """
+    H_tmp = dg - np.einsum("ij,j->i", H, dx)
+    denominator = np.dot(H_tmp, dx)
+    if np.absolute(denominator) < threshold:
+        denominator += threshold
+    return np.outer(H_tmp, H_tmp) / denominator + H
+
+
+def get_PSB(dx, dg, H):
+    """
+    Args:
+        dx ((n,)-numpy.ndarray/list): Change in positions
+        dg ((n,)-numpy.ndarray/list): Change in derivatives
+        H ((n,n)-numpy.ndarray/list): Current Hessian
+
+    Returns:
+        ((n,n)-numpy.ndarray): Updated Hessian matrix
+    """
+    H_tmp = dg - np.einsum("ij,j->i", H, dx)
+    dxdx = np.einsum("i,i->", dx, dx)
+    dH = np.einsum("i,j->ij", H_tmp, dx)
+    dH = (dH + dH.T) / dxdx
+    return (
+        dH - np.einsum("i,i,j,k->jk", dx, H_tmp, dx, dx, optimize="optimal") / dxdx**2
+    ) + H
+
+
+def get_BFGS(dx, dg, H):
+    Hx = H.dot(dx)
+    return np.outer(dg, dg) / dg.dot(dx) - np.outer(Hx, Hx) / dx.dot(Hx) + H
+
+
+get_BFGS.__doc__ = get_PSB.__doc__
+
+
 class QuasiNewtonInteractive:
     """
     Interactive class of Quasi Newton. This class can be used without a pyiron job definition.
@@ -144,42 +190,18 @@ def get_dx(self, g, threshold=1e-4, mode="PSB", update_hessian=True):
             )
         return self.dx
 
-    @staticmethod
-    def _get_SR(dx, dg, H_tmp, threshold=1e-4):
-        denominator = np.dot(H_tmp, dx)
-        if np.absolute(denominator) < threshold:
-            denominator += threshold
-        return np.outer(H_tmp, H_tmp) / denominator
-
-    @staticmethod
-    def _get_PSB(dx, dg, H_tmp):
-        dxdx = np.einsum("i,i->", dx, dx)
-        dH = np.einsum("i,j->ij", H_tmp, dx)
-        dH = (dH + dH.T) / dxdx
-        return (
-            dH
-            - np.einsum("i,i,j,k->jk", dx, H_tmp, dx, dx, optimize="optimal")
-            / dxdx**2
-        )
-
-    @staticmethod
-    def _get_BFGS(dx, dg, H):
-        Hx = H.dot(dx)
-        return np.outer(dg, dg) / dg.dot(dx) - np.outer(Hx, Hx) / dx.dot(Hx)
-
     def update_hessian(self, g, threshold=1e-4, mode="PSB"):
         if self.g_old is None:
             self.g_old = g
             return
         dg = self.get_dg(g).flatten()
         dx = self.dx.flatten()
-        H_tmp = dg - np.einsum("ij,j->i", self.hessian, dx)
         if mode == "SR":
-            self.hessian = self._get_SR(dx, dg, H_tmp) + self.hessian
+            self.hessian = get_SR(dx, dg, self.hessian)
         elif mode == "PSB":
-            self.hessian = self._get_PSB(dx, dg, H_tmp) + self.hessian
+            self.hessian = get_PSB(dx, dg, self.hessian)
         elif mode == "BFGS":
-            self.hessian = self._get_BFGS(dx, dg, self.hessian) + self.hessian
+            self.hessian = get_BFGS(dx, dg, self.hessian)
         else:
             raise ValueError(
                 "Mode not recognized: {}. Choose from `SR`, `PSB` and `BFGS`".format(