Remove global state variables (#89)

* remove all global variables, avoid creating a new unit registry when redefining the set of reduced units

* fix missplaced Kw

pyMBE.py

@@ -40,12 +40,6 @@ class pymbe_library():
         kT(`pint.Quantity`): Thermal energy.
         Kw(`pint.Quantity`): Ionic product of water. Used in the setup of the G-RxMC method.
     """
-    df=None
-    kT=None
-    Kw=None
-    seed=None
-    rng=None
-
 
     class NumpyEncoder(json.JSONEncoder):
         """
@@ -79,8 +73,15 @@ def __init__(self, seed, temperature=None, unit_length=None, unit_charge=None, K
         # Seed and RNG
         self.seed=seed
         self.rng = np.random.default_rng(seed)
-        self.set_reduced_units(unit_length=unit_length, unit_charge=unit_charge,
-                               temperature=temperature, Kw=Kw, verbose=False)
+        self.units=pint.UnitRegistry()
+        self.N_A=scipy.constants.N_A / self.units.mol
+        self.kB=scipy.constants.k * self.units.J / self.units.K
+        self.e=scipy.constants.e * self.units.C
+        self.set_reduced_units(unit_length=unit_length, 
+                               unit_charge=unit_charge,
+                               temperature=temperature, 
+                               Kw=Kw, 
+                               verbose=False)
         self.setup_df()
         return
 
@@ -2708,20 +2709,22 @@ def set_reduced_units(self, unit_length=None, unit_charge=None, temperature=None
             - If no `unit_charge` is given, a value of 1 elementary charge is assumed by default. 
             - If no `Kw` is given, a value of 10^(-14) * mol^2 / l^2 is assumed by default. 
         """
-        self.units=pint.UnitRegistry()
         if unit_length is None:
-            unit_length=0.355*self.units.nm
+            unit_length= 0.355*self.units.nm
         if temperature is None:
-            temperature=298.15 * self.units.K
+            temperature = 298.15 * self.units.K
         if unit_charge is None:
-            unit_charge=self.units.e
+            unit_charge = scipy.constants.e * self.units.C
         if Kw is None:
             Kw = 1e-14
-        self.N_A=scipy.constants.N_A / self.units.mol
-        self.kB=scipy.constants.k * self.units.J / self.units.K
-        self.e=scipy.constants.e * self.units.C
-        self.kT=temperature*self.kB
+        # Sanity check
+        variables=[unit_length,temperature,unit_charge]
+        dimensionalities=["[length]","[temperature]","[charge]"]
+        for variable,dimensionality in zip(variables,dimensionalities):
+            self.check_dimensionality(variable,dimensionality)
         self.Kw=Kw*self.units.mol**2 / (self.units.l**2)
+        self.kT=temperature*self.kB
+        self.units._build_cache()
         self.units.define(f'reduced_energy = {self.kT} ')
         self.units.define(f'reduced_length = {unit_length}')
         self.units.define(f'reduced_charge = {unit_charge}')

testsuite/CTestTestfile.cmake

@@ -48,6 +48,7 @@ pymbe_add_test(PATH gcmc_tests.py LABELS long)
 
 # unit tests
 pymbe_add_test(PATH serialization_test.py)
+pymbe_add_test(PATH test_global_variables.py)
 pymbe_add_test(PATH lj_tests.py)
 pymbe_add_test(PATH set_particle_acidity_test.py)
 pymbe_add_test(PATH bond_tests.py)

testsuite/create_molecule_position_test.py

@@ -22,7 +22,7 @@
 import pyMBE
 pmb = pyMBE.pymbe_library(seed=42)
 
-print("***create_molecule with input position list unit test ***")
+print("*** Create_molecule with input position list unit test ***")
 print("*** Unit test: Check that the positions of the central bead of the first residue in the generated molecules are equal to the input positions ***")
 # Simulation parameters
 pmb.set_reduced_units(unit_length=0.4*pmb.units.nm,

testsuite/test_global_variables.py

@@ -0,0 +1,51 @@
+#
+# Copyright (C) 2024 pyMBE-dev team
+#
+# This file is part of pyMBE.
+#
+# pyMBE 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.
+#
+# pyMBE 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 pyMBE
+import numpy as np
+
+# Test that two different instances of pyMBE do not share the same memory address for their global attributes
+pmb1 = pyMBE.pymbe_library(seed=42)
+pmb2 = pyMBE.pymbe_library(seed=43)
+np.testing.assert_raises(AssertionError,
+                        np.testing.assert_equal,
+                        id(pmb1.units),
+                        id(pmb2.units))
+
+# Test that redefining the system of reduced units does not create a new pint.UnitRegistry
+
+## Define a variables in the old unit registry
+nm_length_1=pmb1.units.Quantity(1,"nm")
+
+## Change the system of reduced units
+pmb1.set_reduced_units(unit_length=0.4*pmb1.units.nm)
+
+## Define variable in the new unit registry
+nm_length_2=pmb1.units.Quantity(2,"nm")
+
+## operations between old and new quantities should work normally 
+np.testing.assert_equal((nm_length_1+nm_length_2).m_as("nm"), 
+                        (pmb1.units.Quantity(3,"nm").m_as("nm")))
+
+
+# Test that set_reduced_units raises a ValueError if the wrong unit is provided
+input_parameters={"unit_length": pmb1.units.Quantity(1, "J")}                      
+
+np.testing.assert_raises(ValueError,
+                        pmb1.set_reduced_units,
+                        **input_parameters)