Allow modifying cutoff values from ASE calculator (#154)

Author: awvwgk

python/dftd3/ase.py

@@ -39,7 +39,8 @@
 ======================== ============ ============================================
  method                   None         Method to calculate dispersion for
  damping                  None         Damping function to use
- params_tweaks            None         Optional dict with the damping parameters
+ params_tweaks            {}           Optional dict with the damping parameters
+ realspace_cutoff         {}           Optional dict to override cutoff values
  cache_api                True         Reuse generate API objects (recommended)
 ======================== ============ ============================================
 
@@ -63,6 +64,20 @@
 for damping parameters of a given method and prefers special two-body only
 damping parameters if available!
 
+The realspace cutoff parameters allow adjusting the distance values for which
+interactions are considered
+
+================== =========== ==========================================
+ Realspace cutoff   Default     Description
+================== =========== ==========================================
+ disp2              60 * Bohr   Pairwise dispersion interactions
+ disp3              40 * Bohr   Triple dispersion interactions
+ cn                 40 * Bohr   Coordination number counting
+================== =========== ==========================================
+
+Values provided in the dict are expected to be in Angstrom. When providing values
+in Bohr multiply the inputs by the `ase.units.Bohr` constant.
+
 Example
 -------
 >>> from ase.build import molecule
@@ -146,6 +161,7 @@ class DFTD3(Calculator):
         "method": None,
         "damping": None,
         "params_tweaks": {},
+        "realspace_cutoff": {},
         "cache_api": True,
     }
 
@@ -237,20 +253,33 @@ def _create_api_calculator(self) -> DispersionModel:
                 _periodic,
             )
 
-        except RuntimeError:
-            raise InputError("Cannot construct dispersion model for dftd3")
+        except RuntimeError as e:
+            raise InputError("Cannot construct dispersion model for dftd3") from e
 
         return disp
 
+    def _apply_realspace_cutoff(self, disp: DispersionModel) -> None:
+        """Apply realspace cutoff parameters to dispersion model"""
+
+        try:
+            if self.parameters.realspace_cutoff:
+                disp2 = self.parameters.realspace_cutoff.get("disp2", 60.0 * Bohr) / Bohr
+                disp3 = self.parameters.realspace_cutoff.get("disp3", 40.0 * Bohr) / Bohr
+                cn = self.parameters.realspace_cutoff.get("cn", 40.0 * Bohr) / Bohr
+
+                disp.set_realspace_cutoff(disp2=disp2, disp3=disp3, cn=cn)
+        except RuntimeError as e:
+            raise InputError("Cannot update realspace cutoff for dftd3") from e
+
     def _create_damping_param(self) -> DampingParam:
         """Create a new API damping parameter object"""
 
         try:
             params_tweaks = self.parameters.params_tweaks if self.parameters.params_tweaks else {"method": self.parameters.get("method")} 
             dpar = _damping_param[self.parameters.get("damping")](**params_tweaks)
 
-        except RuntimeError:
-            raise InputError("Cannot construct damping parameter for dftd3")
+        except RuntimeError as e:
+            raise InputError("Cannot construct damping parameter for dftd3") from e
 
         return dpar
 
@@ -271,12 +300,15 @@ def calculate(
         if self._disp is None:
             self._disp = self._create_api_calculator()
 
+        # Apply realspace cutoff before evaluation (works with cached calculator)
+        self._apply_realspace_cutoff(self._disp)
+
         _dpar = self._create_damping_param()
 
         try:
             _res = self._disp.get_dispersion(param=_dpar, grad=True)
-        except RuntimeError:
-            raise CalculationFailed("dftd3 could not evaluate input")
+        except RuntimeError as e:
+            raise CalculationFailed("dftd3 could not evaluate input") from e
 
         # These properties are garanteed to exist for all implemented calculators
         self.results["energy"] = _res.get("energy") * Hartree

python/dftd3/test_ase.py

@@ -97,3 +97,54 @@ def test_ase_tpssd3():
     assert atoms.get_potential_energy() == approx(4.963774668847532, abs=thr)
     energies = [calc.get_potential_energy() for calc in get_calcs(atoms.calc)]
     assert energies == approx([-0.14230914516094673, 5.106083814008478], abs=thr)
+
+
+@pytest.mark.skipif(ase is None, reason="requires ase")
+def test_ase_realspace_cutoff():
+    """Test that realspace_cutoff parameter works and can be updated with cache_api=True"""
+    from ase.units import Bohr
+    
+    thr = 1.0e-6
+    atoms = molecule("H2O")
+    
+    # Test with default cutoffs
+    calc_default = DFTD3(method="PBE", damping="d3bj", cache_api=True)
+    atoms.calc = calc_default
+    energy_default = atoms.get_potential_energy()
+    forces_default = atoms.get_forces()
+    
+    # Test with very small cutoffs (smaller than H2O bond length ~1 Angstrom) - should give zero or very small interactions
+    calc_custom = DFTD3(
+        method="PBE", 
+        damping="d3bj",
+        realspace_cutoff={"disp2": 0.5, "disp3": 0.5, "cn": 0.5},
+        cache_api=True
+    )
+    atoms.calc = calc_custom
+    energy_custom = atoms.get_potential_energy()
+    forces_custom = atoms.get_forces()
+    
+    # With very small cutoffs, energy should be much smaller (closer to zero)
+    assert abs(energy_custom) < abs(energy_default)
+    assert not np.allclose(forces_custom, forces_default, atol=thr)
+    
+    # Test updating cutoff via set() with cache_api=True
+    # Reset to default-like cutoffs using set()
+    calc_custom.set(realspace_cutoff={"disp2": 60.0 * Bohr, "disp3": 40.0 * Bohr, "cn": 40.0 * Bohr})
+    energy_updated = atoms.get_potential_energy()
+    forces_updated = atoms.get_forces()
+    
+    # Should match the default energy/forces
+    assert energy_updated == approx(energy_default, abs=thr)
+    assert forces_updated == approx(forces_default, abs=thr)
+    
+    # Test with empty dict (should use library defaults)
+    calc_empty = DFTD3(method="PBE", damping="d3bj", realspace_cutoff={}, cache_api=True)
+    atoms.calc = calc_empty
+    energy_empty = atoms.get_potential_energy()
+    
+    # Empty dict should behave like no cutoff override (same as default)
+    assert energy_empty == approx(energy_default, abs=thr)
+
+
+