Use 2019 SI base units (#86)

* Use 2019 SI base units

* Remove SI units global state, add tests

* Rename constant Kb to kB

lib/handy_functions.py

@@ -19,32 +19,31 @@
 def setup_electrostatic_interactions (units, espresso_system, kT, c_salt=None, solvent_permittivity=78.5, method='p3m', tune_p3m=True, accuracy=1e-3,verbose=True):
     """
     Sets up electrostatic interactions in espressomd. 
-    Inputs:
-    system: instance of espressmd system class
-    c_salt: Added salt concentration. If provided, the program outputs the debye screening length. It is a mandatory parameter for the Debye-Huckel method. 
-    solvent_permittivity: Solvent relative permitivity, by default chosen per water at 298.15 K
-    method: method prefered for computing the electrostatic interactions. Currently only P3M (label = p3m) and Debye-Huckel (label = DH) are implemented
-    tune_p3m: If true (default), tunes the p3m parameters to improve efficiency
-    accuracy: desired accuracy for electrostatics, by default 1e-3
-    verbose (`bool`): switch to activate/deactivate verbose. Defaults to True.
+
+    Args:
+        units(`pint.UnitRegistry`): Unit registry.
+        espresso_system: instance of espressmd system class
+        kT(`float`): Thermal energy.
+        c_salt: Added salt concentration. If provided, the program outputs the debye screening length. It is a mandatory parameter for the Debye-Huckel method.
+        solvent_permittivity: Solvent relative permitivity, by default chosen per water at 298.15 K
+        method: method prefered for computing the electrostatic interactions. Currently only P3M (label = p3m) and Debye-Huckel (label = DH) are implemented
+        tune_p3m: If true (default), tunes the p3m parameters to improve efficiency
+        accuracy: desired accuracy for electrostatics, by default 1e-3
+        verbose (`bool`): switch to activate/deactivate verbose. Defaults to True.
     """
     import espressomd.electrostatics
     import numpy as np
-    #Initial checks
+    import scipy.constants
 
+    # Initial checks
     valid_methods_list=['p3m', 'DH']
-
     if method not in valid_methods_list:
-
         raise ValueError('provided an unknown label for method, valid values are', valid_methods_list)
-
     if c_salt is None and method == 'DH':
+        raise ValueError('Please provide the added salt concentration c_salt to setup the Debye-Huckel potential')
 
-        raise ValueError('Please provide the added salt concentration c_salt to settup the Debye-Huckel potential')
-
-
-    e = 1.60217662e-19 *units.C
-    N_A=6.02214076e23    / units.mol
+    e = scipy.constants.e * units.C
+    N_A = scipy.constants.N_A / units.mol
 
     BJERRUM_LENGTH = e.to('reduced_charge')**2 / (4 * units.pi * units.eps0 * solvent_permittivity * kT.to('reduced_energy'))
     if verbose:
@@ -53,19 +52,12 @@ def setup_electrostatic_interactions (units, espresso_system, kT, c_salt=None, s
     COULOMB_PREFACTOR=BJERRUM_LENGTH.to('reduced_length') * kT.to('reduced_energy') 
 
     if c_salt is not None:
-
         if c_salt.check('[substance] [length]**-3'):
-
             KAPPA=1./np.sqrt(8*units.pi*BJERRUM_LENGTH*N_A*c_salt)
-
         elif c_salt.check('[length]**-3'):
-
             KAPPA=1./np.sqrt(8*units.pi*BJERRUM_LENGTH*c_salt)
-
         else:
-
             raise ValueError('Unknown units for c_salt, please provided it in [mol / volume] or [particle / volume]', c_salt)
-
         if verbose:
             print(f"Debye kappa {KAPPA.to('nm')} = {KAPPA.to('reduced_length')}")
     print()

pyMBE.py

@@ -23,6 +23,7 @@
 import pint
 import numpy as np
 import pandas as pd
+import scipy.constants
 import scipy.optimize
 
 
@@ -39,10 +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.
     """
-    units = pint.UnitRegistry()
-    N_A=6.02214076e23    / units.mol
-    Kb=1.38064852e-23    * units.J / units.K
-    e=1.60217662e-19 *units.C
     df=None
     kT=None
     Kw=None
@@ -2165,7 +2162,7 @@ def get_reduced_units(self):
                                        f"{unit_length.to('nm'):.5g} = {unit_length}",
                                        f"{unit_energy.to('J'):.5g} = {unit_energy}",
                                        f"{unit_charge.to('C'):.5g} = {unit_charge}",
-                                       f"Temperature: {(self.kT/self.Kb).to('K'):.5g}"
+                                       f"Temperature: {(self.kT/self.kB).to('K'):.5g}"
                                         ])   
         return reduced_units_text
 
@@ -2720,10 +2717,10 @@ def set_reduced_units(self, unit_length=None, unit_charge=None, temperature=None
             unit_charge=self.units.e
         if Kw is None:
             Kw = 1e-14
-        self.N_A=6.02214076e23 / self.units.mol
-        self.Kb=1.38064852e-23 * self.units.J / self.units.K
-        self.e=1.60217662e-19 *self.units.C
-        self.kT=temperature*self.Kb
+        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
         self.Kw=Kw*self.units.mol**2 / (self.units.l**2)
         self.units.define(f'reduced_energy = {self.kT} ')
         self.units.define(f'reduced_length = {unit_length}')

testsuite/serialization_test.py

@@ -22,6 +22,7 @@
 import pandas as pd
 import pyMBE
 import lib.analysis
+import scipy.constants
 
 
 class Serialization(ut.TestCase):
@@ -60,6 +61,13 @@ def test_pint_units(self):
         pmb = pyMBE.pymbe_library(seed=42)
         reduced_units = pmb.get_reduced_units()
         self.assertEqual(reduced_units, ref_output)
+        np.testing.assert_allclose(
+            [pmb.kB.magnitude, pmb.N_A.magnitude, pmb.e.magnitude],
+            [scipy.constants.k, scipy.constants.N_A, scipy.constants.e],
+            rtol=1e-8, atol=0.)
+        self.assertAlmostEqual((pmb.kT / pmb.kB).magnitude, 298.15, delta=1e-7)
+        self.assertAlmostEqual((pmb.kT / scipy.constants.k).magnitude, 298.15,
+                               delta=1e-7)
 
 
 if __name__ == "__main__":