PhononFromHessian class to save computing time on isotopes

mlcalcdriver/calculators/calculator.py

@@ -1,3 +1,5 @@
+from mlcalcdriver.globals import eVA
+
 r"""
 The :class:`Calculator` is the general class for a machine learning
 calculator. A specific class derived from this one must be implemented
@@ -84,3 +86,9 @@ def units(self, units):
                 raise KeyError("Units key not recognized.")
         else:
             raise TypeError("Units should be given in a dictionary.")
+
+
+class DummyCalculator(Calculator):
+    def __init__(self):
+        properties = ["energy"]
+        super().__init__(available_properties=properties, units=eVA)

mlcalcdriver/workflows/phonon.py

@@ -6,7 +6,7 @@
 
 import numpy as np
 from mlcalcdriver import Job, Posinp
-from mlcalcdriver.calculators import Calculator
+from mlcalcdriver.calculators.calculator import Calculator, DummyCalculator
 from mlcalcdriver.workflows import Geopt
 from copy import deepcopy
 from mlcalcdriver.globals import ANG_TO_B, B_TO_ANG, EV_TO_HA, HA_TO_CMM1, AMU_TO_EMU
@@ -294,7 +294,7 @@ def _compute_hessian(self, job):
             h = (
                 job.results["hessian"].reshape(3 * n_at, 3 * n_at)
                 * EV_TO_HA
-                * B_TO_ANG ** 2
+                * B_TO_ANG**2
             )
             return (h + h.T) / 2.0
         else:
@@ -316,3 +316,52 @@ def _solve_dyn_mat(self):
         eigs, vecs = np.linalg.eigh(self.dyn_mat)
         eigs = np.sign(eigs) * np.sqrt(np.where(eigs < 0, -eigs, eigs))
         return eigs, vecs
+
+
+class PhononFromHessian(Phonon):
+    r"""
+    Similar to the main Phonon class, but can be used when calculating
+    many structures with the same hessian matrix (isotopes study). Saves
+    the time to compute the hessian matrix each time
+    """
+
+    def __init__(self, posinp, hessian):
+        r"""
+        Parameters
+        ----------
+        posinp : mlcaldriver.Posinp
+            Initial positions of the system under consideration.
+        hessian : np.ndarray or str
+            Hessian matrix calculated before instanciating this class.
+            Can be the array or a path to .npy file (created with np.save).
+        """
+        super().__init__(
+            posinp=posinp,
+            calculator=DummyCalculator(),
+            relax=False,
+            finite_difference=False,
+        )
+        self.hessian = hessian
+
+    @property
+    def hessian(self):
+        return self._hessian
+
+    @hessian.setter
+    def hessian(self, hessian):
+        if isinstance(hessian, str):
+            hessian = np.load(hessian)
+
+        if isinstance(hessian, np.ndarray):
+            assert hessian[0].shape == (
+                3 * len(self.posinp),
+                3 * len(self.posinp),
+            ), f"The hessian shape {hessian.shape} does not match the number of atoms {len(self.posinp)}"
+            self._hessian = hessian
+        else:
+            raise TypeError("The hessian matrix should be a numpy array.")
+
+    def run(self):
+        job = Job(posinp=self.posinp, calculator=self.calculator)
+        job.results["hessian"] = self.hessian
+        self._post_proc(job)

tests/test_phonon.py

@@ -3,14 +3,13 @@
 import numpy as np
 from mlcalcdriver import Posinp, Job
 from mlcalcdriver.calculators import SchnetPackCalculator
-from mlcalcdriver.workflows import Phonon
+from mlcalcdriver.workflows.phonon import Phonon, PhononFromHessian
 
 pos_folder = "tests/posinp_files/"
 model_folder = "tests/models/"
 
 
 class TestPhononFinite:
-
     posN2 = Posinp.from_file(os.path.join(pos_folder, "N2_unrelaxed.xyz"))
     calcN2 = SchnetPackCalculator(os.path.join(model_folder, "myN2_model"))
 
@@ -49,7 +48,6 @@ def test_phonon_calc_error(self):
 
 
 class TestPhononAutoGrad:
-
     posH2O = Posinp.from_file(os.path.join(pos_folder, "H2Orelaxed.xyz"))
     calc_ener = SchnetPackCalculator(os.path.join(model_folder, "H2O_model"))
     calc_for = SchnetPackCalculator(os.path.join(model_folder, "H2O_forces_model"))
@@ -65,3 +63,11 @@ def test_ph_h2o_autograd_1st_derivative(self):
         ph1.run()
         ph1.energies.sort()
         assert np.allclose(ph1.energies[6:9], [1589, 3703, 3812], atol=1)
+
+    def test_ph_from_hessian(self):
+        job = Job(posinp=self.posH2O, calculator=self.calc_for)
+        job.run("hessian")
+        ph = PhononFromHessian(posinp=self.posH2O, hessian=job.results["hessian"])
+        ph.run()
+        ph.energies.sort()
+        assert np.allclose(ph.energies[6:9], [1589, 3705, 3814], atol=1)