Sympy potential

[ludogibbs]

Added feature:
Generating tabulated forces and energies using analytical functions
Giving a potential V(r,params),params and number of steps instead of force and energy tables can now also be interpreted correctly.

[espresso-ci]

Your pull request does not meet our code formatting rules. Specifically, I suggest you make the following changes:

diff --git a/src/python/espressomd/interactions.py b/src/python/espressomd/interactions.py
index ce94c72f8..fca4b622a 100644
--- a/src/python/espressomd/interactions.py
+++ b/src/python/espressomd/interactions.py
@@ -19,7 +19,7 @@
 
 import abc
 import enum
-from sympy import sympify, symbols,oo,nan
+from sympy import sympify, symbols, oo, nan
 import numpy as np
 from . import code_features
 from .script_interface import ScriptObjectMap, ScriptInterfaceHelper, script_interface_register
@@ -59,7 +59,6 @@ class NonBondedInteraction(ScriptInterfaceHelper, metaclass=abc.ABCMeta):
         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
@@ -352,7 +351,7 @@ class TabulatedNonBonded(NonBondedInteraction):
 
         """
         return {}
-    
+
     def set_params(self, **kwargs):
         """Set new parameters.
 
@@ -366,16 +365,18 @@ class TabulatedNonBonded(NonBondedInteraction):
             if missing_keys:
                 raise ValueError(
                     f"Missing keys for table generation: {missing_keys}. "
-                    f"To generate energy and force tables, provide {required_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
+                # Zusätzliche Parameter
+                **{k: v for k, v in params.items() if k not in required_keys}
             )
             params["energy"] = energy_tab
             params["force"] = force_tab
@@ -384,9 +385,8 @@ class TabulatedNonBonded(NonBondedInteraction):
         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)
-
-        
+        self.call_method(
+            "set_params", handle_errors_message=err_msg, **filtered_params)
 
     def get_table(self, **kwargs):
         """
diff --git a/testsuite/python/tabulated_sympy.py b/testsuite/python/tabulated_sympy.py
index 7c49543b2..1cd80fec9 100644
--- a/testsuite/python/tabulated_sympy.py
+++ b/testsuite/python/tabulated_sympy.py
@@ -23,6 +23,7 @@ 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
@@ -34,8 +35,10 @@ class TabulatedTest(ut.TestCase):
         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.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])
 
@@ -47,13 +50,13 @@ class TabulatedTest(ut.TestCase):
         # 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_)):
+        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.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])
@@ -61,10 +64,11 @@ class TabulatedTest(ut.TestCase):
     @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)")
+            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()
\ No newline at end of file
+    ut.main()

To apply these changes, please do one of the following:

• You can download a patch with my suggested changes here, inspect it and make changes manually.
• You can directly apply it to your repository by running curl https://gitlab.icp.uni-stuttgart.de/espressomd/espresso/-/jobs/376854/artifacts/raw/style.patch | git apply -.
• You can run maintainer/CI/fix_style.sh to automatically fix your coding style. This is the same command that I have executed to generate the patch above, but it requires certain tools to be installed on your computer.

You can run gitlab-runner exec docker style afterwards to check if your changes worked out properly.

Please note that there are often multiple ways to correctly format code. As I am just a robot, I sometimes fail to identify the most aesthetically pleasing way. So please look over my suggested changes and adapt them where the style does not make sense.

[espresso-ci]

Your pull request does not meet our code style rules. Pylint summary:

<unknown>:246: SyntaxWarning: invalid escape sequence '\*'
************* Module tabulated_sympy
testsuite/python/tabulated_sympy.py:24:0: W0611: Unused symbols imported from sympy (unused-import)
testsuite/python/tabulated_sympy.py:24:0: W0611: Unused sympify imported from sympy (unused-import)

You can generate these warnings with maintainer/CI/fix_style.sh. This is the same command that I have executed to generate the log above.

[RudolfWeeber]

Thanks, Lukas. Plesae address the following:

• Add your test case to the CMakeList.txt in testsuite/python, so it gets actually run in CI.
• translate the comments in your code to English
• pls use a more meaningful arument name than f for the energy expression
• I like the string-based expression input, but it does not cover all cases (e.g. piecewhise functions or other more unusual sympy expressions). So please keep that but also accept a sympy expression in additoin to a string for the argument.
• A furhter argument distance_symbol should be added, so people can use something else than r. Keeping r as default is fine, I guess.

[jngrad]

Thanks! Don't forget to add the test to the CMake build system and to install and run the linters on your workstation. Here are the relevant sections in the developer guide:
Linter and code formatter, Python integration tests.

[jngrad]

Consider writing import sympy or import sympy as sp to limit module namespace pollution.

[jngrad]

Also: the present code makes sympy a hard requirement, even when TABULATED is not compiled in. Consider using a deferred import.

[jngrad]

I'm a bit uneasy about dynamically changing the set of expected arguments of a function depending on whether an argument is missing. From a user experience perspective, it's unclear what the function actually wants, or what happens when forces, energies, and a function are provided simultaneously. This also introduces a maintenance issue: the logic from the parent class set_params() method is duplicated here, thus if the parent class behavior changes, it won't be reflected here and TabulatedNonBonded will diverge from the rest of the code.

What would you think about creating a static class method like make_tables(cls, f, min, max, steps) that would return the energy and force tables, which can then be fed to set_params(self, **kwargs)?

[jngrad]

Here and elsewhere: please write comments in English language.

[jngrad]

The license header should reflect the year range you actually worked on the code, plus the years of the original code if you copied parts of it from somewhere else.

[jngrad]

sympy should be available in some of the CI containers. You can disable TabulatedTest when sympy is missing with fixture @utx.skipIfMissingModules("sympy") above the class name. Don't forget to add a linter rule exception for unused imports.

[jngrad]

"...from a function :math:`f(r)` that evaluates the force magnitude for a given inter-particle distance :math:`r`."