Merge pull request #301 from libAtoms/abort_low_energy

Check and exit when MD energy drops too much

tests/test_md.py

@@ -12,7 +12,7 @@
 
 from wfl.generate import md
 from wfl.configset import ConfigSet, OutputSpec
-from wfl.generate.md.abort import AbortOnCollision
+from wfl.generate.md.abort import AbortOnCollision, AbortOnLowEnergy
 
 
 def select_every_10_steps_for_tests_during(at):
@@ -191,16 +191,28 @@ def test_subselector_function_during(cu_slab):
 def test_md_abort_function(cu_slab):
 
     calc = EMT()
+    autopara_info = autoparainfo.AutoparaInfo(skip_failed=False)
 
+    # test bulk Cu slab for collision
     inputs = ConfigSet(cu_slab)
     outputs = OutputSpec()
-
     md_stopper = AbortOnCollision(collision_radius=2.25)
-    autopara_info = autoparainfo.AutoparaInfo(skip_failed=False)
 
     # why doesn't this throw an raise a RuntimeError even if md failed and `skip_failed` is False?
     atoms_traj = md.md(inputs, outputs, calculator=calc, steps=500, dt=10.0,
                        temperature=2000.0, abort_check=md_stopper, autopara_info=autopara_info,
                        rng=np.random.default_rng(1))
 
     assert len(list(atoms_traj)) < 501
+
+    # test surface Cu slab for low energy
+    cu_slab.set_cell(cu_slab.cell * (1, 1, 2), False)
+    inputs = ConfigSet(cu_slab)
+    outputs = OutputSpec()
+
+    md_stopper = AbortOnLowEnergy(0.0007)
+    atoms_traj = md.md(inputs, outputs, calculator=calc, steps=500, dt=10.0,
+                       temperature=5.0, abort_check=md_stopper, autopara_info=autopara_info,
+                       rng=np.random.default_rng(1))
+
+    assert len(list(atoms_traj)) < 501

wfl/generate/md/abort.py

@@ -2,6 +2,8 @@
 run under specified conditions
 """
 
+import numpy as np
+
 from .abort_base import AbortSimBase
 from ase.neighborlist import neighbor_list
 
@@ -16,18 +18,42 @@ class AbortOnCollision(AbortSimBase):
         distance for atoms to be considered a collision
 
     n_failed_steps: int, default 1
-        how many steps in a row any atom pairs have to be too cloe
+        how many steps in a row any atom pairs have to be too close
     """
 
     def __init__(self, collision_radius, n_failed_steps=3):
         super().__init__(n_failed_steps)
         self.collision_radius = collision_radius
 
 
-    def atoms_ok(self, at):
-        i = neighbor_list('i', at, self.collision_radius)
+    def atoms_ok(self, atoms):
+        i = neighbor_list('i', atoms, self.collision_radius)
 
         if len(i) > 0:
             return False
         else:
             return True
+
+
+class AbortOnLowEnergy(AbortSimBase):
+    """Abort an MD run if the energy drops by too much
+
+    Parameters
+    ----------
+    delta_E_per_atom: float
+        drop in energy per atom to trigger abort
+    """
+
+    def __init__(self, delta_E_per_atom):
+        super().__init__(1)
+        self.delta_E_per_atom = -np.abs(delta_E_per_atom)
+        self.initial_E_per_atom = None
+
+
+    def atoms_ok(self, atoms):
+        E_per_atom = atoms.get_potential_energy() / len(atoms)
+        if self.initial_E_per_atom is None:
+            self.initial_E_per_atom = E_per_atom
+            return True
+        else:
+            return (E_per_atom - self.initial_E_per_atom) >= self.delta_E_per_atom