From 7d56506a57a2e71ed2d0fb5b6cc50100b04ee724 Mon Sep 17 00:00:00 2001 From: Lukas Dorschel Date: Mon, 25 Nov 2024 08:26:09 +0100 Subject: [PATCH 1/4] Updated interaction File --- src/python/espressomd/interactions.py | 38 ++++++++++++++++++++++++++- 1 file changed, 37 insertions(+), 1 deletion(-) diff --git a/src/python/espressomd/interactions.py b/src/python/espressomd/interactions.py index 9b95b78fef..196303c9a2 100644 --- a/src/python/espressomd/interactions.py +++ b/src/python/espressomd/interactions.py @@ -19,6 +19,8 @@ import abc import enum +from sympy import sympify, symbols +import numpy as np from . import code_features from .script_interface import ScriptObjectMap, ScriptInterfaceHelper, script_interface_register @@ -57,6 +59,41 @@ 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) + def get_table(self, **kwargs): + min_val = kwargs.pop("min") + max_val = kwargs.pop("max") + steps = kwargs.pop("steps") + expression = kwargs.pop("f") + r = symbols("r") + energy = sympify(expression) + force = -1 * energy.diff(r) + + for i in kwargs.keys(): + if str(i) in expression: + j = symbols(i) + energy = energy.subs(j, kwargs[i]) + force = force.subs(j, kwargs[i]) + + energy_tab = [] + force_tab = [] + + for x in np.linspace(min_val, max_val, steps): + try: + energy_value = energy.subs(r, x / steps) + force_value = force.subs(r, x / steps) + if energy_value.is_real and energy_value == energy_value: + energy_tab.append(energy_value) + if force_value.is_real and force_value == force_value: + force_tab.append(force_value) + + except (ZeroDivisionError, ValueError): + continue + if "shift" in kwargs.keys(): + j = kwargs.pop("shift") + energy_tab = [float(x*j) for x in energy_tab] + force_tab = [float(x*j) for x in force_tab] + return energy_tab, force_tab + @abc.abstractmethod def default_params(self): pass @@ -341,7 +378,6 @@ class TabulatedNonBonded(NonBondedInteraction): The force table. """ - _so_name = "Interactions::InteractionTabulated" _so_feature = "TABULATED" From 6ec719e5cead0977cf4262fa4597ef6d109d99b1 Mon Sep 17 00:00:00 2001 From: Lukas Dorschel Date: Mon, 25 Nov 2024 08:52:26 +0100 Subject: [PATCH 2/4] Original file --- src/python/espressomd/interactions.py | 38 +-------------------------- 1 file changed, 1 insertion(+), 37 deletions(-) diff --git a/src/python/espressomd/interactions.py b/src/python/espressomd/interactions.py index 196303c9a2..9b95b78fef 100644 --- a/src/python/espressomd/interactions.py +++ b/src/python/espressomd/interactions.py @@ -19,8 +19,6 @@ import abc import enum -from sympy import sympify, symbols -import numpy as np from . import code_features from .script_interface import ScriptObjectMap, ScriptInterfaceHelper, script_interface_register @@ -59,41 +57,6 @@ 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) - def get_table(self, **kwargs): - min_val = kwargs.pop("min") - max_val = kwargs.pop("max") - steps = kwargs.pop("steps") - expression = kwargs.pop("f") - r = symbols("r") - energy = sympify(expression) - force = -1 * energy.diff(r) - - for i in kwargs.keys(): - if str(i) in expression: - j = symbols(i) - energy = energy.subs(j, kwargs[i]) - force = force.subs(j, kwargs[i]) - - energy_tab = [] - force_tab = [] - - for x in np.linspace(min_val, max_val, steps): - try: - energy_value = energy.subs(r, x / steps) - force_value = force.subs(r, x / steps) - if energy_value.is_real and energy_value == energy_value: - energy_tab.append(energy_value) - if force_value.is_real and force_value == force_value: - force_tab.append(force_value) - - except (ZeroDivisionError, ValueError): - continue - if "shift" in kwargs.keys(): - j = kwargs.pop("shift") - energy_tab = [float(x*j) for x in energy_tab] - force_tab = [float(x*j) for x in force_tab] - return energy_tab, force_tab - @abc.abstractmethod def default_params(self): pass @@ -378,6 +341,7 @@ class TabulatedNonBonded(NonBondedInteraction): The force table. """ + _so_name = "Interactions::InteractionTabulated" _so_feature = "TABULATED" From 2cd586d2b8a1dff08987456fb595fdb2592ec5e9 Mon Sep 17 00:00:00 2001 From: Lukas Dorschel Date: Sun, 15 Dec 2024 13:50:33 +0100 Subject: [PATCH 3/4] Include test script --- testsuite/python/tabulated_sympy.py | 70 +++++++++++++++++++++++++++++ 1 file changed, 70 insertions(+) create mode 100644 testsuite/python/tabulated_sympy.py diff --git a/testsuite/python/tabulated_sympy.py b/testsuite/python/tabulated_sympy.py new file mode 100644 index 0000000000..7c49543b26 --- /dev/null +++ b/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 . +# +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() \ No newline at end of file From ac5d4ace5badded490903a04663b41843a727111 Mon Sep 17 00:00:00 2001 From: Lukas Dorschel Date: Sun, 15 Dec 2024 13:53:52 +0100 Subject: [PATCH 4/4] Includes interpretation of analytical functions --- src/python/espressomd/interactions.py | 79 ++++++++++++++++++++++++++- 1 file changed, 78 insertions(+), 1 deletion(-) diff --git a/src/python/espressomd/interactions.py b/src/python/espressomd/interactions.py index 9b95b78fef..ce94c72f86 100644 --- a/src/python/espressomd/interactions.py +++ b/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):