Merge pull request #933 from pyiron/correct_qha

replace linear regression by ridge regression

pyiron_contrib/atomistics/atomistics/master/qha.py

@@ -1,7 +1,8 @@
 import numpy as np
 import pint
-from sklearn.linear_model import LinearRegression
+from sklearn.linear_model import RidgeCV
 from pyiron_atomistics.atomistics.master.parallel import AtomisticParallelMaster
+import warnings
 from pyiron_base import JobGenerator
 
 
@@ -148,6 +149,7 @@ def __init__(
         symprec=1.0e-2,
         include_zero_displacement=True,
         second_order=True,
+        alphas=np.logspace(-20, 20, 50),
     ):
         """
         Calculation of the Hessian matrix
@@ -174,6 +176,7 @@ def __init__(
         self._hessian = None
         self.include_zero_displacement = include_zero_displacement
         self.second_order = second_order
+        self.alphas = alphas
 
     @property
     def symmetry(self):
@@ -288,12 +291,21 @@ def displacements(self, d):
         self._displacements = np.asarray(d).tolist()
         self._inequivalent_displacements = None
 
+    @property
+    def E_fx(self):
+        return np.repeat(
+            self.energy, len(self.symmetry.rotations)
+        )[self.inequivalent_indices]
+
     @property
     def _fit(self):
-        E = self.energy + 0.5 * np.einsum("nij,nij->n", self.forces, self.displacements)
-        E = np.repeat(E, len(self.symmetry.rotations))[self.inequivalent_indices]
-        reg = LinearRegression()
-        reg.fit(self.inequivalent_forces, E)
+        reg = RidgeCV(alphas=self.alphas)
+        reg.fit(self.inequivalent_forces, self.E_fx)
+        if np.linalg.norm(reg.coef_) > self.dx:
+            warnings.warn(
+                "Predicted minimum energy structure might be too far from the"
+                " original structure"
+            )
         return reg
 
     @property
@@ -304,7 +316,7 @@ def min_energy(self):
     def origin(self):
         if self.energy is None or self.forces is None:
             return np.zeros_like(self.structure.positions)
-        return 2 * self.symmetry.symmetrize_vectors(self._fit.coef_.reshape(-1, 3))
+        return self.symmetry.symmetrize_vectors(self._fit.coef_.reshape(-1, 3))
 
     @property
     def volume(self):