Sympy potential

src/python/espressomd/interactions.py

@@ -19,6 +19,8 @@
 
 import abc
 import enum
+from sympy import sympify, symbols,oo,nan
+import numpy as np
 from . import code_features
 from .script_interface import ScriptObjectMap, ScriptInterfaceHelper, script_interface_register
 
@@ -57,6 +59,7 @@ def set_params(self, **kwargs):
         err_msg = f"setting {self.__class__.__name__} raised an error"
         self.call_method("set_params", handle_errors_message=err_msg, **params)
 
+
     @abc.abstractmethod
     def default_params(self):
         pass
@@ -341,7 +344,6 @@ class TabulatedNonBonded(NonBondedInteraction):
             The force table.
 
     """
-
     _so_name = "Interactions::InteractionTabulated"
     _so_feature = "TABULATED"
 
@@ -350,6 +352,81 @@ def default_params(self):
 
         """
         return {}
+
+    def set_params(self, **kwargs):
+        """Set new parameters.
+
+        """
+        params = self.default_params()
+        params.update(kwargs)
+
+        if "energy" not in params or "force" not in params:
+            required_keys = {"min", "max", "steps", "f"}
+            missing_keys = required_keys - params.keys()
+            if missing_keys:
+                raise ValueError(
+                    f"Missing keys for table generation: {missing_keys}. "
+                    f"To generate energy and force tables, provide {required_keys}."
+                )
+
+            # Berechnung von Energie- und Krafttabellen mit `get_table`
+            energy_tab, force_tab = self.get_table(
+                min=params["min"],
+                max=params["max"],
+                steps=params["steps"],
+                f=params["f"],
+                **{k: v for k, v in params.items() if k not in required_keys}  # Zusätzliche Parameter
+            )
+            params["energy"] = energy_tab
+            params["force"] = force_tab
+
+        relevant_keys = {"min", "max", "energy", "force"}
+        filtered_params = {k: params[k] for k in relevant_keys if k in params}
+
+        err_msg = f"setting {self.__class__.__name__} raised an error"
+        self.call_method("set_params", handle_errors_message=err_msg, **filtered_params)
+
+
+
+    def get_table(self, **kwargs):
+        """
+        Generate energy and force tables.
+        """
+        # Erforderliche Parameter extrahieren
+        try:
+            min_val = kwargs["min"]
+            max_val = kwargs["max"]
+            steps = kwargs["steps"]
+            expression = kwargs["f"]
+        except KeyError as e:
+            raise ValueError(f"Missing required parameter: {e.args[0]}") from e
+
+        # Symbole und Ausdrücke vorbereiten
+        r = symbols("r")
+        energy = sympify(expression)
+        force = -energy.diff(r)
+
+        # Zusätzliche Variablen ersetzen
+        for var, value in kwargs.items():
+            if str(var) in expression:
+                symbol = symbols(var)
+                energy = energy.subs(symbol, value)
+                force = force.subs(symbol, value)
+
+        # Tabellen berechnen
+        x_values = np.linspace(min_val, max_val, steps)
+        energy_tab = [
+            float(energy.subs(r, x))
+            for x in x_values
+            if energy.subs(r, x).is_real and not energy.subs(r, x).has(oo, nan)
+        ]
+        force_tab = [
+            float(force.subs(r, x))
+            for x in x_values
+            if force.subs(r, x).is_real and not force.subs(r, x).has(oo, nan)
+        ]
+
+        return energy_tab, force_tab
 
     @property
     def cutoff(self):

testsuite/python/tabulated_sympy.py

@@ -0,0 +1,70 @@
+#
+# Copyright (C) 2013-2022 The ESPResSo project
+#
+# This file is part of ESPResSo.
+#
+# ESPResSo is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# ESPResSo is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+import unittest as ut
+import unittest_decorators as utx
+import espressomd
+import espressomd.interactions
+import numpy as np
+from sympy import symbols, sympify
+
+class TabulatedTest(ut.TestCase):
+    system = espressomd.System(box_l=3 * [10.])
+    system.time_step = 0.01
+    system.cell_system.skin = 0.4
+
+    def setUp(self):
+        self.min_ = 1.
+        self.max_ = 2.
+        self.eps_ = 1.
+        self.sig_ = 2.
+        self.steps_ = 100
+        self.force = [-4*self.eps_*(12/r*(self.sig_/r)**12-6/r*(self.sig_/r)**6) for r in np.linspace(self.min_,self.max_,self.steps_)]
+        self.energy = [4*self.eps_*((self.sig_/r)**12-(self.sig_/r)**6) for r in np.linspace(self.min_,self.max_,self.steps_)]
+        self.system.part.add(type=0, pos=[5., 5., 5.0])
+        self.system.part.add(type=0, pos=[5., 5., 5.5])
+
+    def tearDown(self):
+        self.system.part.clear()
+
+    def check(self):
+        p0, p1 = self.system.part.all()
+        # Below cutoff
+        np.testing.assert_allclose(np.copy(self.system.part.all().f), 0.0)
+
+        for i,z in enumerate(np.linspace(0, self.max_ - self.min_, self.steps_)):
+            if z >= self.max_ - self.min_:
+                continue
+            p1.pos = [5., 5., 6. + z]
+            self.system.integrator.run(0)
+            np.testing.assert_allclose(
+                np.copy(p0.f), [0., 0.,self.force[i]])
+            np.testing.assert_allclose(np.copy(p0.f), -np.copy(p1.f))
+            self.assertAlmostEqual(
+                self.system.analysis.energy()['total'], self.energy[i])
+
+    @utx.skipIfMissingFeatures("TABULATED")
+    def test_tabulated_sympy(self):
+        self.system.non_bonded_inter[0, 0].tabulated.set_params(
+            min=self.min_, max=self.max_,steps=self.steps_,sig=self.sig_,eps=self.eps_, f="4*eps*((sig/r)**12-(sig/r)**6)")
+        self.assertEqual(
+            self.system.non_bonded_inter[0, 0].tabulated.cutoff, self.max_)
+        self.check()
+
+if __name__ == "__main__":
+    ut.main()