improved units managments

Author: simongravelle

docs/source/chapters/chapter1.rst

@@ -1,8 +1,13 @@
 .. _chapter1-label:
 
-Start coding
+Start Coding
 ============
 
+Let's start with the Python script. During this first chapter, some Python
+files will be created and filled in a minimal fashion. At the end of this
+chapter, a small test will be set up to ensure that the files were correctly
+created.
+
 Presentation
 ------------
 
@@ -38,7 +43,8 @@ containing either Python functions or classes:
      - *Prepare* class: Methods for preparing the non-dimensionalization of the
        units
    * - *Utilities.py* 
-     - *Utilities* class: General-purpose methods, inherited by all other classes
+     - *Utilities* class: General-purpose methods, inherited by all other
+       classes
    * - *InitializeSimulation.py*
      - *InitializeSimulation* class: Methods necessary to set up the system and
        prepare the simulation, inherited by all the classes below
@@ -54,21 +60,27 @@ containing either Python functions or classes:
      - Functions for performing specific measurements on the system
    * - *potentials.py* 
      - Functions for calculating the potentials and forces between atoms
-   * - *tools.py*
+   * - *logger.py*
      - Functions for outputting data into text files
+   * - *dumper.py*
+     - Functions for outputting data into trajectory files for visualization
+   * - *reader.py*
+     - Functions for importing data from text files
 
+Some of these files are created in this chapter; others will be created later
+on. All of these files must be created within the same folder.
 
-Potential for inter-atomic interaction
+Potential for Inter-Atomic Interaction
 --------------------------------------
 
 In molecular simulations, potential functions are used to mimic the interaction
 between atoms. Although more complicated options exist, potentials are usually
 defined as functions of the distance :math:`r` between two atoms.
 
-Within a dedicated folder, create the first file named *potentials.py*. This
-file will contain a function called *potentials*. Two types of potential can 
-be returned by this function: the Lennard-Jones potential (LJ), and the
-hard-sphere potential.
+Create the first file named *potentials.py*. This file will contain a function
+called *potentials*. For the moment, the only potential that can be returned by
+this function is the Lennard-Jones potential (LJ). This may change in the
+future.
 
 Copy the following lines into *potentials.py*:
 
@@ -84,44 +96,33 @@ Copy the following lines into *potentials.py*:
                 return 48 * epsilon * ((sigma / r) ** 12 - 0.5 * (sigma / r) ** 6) / r
             else:
                 return 4 * epsilon * ((sigma / r) ** 12 - (sigma / r) ** 6)
-        elif potential_type == "Hard-Sphere":
-            if derivative:
-                # Derivative is not defined for Hard-Sphere potential.
-                # --> return 0
-                return np.zeros(len(r))
-            else:
-                return np.where(r > sigma, 0, 1000)
         else:
             raise ValueError(f"Unknown potential type: {potential_type}")
 
 .. label:: end_potentials_class
 
-The hard-sphere potential either returns a value of 0 when the distance between
-the two particles is larger than the parameter, :math:`r > \sigma`, or 1000 when
-:math:`r < \sigma`. The value of *1000* was chosen to be large enough to ensure
-that any Monte Carlo move that would lead to the two particles to overlap will
-be rejected.
-
-In the case of the LJ potential, depending on the value of the optional
-argument *derivative*, which can be either *False* or *True*, the *LJ_potential*
-function will return the force:
+Depending on the value of the optional argument *derivative*, which can be
+either *False* or *True*, this function returns the derivative of the potential,
+i.e., the force, :math:`F_\text{LJ} = - \mathrm{d} U_\text{LJ} / \mathrm{d} r`:
 
 .. math::
 
-    F_\text{LJ} = 48 \dfrac{\epsilon}{r} \left[ \left( \frac{\sigma}{r} \right)^{12} - \frac{1}{2} \left( \frac{\sigma}{r} \right)^6 \right],
+    F_\text{LJ} = 48 \dfrac{\epsilon}{r} \left[ \left( \frac{\sigma}{r} \right)^{12}
+    - \frac{1}{2} \left( \frac{\sigma}{r} \right)^6 \right],
 
 or the potential energy:
 
 .. math::
 
-    U_\text{LJ} = 4 \epsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r} \right)^6 \right].
+    U_\text{LJ} = 4 \epsilon \left[ \left( \frac{\sigma}{r} \right)^{12}
+    - \left( \frac{\sigma}{r} \right)^6 \right].
 
 Create the Classes
 ------------------
 
-Let's create the files with the minimal information about the classes and
-their inheritance. The classes will be developed progressively in the
-following chapters.
+Let's create the files with some minimal details about the classes and their
+inheritance. The classes will be developed progressively in the following
+chapters.
 
 The first class is the *Prepare* class, which will be used for the
 nondimensionalization of the parameters. In the same folder as *potentials.py*,
@@ -157,8 +158,8 @@ copy the following lines:
 
 .. label:: end_Utilities_class
 
-The line *from potentials import LJ_potential* is used to import the
-*LJ_potential* function.
+The line *from potentials import potentials* is used to import the
+previously created *potentials* function.
 
 Within the *InitializeSimulation.py* file, copy the following lines:
 
@@ -168,9 +169,10 @@ Within the *InitializeSimulation.py* file, copy the following lines:
 
     import numpy as np
     from Prepare import Prepare
+    from Utilities import Utilities
 
 
-    class InitializeSimulation(Prepare):
+    class InitializeSimulation(Prepare, Utilities):
         def __init__(self,
                     *args,
                     **kwargs,
@@ -180,7 +182,8 @@ Within the *InitializeSimulation.py* file, copy the following lines:
 .. label:: end_InitializeSimulation_class
 
 The *InitializeSimulation* class inherits from the previously created
-*Prepare* class. Additionally, we anticipate that *NumPy* will be required.
+*Prepare* and Utilities classes. Additionally, we anticipate that *NumPy* will
+be required.
 
 Within the *Measurements.py* file, copy the following lines:
 
@@ -189,10 +192,9 @@ Within the *Measurements.py* file, copy the following lines:
 .. code-block:: python
 
     from InitializeSimulation import InitializeSimulation
-    from Utilities import Utilities
 
 
-    class Measurements(InitializeSimulation, Utilities):
+    class Measurements(InitializeSimulation):
         def __init__(self,
                     *args,
                     **kwargs):
@@ -206,7 +208,9 @@ The *Measurements* class inherits both the *InitializeSimulation* and
 Finally, let us create the three remaining classes, named *MinimizeEnergy*,
 *MonteCarlo*, and *MolecularDynamics*. Each of these three classes inherits
 from the *Measurements* class, and thus from the classes inherited by
-*Measurements*. Within the *MinimizeEnergy.py* file, copy the following lines:
+*Measurements*.
+
+Within the *MinimizeEnergy.py* file, copy the following lines:
 
 .. label:: start_MinimizeEnergy_class
 
@@ -317,7 +321,7 @@ any *AssertionError*:
     Utilities does not inherit from MonteCarlo, as expected
     MonteCarlo correctly inherits from Utilities
 
-Alternatively, this test can also be launched using Pytest by typing in a terminal:
+Alternatively, this test can also be launched using *Pytest* by typing in a terminal:
 
 .. code-block:: bash
 

docs/source/chapters/chapter2.rst

@@ -95,13 +95,15 @@ Modify the *Prepare* class as follows:
 
     class Prepare:
         def __init__(self,
-                    number_atoms=[10],  # List - no unit
-                    epsilon=[0.1],  # List - Kcal/mol
-                    sigma=[3],  # List - Angstrom
-                    atom_mass=[10],  # List - g/mol
+                    ureg, # Pint unit registry
+                    number_atoms, # List - no unit
+                    epsilon, # List - Kcal/mol
+                    sigma, # List - Angstrom
+                    atom_mass,  # List - g/mol
                     potential_type="Lennard-Jones",
                     *args,
                     **kwargs):
+            self.ureg = ureg
             self.number_atoms = number_atoms
             self.epsilon = epsilon
             self.sigma = sigma
@@ -116,7 +118,7 @@ Here, the four lists *number_atoms* :math:`N`, *epsilon* :math:`\epsilon`,
 :math:`10`, :math:`0.1~\text{[Kcal/mol]}`, :math:`3~\text{[Å]}`, and
 :math:`10~\text{[g/mol]}`, respectively.
 
-The type of potential is also specified, with Lennard-Jones being closen as
+The type of potential is also specified, with Lennard-Jones being chosen as
 the default option.
 
 All the parameters are assigned to *self*, allowing other methods to access
@@ -136,25 +138,33 @@ unit system to the *LJ* unit system:
 
     def calculate_LJunits_prefactors(self):
         """Calculate the Lennard-Jones units prefactors."""
-        # Define the reference distance, energy, and mass
-        self.reference_distance = self.sigma[0]  # Angstrom
-        self.reference_energy = self.epsilon[0]  # kcal/mol
-        self.reference_mass = self.atom_mass[0]  # g/mol
-
-        # Calculate the prefactor for the time
-        mass_kg = self.atom_mass[0]/cst.kilo/cst.Avogadro  # kg
-        epsilon_J = self.epsilon[0]*cst.calorie*cst.kilo/cst.Avogadro  # J
-        sigma_m = self.sigma[0]*cst.angstrom  # m
-        time_s = np.sqrt(mass_kg*sigma_m**2/epsilon_J)  # s
-        self.reference_time = time_s / cst.femto  # fs
-
-        # Calculate the prefactor for the temperature
+        # First define constants
         kB = cst.Boltzmann*cst.Avogadro/cst.calorie/cst.kilo  # kcal/mol/K
-        self.reference_temperature = self.epsilon[0]/kB  # K
-
-        # Calculate the prefactor for the pressure
-        pressure_pa = epsilon_J/sigma_m**3  # Pa
-        self.reference_pressure = pressure_pa/cst.atm  # atm
+        kB *= self.ureg.kcal/self.ureg.mol/self.ureg.kelvin
+        Na = cst.Avogadro/self.ureg.mol
+        # Define the reference distance, energy, and mass
+        self.ref_length = self.sigma[0]  # Angstrom
+        self.ref_energy = self.epsilon[0]  # kcal/mol
+        self.ref_mass = self.atom_mass[0]  # g/mol
+        # Optional: assert that units were correctly provided by users
+        assert self.ref_length.units == self.ureg.angstrom, \
+            f"Error: Provided sigma has wrong units, should be angstrom"
+        assert self.ref_energy.units == self.ureg.kcal/self.ureg.mol, \
+            f"Error: Provided epsilon has wrong units, should be kcal/mol"
+        assert self.ref_mass.units == self.ureg.g/self.ureg.mol, \
+            f"Error: Provided mass has wrong units, should be g/mol"
+        # Calculate the prefactor for the time (in femtosecond)
+        self.ref_time = np.sqrt(self.ref_mass \
+            *self.ref_length**2/self.ref_energy).to(self.ureg.femtosecond)
+        # Calculate the prefactor for the temperature (in Kelvin)
+        self.ref_temperature = self.ref_energy/kB  # Kelvin
+        # Calculate the prefactor for the pressure (in Atmosphere)
+        self.ref_pressure = (self.ref_energy \
+            /self.ref_length**3/Na).to(self.ureg.atmosphere)
+        # Regroup all the reference quantities in list, for practicality
+        self.ref_quantities = [self.ref_length, self.ref_energy,
+            self.ref_mass, self.ref_time, self.ref_pressure]
+        self.ref_units = [ref.units for ref in self.ref_quantities]
 
 .. label:: end_Prepare_class
 
@@ -192,23 +202,26 @@ Let us take advantage of the calculated reference values and normalize the
 three inputs of the *Prepare* class that have physical dimensions, i.e.,
 *epsilon*, *sigma*, and *atom_mass*.
 
-Create a new method called *nondimensionalize_units_0* within the *Prepare*
+Create a new method called *nondimensionalize_units* within the *Prepare*
 class:
 
 .. label:: start_Prepare_class
 
 .. code-block:: python
 
-    def nondimensionalize_units_0(self):
-        # Normalize LJ properties
-        epsilon, sigma, atom_mass = [], [], []
-        for e0, s0, m0 in zip(self.epsilon, self.sigma, self.atom_mass):
-            epsilon.append(e0/self.reference_energy)
-            sigma.append(s0/self.reference_distance)
-            atom_mass.append(m0/self.reference_mass)
-        self.epsilon = epsilon
-        self.sigma = sigma
-        self.atom_mass = atom_mass
+    def nondimensionalize_units(self, quantities_to_normalise):
+        for quantity in quantities_to_normalise:
+            if isinstance(quantity, list):
+                for i, element in enumerate(quantity):
+                    assert element.units in self.ref_units, \
+                        f"Error: Units not part of the reference units"
+                    ref_value = self.ref_quantities[self.ref_units.index(element.units)]
+                    quantity[i] = element/ref_value
+                    assert quantity[i].units == self.ureg.dimensionless, \
+                        f"Error: Quantities are not properly nondimensionalized"
+                    quantity[i] = quantity[i].magnitude # get rid of ureg
+            else:
+                print("NON ANTICIPATED!")
 
 .. label:: end_Prepare_class
 
@@ -218,7 +231,7 @@ will be used to nondimensionalize units in future classes. We anticipate that
 each element is normalized with the corresponding reference value. The
 *zip()* function allows us to loop over all three lists at once.
 
-Let us also call the *nondimensionalize_units_0* from the *__init__()* method
+Let us also call the *nondimensionalize_units* from the *__init__()* method
 of the *Prepare* class:
 
 .. label:: start_Prepare_class
@@ -228,7 +241,7 @@ of the *Prepare* class:
     def __init__(self,
         (...)
         self.calculate_LJunits_prefactors()
-        self.nondimensionalize_units_0()
+        self.nondimensionalize_units([self.epsilon, self.sigma, self.atom_mass])
 
 .. label:: end_Prepare_class
 
@@ -288,7 +301,7 @@ Let us call the *identify_atom_properties* from the *__init__()* method:
 
     def __init__(self,
         (...)
-        self.nondimensionalize_units_0()
+        self.nondimensionalize_units([self.epsilon, self.sigma, self.atom_mass])
         self.identify_atom_properties()
 
 .. label:: end_Prepare_class
@@ -306,13 +319,25 @@ type 1, and 3 atoms of type 2:
 
     import numpy as np
     from Prepare import Prepare
-
-    # Initialize the Prepare object
+    from pint import UnitRegistry
+    ureg = UnitRegistry()
+
+    # Define atom number of each group
+    nmb_1, nmb_2= [2, 3]
+    # Define LJ parameters (sigma)
+    sig_1, sig_2 = [3, 4]*ureg.angstrom
+    # Define LJ parameters (epsilon)
+    eps_1, eps_2 = [0.2, 0.4]*ureg.kcal/ureg.mol
+    # Define atom mass
+    mss_1, mss_2 = [10, 20]*ureg.gram/ureg.mol
+
+    # Initialize the prepare object
     prep = Prepare(
-        number_atoms=[2, 3],
-        epsilon=[0.2, 0.4], # kcal/mol
-        sigma=[3, 4], # A
-        atom_mass=[10, 20], # g/mol
+        ureg = ureg,
+        number_atoms=[nmb_1, nmb_2],
+        epsilon=[eps_1, eps_2], # kcal/mol
+        sigma=[sig_1, sig_2], # A
+        atom_mass=[mss_1, mss_2], # g/mol
     )
 
     # Test function using pytest