From 717f020d7fc2a42ebecd7ba91fea3767e3d8c2d6 Mon Sep 17 00:00:00 2001
From: Simon Gravelle <simon.gravelle@live.fr>
Date: Sat, 31 Aug 2024 18:40:11 +0200
Subject: [PATCH 01/13] improved units managments

---
 docs/source/chapters/chapter1.rst |  78 ++++++-------
 docs/source/chapters/chapter2.rst | 111 +++++++++++--------
 docs/source/chapters/chapter3.rst | 128 +++++++++++++++++++--
 docs/source/chapters/chapter4.rst |  90 +--------------
 docs/source/chapters/chapter6.rst |   9 +-
 docs/source/chapters/chapter8.rst | 178 +++++++++++++++++++++---------
 docs/source/index.rst             |   1 +
 tests/build-documentation.py      |   2 +-
 8 files changed, 360 insertions(+), 237 deletions(-)

diff --git a/docs/source/chapters/chapter1.rst b/docs/source/chapters/chapter1.rst
index 92bbed6..3fd1b8a 100644
--- a/docs/source/chapters/chapter1.rst
+++ b/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
 
diff --git a/docs/source/chapters/chapter2.rst b/docs/source/chapters/chapter2.rst
index 5752c44..842c0ff 100644
--- a/docs/source/chapters/chapter2.rst
+++ b/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
diff --git a/docs/source/chapters/chapter3.rst b/docs/source/chapters/chapter3.rst
index 4731e76..b67141b 100644
--- a/docs/source/chapters/chapter3.rst
+++ b/docs/source/chapters/chapter3.rst
@@ -28,8 +28,10 @@ Let us improve the previously created *InitializeSimulation* class:
                     box_dimensions=[10, 10, 10],  # List - Angstroms
                     seed=None,  # Int
                     initial_positions=None,  # Array - Angstroms
-                    thermo_period=None,
-                    dumping_period=None,
+                    #thermo_period=None, # TOFIX to be added later
+                    #dumping_period=None, # TOFIX to be added later
+                    neighbor=1,
+                    cut_off=10,
                     *args,
                     **kwargs,
                     ):
@@ -38,9 +40,12 @@ Let us improve the previously created *InitializeSimulation* class:
             # If a box dimension was entered as 0, dimensions will be 2
             self.dimensions = len(list(filter(lambda x: x > 0, box_dimensions)))
             self.seed = seed
+            self.neighbor = neighbor
+            self.cut_off = cut_off
+            self.step = 0 # initialize simulation step
             self.initial_positions = initial_positions
-            self.thermo_period = thermo_period
-            self.dumping_period = dumping_period
+            #self.thermo_period = thermo_period # TOFIX to be added later
+            #self.dumping_period = dumping_period # TOFIX to be added later
 
 .. label:: end_InitializeSimulation_class
 
@@ -57,9 +62,16 @@ Several parameters are provided to the *InitializeSimulation* class:
   of length corresponding to the number of atoms. If *initial_positions* 
   is left equal to *None*, positions will be randomly assigned to the
   atoms (see below).
-- Optionally, a thermo period and a dumping_period can be provided to
-  control the outputs from the simulation (it will be implemented
-  in :ref:`chapter5-label`).
+- A cut off value with a default of 10 Ångströms is defined.
+- The *neighbor* parameter determines the interval between recalculations of
+  the neighbor lists. By default, a value of 1 is used, meaning that neighbor
+  list will be reconstructed every steps. In certain cases, the number will be
+  increased to recude the computation time.
+
+..
+    - Optionally, a thermo period and a dumping_period can be provided to  # TOFIX to be added later
+    control the outputs from the simulation (it will be implemented  # TOFIX to be added later
+    in :ref:`chapter5-label`).  # TOFIX to be added later
 
 Finally, the *dimensions* of the system are calculated as the length of
 *box_dimensions*.
@@ -110,6 +122,7 @@ parameters, here the *box_dimensions* as well as the *initial_positions*:
         # Normalize the box dimensions
         if self.initial_positions is not None:
             self.initial_positions = self.initial_positions/self.reference_distance
+        self.cut_off /= self.reference_distance
 
 .. label:: end_InitializeSimulation_class
 
@@ -219,6 +232,104 @@ Let us call the *populate_box* method from the *__init__* class:
         self.populate_box()
 
 .. label:: end_InitializeSimulation_class
+
+Build neighbor lists
+--------------------
+
+In molecular simulations, it is common practice to identify neighboring atoms
+to save computational time. By focusing only on interactions between
+neighboring atoms, the simulation becomes more efficient. Add the following
+*update_neighbor_lists()* method to the *Utilities* class:
+
+.. label:: start_Utilities_class
+
+.. code-block:: python
+
+    def update_neighbor_lists(self):
+        if (self.step % self.neighbor == 0):
+            matrix = distances.contact_matrix(self.atoms_positions,
+                cutoff=self.cut_off, #+2,
+                returntype="numpy",
+                box=self.box_size)
+            neighbor_lists = []
+            for cpt, array in enumerate(matrix[:-1]):
+                list = np.where(array)[0].tolist()
+                list = [ele for ele in list if ele > cpt]
+                neighbor_lists.append(list)
+            self.neighbor_lists = neighbor_lists
+
+.. label:: end_Utilities_class
+
+The *update_neighbor_lists()* method generates neighbor lists for each
+atom, ensuring that only relevant interactions are considered in the
+calculations. These lists will be recalculated at intervals specified by
+the *neighbor* input parameter.
+
+Update cross coefficients
+-------------------------
+
+At the same time as the neighbor lists are getting build up, let us also
+pre-calculate the cross coefficients. This will make the force calculation
+more practical (see below).
+
+.. label:: start_Utilities_class
+
+.. code-block:: python
+
+    def update_cross_coefficients(self):
+        if (self.step % self.neighbor == 0):
+            # Precalculte LJ cross-coefficients
+            sigma_ij_list = []
+            epsilon_ij_list = []
+            for Ni in np.arange(self.total_number_atoms-1): # tofix error for GCMC
+                # Read information about atom i
+                sigma_i = self.atoms_sigma[Ni]
+                epsilon_i = self.atoms_epsilon[Ni]
+                neighbor_of_i = self.neighbor_lists[Ni]
+                # Read information about neighbors j
+                sigma_j = self.atoms_sigma[neighbor_of_i]
+                epsilon_j = self.atoms_epsilon[neighbor_of_i]
+                # Calculare cross parameters
+                sigma_ij_list.append((sigma_i+sigma_j)/2)
+                epsilon_ij_list.append((epsilon_i+epsilon_j)/2)
+            self.sigma_ij_list = sigma_ij_list
+            self.epsilon_ij_list = epsilon_ij_list
+
+.. label:: end_Utilities_class
+
+Here, the values of the cross coefficients between atom of type 1 and 2,
+:math:`\sigma_{12}` and :math:`\epsilon_{12}`, are assumed to follow the arithmetic mean:
+
+.. math::
+
+    \sigma_{12} = (\sigma_{11}+\sigma_{22})/2 \\
+    \epsilon_{12} = (\epsilon_{11}+\epsilon_{22})/2
+
+Finally, import the following library in the *Utilities.py* file:
+
+.. label:: start_Utilities_class
+
+.. code-block:: python
+
+    import numpy as np
+    from MDAnalysis.analysis import distances
+
+.. label:: end_Utilities_class
+
+Let us call the *update_neighbor_lists* and *update_cross_coefficients*
+methods from the *__init__* class:
+
+.. label:: start_InitializeSimulation_class
+
+.. code-block:: python
+
+    def __init__(self,
+        (...)
+        self.populate_box()
+        self.update_neighbor_lists()
+        self.update_cross_coefficients()
+
+.. label:: end_InitializeSimulation_class
         
 Test the code
 -------------
@@ -249,7 +360,8 @@ boundaries along all 3 dimensions of space:
         atoms_positions = init.atoms_positions
         for atom_position in atoms_positions:
             for x, boundary in zip(atom_position, box_boundaries):
-                assert (x >= boundary[0]) and (x <= boundary[1]), f"Test failed: Atoms outside of the box at position {atom_position}"
+                assert (x >= boundary[0]) and (x <= boundary[1]), \
+                f"Test failed: Atoms outside of the box at position {atom_position}"
         print("Test passed")
 
     # If the script is run directly, execute the tests
diff --git a/docs/source/chapters/chapter4.rst b/docs/source/chapters/chapter4.rst
index e460b16..63c3b1a 100644
--- a/docs/source/chapters/chapter4.rst
+++ b/docs/source/chapters/chapter4.rst
@@ -75,8 +75,6 @@ Then, let us fill the *__init__()* method:
                     data_folder=None,
                     *args,
                     **kwargs):
-            self.neighbor = neighbor
-            self.cut_off = cut_off
             self.displacement = displacement
             self.maximum_steps = maximum_steps
             self.thermo_outputs = thermo_outputs
@@ -89,9 +87,7 @@ Then, let us fill the *__init__()* method:
 .. label:: end_MinimizeEnergy_class
 
 An important parameter is *maximum_steps*, which sets the maximum number
-of steps for the energy minimization process. A *cut_off* value with a
-default of 9 Ångströms is also defined. The *neighbor* parameter determines
-the interval between recalculations of the neighbor lists, and the *displacement*
+of steps for the energy minimization process. The *displacement*
 parameter, with a default value of 0.01 Ångström, sets the initial atom
 displacement value.
 
@@ -111,7 +107,6 @@ method to the *MinimizeEnergy* class:
 
     def nondimensionalize_units_2(self):
         """Use LJ prefactors to convert units into non-dimensional."""
-        self.cut_off = self.cut_off/self.reference_distance
         self.displacement = self.displacement/self.reference_distance
 
 .. label:: end_MinimizeEnergy_class
@@ -175,89 +170,6 @@ The displacement, which has an initial value of 0.01, is adjusted through energy
 minimization. When the trial is successful, its value is multiplied by 1.2. When
 the trial is rejected, its value is multiplied by 0.2.
 
-Build neighbor lists
---------------------
-
-In molecular simulations, it is common practice to identify neighboring atoms
-to save computational time. By focusing only on interactions between
-neighboring atoms, the simulation becomes more efficient. Add the following
-*update_neighbor_lists()* method to the *Utilities* class:
-
-.. label:: start_Utilities_class
-
-.. code-block:: python
-
-    def update_neighbor_lists(self):
-        if (self.step % self.neighbor == 0):
-            matrix = distances.contact_matrix(self.atoms_positions,
-                cutoff=self.cut_off, #+2,
-                returntype="numpy",
-                box=self.box_size)
-            neighbor_lists = []
-            for cpt, array in enumerate(matrix[:-1]):
-                list = np.where(array)[0].tolist()
-                list = [ele for ele in list if ele > cpt]
-                neighbor_lists.append(list)
-            self.neighbor_lists = neighbor_lists
-
-.. label:: end_Utilities_class
-
-The *update_neighbor_lists()* method generates neighbor lists for each
-atom, ensuring that only relevant interactions are considered in the
-calculations. These lists will be recalculated at intervals specified by
-the *neighbor* input parameter.
-
-Update cross coefficients
--------------------------
-
-At the same time as the neighbor lists are getting build up, let us also
-pre-calculate the cross coefficients. This will make the force calculation
-more practical (see below).
-
-.. label:: start_Utilities_class
-
-.. code-block:: python
-
-    def update_cross_coefficients(self):
-        if (self.step % self.neighbor == 0):
-            # Precalculte LJ cross-coefficients
-            sigma_ij_list = []
-            epsilon_ij_list = []
-            for Ni in np.arange(self.total_number_atoms-1): # tofix error for GCMC
-                # Read information about atom i
-                sigma_i = self.atoms_sigma[Ni]
-                epsilon_i = self.atoms_epsilon[Ni]
-                neighbor_of_i = self.neighbor_lists[Ni]
-                # Read information about neighbors j
-                sigma_j = self.atoms_sigma[neighbor_of_i]
-                epsilon_j = self.atoms_epsilon[neighbor_of_i]
-                # Calculare cross parameters
-                sigma_ij_list.append((sigma_i+sigma_j)/2)
-                epsilon_ij_list.append((epsilon_i+epsilon_j)/2)
-            self.sigma_ij_list = sigma_ij_list
-            self.epsilon_ij_list = epsilon_ij_list
-
-.. label:: end_Utilities_class
-
-Here, the values of the cross coefficients between atom of type 1 and 2,
-:math:`\sigma_{12}` and :math:`\epsilon_{12}`, are assumed to follow the arithmetic mean:
-
-.. math::
-
-    \sigma_{12} = (\sigma_{11}+\sigma_{22})/2 \\
-    \epsilon_{12} = (\epsilon_{11}+\epsilon_{22})/2
-
-Finally, import the following library in the *Utilities.py* file:
-
-.. label:: start_Utilities_class
-
-.. code-block:: python
-
-    import numpy as np
-    from MDAnalysis.analysis import distances
-
-.. label:: end_Utilities_class
-
 Compute_potential
 -----------------
 
diff --git a/docs/source/chapters/chapter6.rst b/docs/source/chapters/chapter6.rst
index 52fe593..7daaaf1 100644
--- a/docs/source/chapters/chapter6.rst
+++ b/docs/source/chapters/chapter6.rst
@@ -47,6 +47,10 @@ Let us add a method named *monte_carlo_move* to the *MonteCarlo* class:
     def monte_carlo_move(self):
         """Monte Carlo move trial."""
         if self.displace_mc is not None: # only trigger if displace_mc was provided by the user
+            # If needed, recalculate neighbor/coeff lists
+            self.update_neighbor_lists()
+            self.identify_atom_properties()
+            self.update_cross_coefficients()
             try: # try using the last saved Epot, if it exists
                 initial_Epot = self.Epot
             except: # If self.Epot does not exists yet, calculate it
@@ -146,16 +150,13 @@ perform a loop over the desired number of steps *maximum_steps*:
     def run(self):
         """Perform the loop over time."""
         for self.step in range(0, self.maximum_steps+1):
-            self.update_neighbor_lists()
-            self.update_cross_coefficients()
             self.monte_carlo_move()
             self.wrap_in_box()
 
 .. label:: end_MonteCarlo_class
 
 At each step, the *monte_carlo_move* method is called. The previously defined
-methods *update_neighbor_lists* and *wrap_in_box* are also called to ensure that
-the neighbor lists are kept up to date despite the motion of the atoms, and that
+mthe *wrap_in_box* method is also called to ensure that
 the atoms remain inside the box, respectively.
 
 Let us call *update_log_md_mc* from the run method of the MonteCarlo class.
diff --git a/docs/source/chapters/chapter8.rst b/docs/source/chapters/chapter8.rst
index ba2926b..9982fc4 100644
--- a/docs/source/chapters/chapter8.rst
+++ b/docs/source/chapters/chapter8.rst
@@ -3,64 +3,131 @@
 Monte Carlo insert
 ==================
 
-In the Grand Canonical ensemble, the number of atom in the
-simulation box is not constant. 
+Here, a *Monte Carlo* simulation is implemented in the Grand Canonical ensemble,
+where the number of atoms in the system fluctuates. The principle of the
+simulation resemble the Monte Carlo move from :ref:`chapter6-label`:
+
+- 1) We start from a given intial configuration, and measure the potential
+  energy, :math:`E_\text{pot}^\text{initial}`.
+- 2) A random number is chosen, and depending on its value, either a new particle
+  will tried to be inserted, or an existing particle will tried to be deleted.
+- 3) The energy of the system after the insersion/deletion,
+  :math:`E_\text{pot}^\text{trial}`, is measured.
+- 4) We then decide to keep or reject the move by calculating
+  the difference in energy between the trial and the initial configurations:
+  :math:`\Delta E = E_\text{pot}^\text{trial} - E_\text{pot}^\text{initial}`.
+  
+  - If :math:`\Delta E < 0`, then the move is automatically accepted. 
+  - If :math:`\Delta E > 0`, then the move is accepted with a probability given
+    by the Boltzmann factor :math:`\exp{- \beta \Delta E}`, where
+    :math:`\beta = 1 / k_\text{B} T` and :math:`T` is the imposed temperature.
+
+- 5) Steps 1-4 are repeated a large number of times, generating a broad range of
+     possible configurations.
+
+Implementation
+--------------
+
+Let us add the following method called *monte_carlo_exchange* to the *MonteCarlo* class:
 
-Let us add the following method to the MonteCarlo class:
+.. label:: start_MonteCarlo_class
+
+.. code-block:: python
+
+    def monte_carlo_exchange(self):
+        # The first step is to make a copy of the previous state
+        # Since atoms numbers are evolving, its also important to store the
+        # neighbor, sigma, and epsilon lists
+        self.Epot = self.compute_potential()
+        initial_positions = copy.deepcopy(self.atoms_positions)
+        initial_number_atoms = copy.deepcopy(self.number_atoms)
+        initial_neighbor_lists = copy.deepcopy(self.neighbor_lists)
+        initial_sigma_lists = copy.deepcopy(self.sigma_ij_list)
+        initial_epsilon_lists = copy.deepcopy(self.epsilon_ij_list)
+        # Apply a 50-50 probability to insert or delete
+        if np.random.random() < 0.5:
+            self.monte_carlo_insert()
+        else:
+            self.monte_carlo_delete()
+        # If np.random.random() < self.acceptation_probability, do nothing
+        if np.random.random() > self.acceptation_probability:
+            # Reject the new position, revert to inital position
+            self.neighbor_lists = initial_neighbor_lists
+            self.sigma_ij_list = initial_sigma_lists
+            self.epsilon_ij_list = initial_epsilon_lists
+            self.atoms_positions = initial_positions
+            self.number_atoms = initial_number_atoms
+
+.. label:: end_MonteCarlo_class
+
+First, the potential energy is calculated, and the current state of the
+simulations, such as positions and atom numbers, are saved. Then, an insertion
+try or a deletion trial is made, each with a probability of 0.5. The
+*monte_carlo_insert* and *monte_carlo_delete*, that are implemented here below,
+are both calculting the *acceptation_probability*. A random number is again selected,
+and if that number is larger than the acceptation probability, then the system
+revert to its previous position. If the random number is lower than the acceptation
+probability, nothing appends, which means that the last trial is accepted.
+
+Then, let us write the *monte_carlo_insert* method:
+
+.. label:: start_MonteCarlo_class
+
+.. code-block:: python
+
+    def monte_carlo_insert(self):
+        self.number_atoms[self.inserted_type] += 1
+        new_atom_position = np.random.random(3)*np.diff(self.box_boundaries).T \
+            - np.diff(self.box_boundaries).T/2
+        shift_id = 0 
+        for N in self.number_atoms[:self.inserted_type]:
+            shift_id += N
+        self.atoms_positions = np.vstack([self.atoms_positions[:shift_id],
+                                        new_atom_position,
+                                        self.atoms_positions[shift_id:]])
+        self.total_number_atoms = np.sum(self.number_atoms)
+        self.update_neighbor_lists()
+        self.identify_atom_properties()
+        self.update_cross_coefficients()
+        trial_Epot = self.compute_potential()
+        Lambda = self.calculate_Lambda(self.atom_mass[self.inserted_type])
+        beta =  1/self.desired_temperature
+        Nat = np.sum(self.number_atoms) # NUmber atoms, should it relly be N? of N (type) ?
+        Vol = np.prod(self.box_size[:3]) # box volume
+        # dimension of 3 is enforced in the power of the Lambda
+        self.acceptation_probability = np.min([1, Vol/(Lambda**3*Nat) \
+            *np.exp(beta*(self.desired_mu-trial_Epot+self.Epot))])
+
+.. label:: end_MonteCarlo_class
+
+as well as the *monte_carlo_delete* method:
 
 .. label:: start_MonteCarlo_class
 
 .. code-block:: python
 
-    def monte_carlo_insert_delete(self):
-        if self.desired_mu is not None: # trigger method if desired_mu was entered
-            initial_Epot = self.compute_potential(output="potential")
-            initial_positions = copy.deepcopy(self.atoms_positions)
-            initial_number_atoms = copy.deepcopy(self.number_atoms)
-            initial_neighbor_lists = copy.deepcopy(self.neighbor_lists)
-            volume = np.prod(np.diff(self.box_boundaries))
-            if np.random.random() < 0.5:
-                # Try adding an atom
-                self.number_atoms[self.inserted_type] += 1
-                atom_position = np.zeros((1, self.dimensions))
-                for dim in np.arange(self.dimensions):
-                    atom_position[:, dim] = np.random.random(1)*np.diff(self.box_boundaries[dim]) - np.diff(self.box_boundaries[dim])/2
-                shift_id = 0
-                for N in self.number_atoms[:self.inserted_type]:
-                    shift_id += N
-                self.atoms_positions = np.vstack([self.atoms_positions[:shift_id], atom_position, self.atoms_positions[shift_id:]])
-                self.update_neighbor_lists()
-                self.calculate_cross_coefficients()
-                trial_Epot = self.compute_potential(output="potential")
-                Lambda = self.calculate_Lambda(self.atom_mass[self.inserted_type])
-                beta =  1/self.desired_temperature
-                acceptation_probability = np.min([1, volume/(Lambda**self.dimensions*(self.total_number_atoms))*np.exp(beta*(self.desired_mu-trial_Epot+initial_Epot))])
-            else:
-                # Pick one atom to delete randomly
-                atom_id = np.random.randint(self.number_atoms[self.inserted_type])
-                self.number_atoms[self.inserted_type] -= 1
-                if self.number_atoms[self.inserted_type] > 0:
-                    shift_id = 0
-                    for N in self.number_atoms[:self.inserted_type]:
-                        shift_id += N
-                    self.atoms_positions = np.delete(self.atoms_positions, shift_id+atom_id, axis=0)
-                    self.update_neighbor_lists()
-                    self.calculate_cross_coefficients()
-                    trial_Epot = self.compute_potential(output="potential")
-                    Lambda = self.calculate_Lambda(self.atom_mass[self.inserted_type])
-                    beta =  1/self.desired_temperature
-                    acceptation_probability = np.min([1, (Lambda**self.dimensions*(self.total_number_atoms-1)/volume)*np.exp(-beta*(self.desired_mu+trial_Epot-initial_Epot))])
-                else:
-                    acceptation_probability = 0
-            if np.random.random() < acceptation_probability:
-                # Accept the new position
-                pass
-            else:
-                # Reject the new position
-                self.neighbor_lists = initial_neighbor_lists
-                self.atoms_positions = initial_positions
-                self.number_atoms = initial_number_atoms
-                self.calculate_cross_coefficients()
+    def monte_carlo_delete(self):
+        # Pick one atom to delete randomly
+        atom_id = np.random.randint(self.number_atoms[self.inserted_type])
+        self.number_atoms[self.inserted_type] -= 1
+        if self.number_atoms[self.inserted_type] > 0:
+            shift_id = 0
+            for N in self.number_atoms[:self.inserted_type]:
+                shift_id += N
+            self.atoms_positions = np.delete(self.atoms_positions, shift_id+atom_id, axis=0)
+            self.update_neighbor_lists()
+            self.identify_atom_properties()
+            self.update_cross_coefficients()
+            trial_Epot = self.compute_potential()
+            Lambda = self.calculate_Lambda(self.atom_mass[self.inserted_type])
+            beta =  1/self.desired_temperature
+            Nat = np.sum(self.number_atoms) # NUmber atoms, should it relly be N? of N (type) ?
+            Vol = np.prod(self.box_size[:3]) # box volume
+            # dimension of 3 is enforced in the power of the Lambda
+            self.acceptation_probability = np.min([1, (Lambda**3 *(Nat-1)/Vol) \
+                *np.exp(-beta*(self.desired_mu+trial_Epot-self.Epot))])
+        else:
+            print("Error: no more atoms to delete")
 
 .. label:: end_MonteCarlo_class
 
@@ -113,7 +180,7 @@ Let us non-dimentionalize desired_mu by adding:
 
 .. label:: end_MonteCarlo_class
 
-Finally, *monte_carlo_insert_delete* must be included in the run:
+Finally, the *monte_carlo_insert_delete()* method must be included in the run:
 
 .. label:: start_MonteCarlo_class
 
@@ -121,8 +188,8 @@ Finally, *monte_carlo_insert_delete* must be included in the run:
 
     def run(self):
         (...)
-            self.monte_carlo_displacement()
-            self.monte_carlo_insert_delete()
+            self.monte_carlo_move()
+            self.monte_carlo_exchange()
             self.wrap_in_box()
         (...)
 
@@ -172,6 +239,7 @@ One can use a similar test as previously. Let us use a displace distance of
         sigma=[3], # A
         atom_mass=[1], # g/mol
         box_dimensions=[20, 20, 20], # A
+        data_folder = "outputs/",
         )
     mc.run()
 
diff --git a/docs/source/index.rst b/docs/source/index.rst
index cf64b19..f094c45 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -57,6 +57,7 @@ Learn Molecular Simulations with Python can be followed by anyone with a compute
     chapters/chapter5
     chapters/chapter6
     chapters/chapter7
+    chapters/chapter8
 
 .. toctree::
     :maxdepth: 2
diff --git a/tests/build-documentation.py b/tests/build-documentation.py
index 247c50a..fbaa9bc 100644
--- a/tests/build-documentation.py
+++ b/tests/build-documentation.py
@@ -18,7 +18,7 @@
     os.mkdir("generated-codes/")
 
 # loop on the different chapter
-for chapter_id in [1, 2, 3, 4, 5, 6, 7]:
+for chapter_id in [1, 2]:
     # for each chapter, create the corresponding code
     RST_EXISTS, created_tests, folder = sphinx_to_python(path_to_docs, chapter_id)
     if RST_EXISTS:

From ff5a5983b44f2180fceb51a22f3f9dc261c8054f Mon Sep 17 00:00:00 2001
From: Simon Gravelle <simon.gravelle@live.fr>
Date: Sat, 31 Aug 2024 21:00:32 +0200
Subject: [PATCH 02/13] improved units normalizing protocole

---
 docs/source/chapters/chapter2.rst |  29 +++++---
 docs/source/chapters/chapter3.rst | 115 ++++++++++--------------------
 tests/build-documentation.py      |   2 +-
 3 files changed, 60 insertions(+), 86 deletions(-)

diff --git a/docs/source/chapters/chapter2.rst b/docs/source/chapters/chapter2.rst
index 842c0ff..468edb4 100644
--- a/docs/source/chapters/chapter2.rst
+++ b/docs/source/chapters/chapter2.rst
@@ -210,7 +210,8 @@ class:
 .. code-block:: python
 
     def nondimensionalize_units(self, quantities_to_normalise):
-        for quantity in quantities_to_normalise:
+        for name in quantities_to_normalise:
+            quantity = getattr(self, name)  # Get the attribute by name
             if isinstance(quantity, list):
                 for i, element in enumerate(quantity):
                     assert element.units in self.ref_units, \
@@ -220,8 +221,19 @@ class:
                     assert quantity[i].units == self.ureg.dimensionless, \
                         f"Error: Quantities are not properly nondimensionalized"
                     quantity[i] = quantity[i].magnitude # get rid of ureg
+                setattr(self, name, quantity)
             else:
-                print("NON ANTICIPATED!")
+                if quantity is not None:
+                    assert np.shape(quantity) == (), \
+                        f"Error: The quantity is a list or an array"
+                    assert quantity.units in self.ref_units, \
+                        f"Error: Units not part of the reference units"
+                    ref_value = self.ref_quantities[self.ref_units.index(quantity.units)]
+                    quantity = quantity/ref_value
+                    assert quantity.units == self.ureg.dimensionless, \
+                        f"Error: Quantities are not properly nondimensionalized"
+                    quantity = quantity.magnitude # get rid of ureg
+                    setattr(self, name, quantity)
 
 .. label:: end_Prepare_class
 
@@ -241,7 +253,7 @@ of the *Prepare* class:
     def __init__(self,
         (...)
         self.calculate_LJunits_prefactors()
-        self.nondimensionalize_units([self.epsilon, self.sigma, self.atom_mass])
+        self.nondimensionalize_units(["epsilon", "sigma", "atom_mass"])
 
 .. label:: end_Prepare_class
 
@@ -270,8 +282,7 @@ within the *Prepare* class:
 .. code-block:: python
 
     def identify_atom_properties(self):
-        """Identify the atom properties for each atom."""
-        self.total_number_atoms = np.sum(self.number_atoms)
+        """Identify the properties for each atom."""
         atoms_sigma = []
         atoms_epsilon = []
         atoms_mass = []
@@ -301,7 +312,7 @@ Let us call the *identify_atom_properties* from the *__init__()* method:
 
     def __init__(self,
         (...)
-        self.nondimensionalize_units([self.epsilon, self.sigma, self.atom_mass])
+        self.nondimensionalize_units(["epsilon", "sigma", "atom_mass"])
         self.identify_atom_properties()
 
 .. label:: end_Prepare_class
@@ -309,9 +320,9 @@ Let us call the *identify_atom_properties* from the *__init__()* method:
 Test the code
 -------------
 
-Let us test the *Prepare* class to make sure that it does what is expected.
-Just like in the previous example, let us initiates a system with 2 atoms of
-type 1, and 3 atoms of type 2:
+Let's test the *Prepare* class to make sure that it does what is expected.
+Here, a system containing 2 atoms of type 1, and 3 atoms of type 2 is
+prepared. LJs parameters and masses for each groups are also defined.
 
 .. label:: start_test_2a_class
 
diff --git a/docs/source/chapters/chapter3.rst b/docs/source/chapters/chapter3.rst
index b67141b..1fe44b2 100644
--- a/docs/source/chapters/chapter3.rst
+++ b/docs/source/chapters/chapter3.rst
@@ -25,27 +25,19 @@ Let us improve the previously created *InitializeSimulation* class:
 
     class InitializeSimulation(Prepare):
         def __init__(self,
-                    box_dimensions=[10, 10, 10],  # List - Angstroms
-                    seed=None,  # Int
+                    box_dimensions,  # List - Angstroms
+                    cut_off, # Angstroms
                     initial_positions=None,  # Array - Angstroms
-                    #thermo_period=None, # TOFIX to be added later
-                    #dumping_period=None, # TOFIX to be added later
-                    neighbor=1,
-                    cut_off=10,
+                    neighbor=1, # Integer
                     *args,
                     **kwargs,
                     ):
             super().__init__(*args, **kwargs)
             self.box_dimensions = box_dimensions
-            # If a box dimension was entered as 0, dimensions will be 2
-            self.dimensions = len(list(filter(lambda x: x > 0, box_dimensions)))
-            self.seed = seed
-            self.neighbor = neighbor
             self.cut_off = cut_off
+            self.neighbor = neighbor
             self.step = 0 # initialize simulation step
             self.initial_positions = initial_positions
-            #self.thermo_period = thermo_period # TOFIX to be added later
-            #self.dumping_period = dumping_period # TOFIX to be added later
 
 .. label:: end_InitializeSimulation_class
 
@@ -76,57 +68,11 @@ Several parameters are provided to the *InitializeSimulation* class:
 Finally, the *dimensions* of the system are calculated as the length of
 *box_dimensions*.
 
-Set a random seed
------------------
-
-Providing a *seed* allows for reproducible simulations. When a seed is
-specified, the simulation will produce the exact same results each time it
-is run, including identical atom positions and velocities. This can be
-particularly useful for debugging purposes.
-
-Add the following *if* condition to the *__init__()* method:
-
-.. label:: start_InitializeSimulation_class
-
-.. code-block:: python
-
-    def __init__(self,
-        (...)
-        self.initial_positions = initial_positions
-        if self.seed is not None:
-            np.random.seed(self.seed)
-
-.. label:: end_InitializeSimulation_class
-
-If a *seed* is provided, it is used to initialize the *random.seed()* function
-from *NumPy*. If *seed* is set to its default value of *None*, the simulation
-will proceed with randomization.
-
 Nondimensionalize units
 -----------------------
 
 Just like we did in :ref:`chapter2-label`, let us nondimensionalize the provided
-parameters, here the *box_dimensions* as well as the *initial_positions*:
-
-.. label:: start_InitializeSimulation_class
-
-.. code-block:: python
-
-    def nondimensionalize_units_1(self):
-        """Use LJ prefactors to convert units into non-dimensional."""
-        # Normalize box dimensions
-        box_dimensions = []
-        for L in self.box_dimensions:
-            box_dimensions.append(L/self.reference_distance)
-        self.box_dimensions = box_dimensions # errase the previously defined box_dimensions
-        # Normalize the box dimensions
-        if self.initial_positions is not None:
-            self.initial_positions = self.initial_positions/self.reference_distance
-        self.cut_off /= self.reference_distance
-
-.. label:: end_InitializeSimulation_class
-
-Let us call the *nondimensionalize_units_1* method from the *__init__* class:
+parameters, here the *box_dimensions*, the *cut_off*, as well as the *initial_positions*:
 
 .. label:: start_InitializeSimulation_class
 
@@ -134,9 +80,9 @@ Let us call the *nondimensionalize_units_1* method from the *__init__* class:
 
     def __init__(self,
         (...)
-        if self.seed is not None:
-            np.random.seed(self.seed)
-        self.nondimensionalize_units_1()
+        self.initial_positions = initial_positions
+        self.nondimensionalize_units(["box_dimensions", "cut_off",
+                                      "initial_positions"])
 
 .. label:: end_InitializeSimulation_class
 
@@ -151,9 +97,7 @@ method to the *InitializeSimulation* class:
 .. code-block:: python
 
     def define_box(self):
-        """Define the simulation box.
-        For 2D simulations, the third dimensions only contains 0.
-        """
+        """Define the simulation box. Only 3D boxes are supported."""
         box_boundaries = np.zeros((3, 2))
         dim = 0
         for L in self.box_dimensions:
@@ -182,7 +126,8 @@ Let us call the *define_box* method from the *__init__* class:
 
     def __init__(self,
         (...)
-        self.nondimensionalize_units_1()
+        self.nondimensionalize_units(["box_dimensions", "cut_off",
+                                      "initial_positions"])
         self.define_box()
 
 .. label:: end_InitializeSimulation_class
@@ -202,10 +147,10 @@ case, the array must be of size 'number of atoms' times 'number of dimensions'.
 
     def populate_box(self):
         if self.initial_positions is None:
-            atoms_positions = np.zeros((self.total_number_atoms, 3))
+            atoms_positions = np.zeros((np.sum(self.number_atoms), 3))
             for dim in np.arange(3):
                 diff_box = np.diff(self.box_boundaries[dim])
-                random_pos = np.random.random(self.total_number_atoms)
+                random_pos = np.random.random(np.sum(self.number_atoms))
                 atoms_positions[:, dim] = random_pos*diff_box-diff_box/2
             self.atoms_positions = atoms_positions
         else:
@@ -247,6 +192,7 @@ neighboring atoms, the simulation becomes more efficient. Add the following
 
     def update_neighbor_lists(self):
         if (self.step % self.neighbor == 0):
+            print(self.cut_off)
             matrix = distances.contact_matrix(self.atoms_positions,
                 cutoff=self.cut_off, #+2,
                 returntype="numpy",
@@ -281,7 +227,7 @@ more practical (see below).
             # Precalculte LJ cross-coefficients
             sigma_ij_list = []
             epsilon_ij_list = []
-            for Ni in np.arange(self.total_number_atoms-1): # tofix error for GCMC
+            for Ni in np.arange(np.sum(self.number_atoms)-1): # tofix error for GCMC
                 # Read information about atom i
                 sigma_i = self.atoms_sigma[Ni]
                 epsilon_i = self.atoms_epsilon[Ni]
@@ -344,14 +290,31 @@ boundaries along all 3 dimensions of space:
 
     import numpy as np
     from InitializeSimulation import InitializeSimulation
-
-    # Initialize the InitializeSimulation 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
+    # Define box size
+    L = 20*ureg.angstrom
+    # Define a cut off
+    rc = 2.5*sig_1
+
+    # Initialize the prepare object
     init = InitializeSimulation(
-        number_atoms=[2, 3],
-        epsilon=[0.2, 0.4], # kcal/mol
-        sigma=[3, 4], # A
-        atom_mass=[10, 20], # g/mol
-        box_dimensions=[20, 20, 20], # A
+        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
+        box_dimensions=[L, L, L], # A
+        cut_off=rc,
     )
 
     # Test function using pytest
diff --git a/tests/build-documentation.py b/tests/build-documentation.py
index fbaa9bc..79cc619 100644
--- a/tests/build-documentation.py
+++ b/tests/build-documentation.py
@@ -18,7 +18,7 @@
     os.mkdir("generated-codes/")
 
 # loop on the different chapter
-for chapter_id in [1, 2]:
+for chapter_id in [1, 2, 3]:
     # for each chapter, create the corresponding code
     RST_EXISTS, created_tests, folder = sphinx_to_python(path_to_docs, chapter_id)
     if RST_EXISTS:

From a612c417b02149e0d146475bb882aeddda19c967 Mon Sep 17 00:00:00 2001
From: Simon Gravelle <simon.gravelle@live.fr>
Date: Sat, 31 Aug 2024 22:11:43 +0200
Subject: [PATCH 03/13] continue simplifying unit conversion system

---
 docs/source/chapters/chapter3.rst |  3 +-
 docs/source/chapters/chapter4.rst | 88 ++++++++++++-------------------
 docs/source/chapters/chapter5.rst | 88 +++++++++++++++++++++++--------
 tests/build-documentation.py      |  2 +-
 4 files changed, 103 insertions(+), 78 deletions(-)

diff --git a/docs/source/chapters/chapter3.rst b/docs/source/chapters/chapter3.rst
index 1fe44b2..8f40ffd 100644
--- a/docs/source/chapters/chapter3.rst
+++ b/docs/source/chapters/chapter3.rst
@@ -23,7 +23,7 @@ Let us improve the previously created *InitializeSimulation* class:
 
 .. code-block:: python
 
-    class InitializeSimulation(Prepare):
+    class InitializeSimulation(Prepare, Utilities):
         def __init__(self,
                     box_dimensions,  # List - Angstroms
                     cut_off, # Angstroms
@@ -192,7 +192,6 @@ neighboring atoms, the simulation becomes more efficient. Add the following
 
     def update_neighbor_lists(self):
         if (self.step % self.neighbor == 0):
-            print(self.cut_off)
             matrix = distances.contact_matrix(self.atoms_positions,
                 cutoff=self.cut_off, #+2,
                 returntype="numpy",
diff --git a/docs/source/chapters/chapter4.rst b/docs/source/chapters/chapter4.rst
index 63c3b1a..e1cf34c 100644
--- a/docs/source/chapters/chapter4.rst
+++ b/docs/source/chapters/chapter4.rst
@@ -68,22 +68,17 @@ Then, let us fill the *__init__()* method:
     class MinimizeEnergy(Measurements):
         def __init__(self,
                     maximum_steps,
-                    cut_off=9,
-                    neighbor=1,
-                    displacement=0.01,
                     thermo_outputs="MaxF",
-                    data_folder=None,
+                    data_folder="Outputs/",
                     *args,
                     **kwargs):
-            self.displacement = displacement
             self.maximum_steps = maximum_steps
             self.thermo_outputs = thermo_outputs
             self.data_folder = data_folder
-            if self.data_folder is not None:
-                if os.path.exists(self.data_folder) is False:
-                    os.mkdir(self.data_folder)
+            if os.path.exists(self.data_folder) is False:
+                os.mkdir(self.data_folder)
             super().__init__(*args, **kwargs)
-
+            
 .. label:: end_MinimizeEnergy_class
 
 An important parameter is *maximum_steps*, which sets the maximum number
@@ -94,37 +89,6 @@ displacement value.
 The *thermo_outputs* and *data_folder* parameters are used for printing data
 to files. These two parameters will be useful in the next chapter, :ref:`chapter5-label`.
 
-Nondimensionalize units
------------------------
-
-As was done previously, some parameters from the *MinimizeEnergy* class
-must be non-dimensionalized: *cut_off* and *displacement*. Add the following
-method to the *MinimizeEnergy* class:
-
-.. label:: start_MinimizeEnergy_class
-
-.. code-block:: python
-
-    def nondimensionalize_units_2(self):
-        """Use LJ prefactors to convert units into non-dimensional."""
-        self.displacement = self.displacement/self.reference_distance
-
-.. label:: end_MinimizeEnergy_class
-
-Let us call the *nondimensionalize_units_2()* method from the *__init__()*
-method:
-
-.. label:: start_MinimizeEnergy_class
-
-.. code-block:: python
-
-    def __init__(self,
-        (...)
-        super().__init__(*args, **kwargs)
-        self.nondimensionalize_units_2()
-
-.. label:: end_MinimizeEnergy_class
-
 Energy minimizer
 ----------------
 
@@ -136,6 +100,7 @@ following *run()* method to the *MinimizeEnergy* class:
 .. code-block:: python
 
     def run(self):
+        self.displacement = 0.01 # pick a random initial displacement (dimentionless)
         # *step* loops for 0 to *maximum_steps*+1
         for self.step in range(0, self.maximum_steps+1):
             # First, meevaluate the initial energy and max force
@@ -185,7 +150,7 @@ class:
     def compute_potential(self):
         """Compute the potential energy by summing up all pair contributions."""
         energy_potential = 0
-        for Ni in np.arange(self.total_number_atoms-1):
+        for Ni in np.arange(np.sum(self.number_atoms)-1):
             # Read neighbor list
             neighbor_of_i = self.neighbor_lists[Ni]
             # Measure distance
@@ -234,11 +199,11 @@ let us create a new method that is dedicated solely to measuring forces:
 
     def compute_force(self, return_vector = True):
         if return_vector: # return a N-size vector
-            force_vector = np.zeros((self.total_number_atoms,3))
+            force_vector = np.zeros((np.sum(self.number_atoms),3))
         else: # return a N x N matrix
-            force_matrix = np.zeros((self.total_number_atoms,
-                                    self.total_number_atoms,3))
-        for Ni in np.arange(self.total_number_atoms-1):
+            force_matrix = np.zeros((np.sum(self.number_atoms),
+                                    np.sum(self.number_atoms),3))
+        for Ni in np.arange(np.sum(self.number_atoms)-1):
             # Read neighbor list
             neighbor_of_i = self.neighbor_lists[Ni]
             # Measure distance
@@ -282,7 +247,7 @@ within the *Utilities* class:
 .. code-block:: python
 
     def wrap_in_box(self):
-        for dim in np.arange(self.dimensions):
+        for dim in np.arange(3):
             out_ids = self.atoms_positions[:, dim] \
                 > self.box_boundaries[dim][1]
             self.atoms_positions[:, dim][out_ids] \
@@ -306,15 +271,32 @@ typically negative.
 .. code-block:: python
 
     from MinimizeEnergy import MinimizeEnergy
-
-    # Initialize the MinimizeEnergy object and run the minimization
+    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
+    # Define box size
+    L = 20*ureg.angstrom
+    # Define a cut off
+    rc = 2.5*sig_1
+
+    # Initialize the prepare object
     minimizer = MinimizeEnergy(
+        ureg = ureg,
         maximum_steps=100,
-        number_atoms=[2, 3],
-        epsilon=[0.2, 0.4], # kcal/mol
-        sigma=[3, 4], # A
-        atom_mass=[10, 20], # g/mol
-        box_dimensions=[20, 20, 20], # A
+        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
+        box_dimensions=[L, L, L], # A
+        cut_off=rc,
     )
     minimizer.run()
 
diff --git a/docs/source/chapters/chapter5.rst b/docs/source/chapters/chapter5.rst
index b15b715..b1c7f97 100644
--- a/docs/source/chapters/chapter5.rst
+++ b/docs/source/chapters/chapter5.rst
@@ -85,9 +85,9 @@ All quantities are re-dimensionalized before getting outputed.
             if code.step % code.thermo_period == 0:
                 if code.step == 0:
                     Epot = code.compute_potential() \
-                        * code.reference_energy  # kcal/mol
+                        * code.ref_energy  # kcal/mol
                 else:
-                    Epot = code.Epot * code.reference_energy  # kcal/mol
+                    Epot = code.Epot * code.ref_energy  # kcal/mol
                 if code.step == 0:
                     if code.thermo_outputs == "Epot":
                         logger.info(f"step Epot")
@@ -101,7 +101,7 @@ All quantities are re-dimensionalized before getting outputed.
                     logger.info(f"{code.step} {Epot:.2f} {code.MaxF:.2f}")
                 elif code.thermo_outputs == "Epot-press":
                     code.calculate_pressure()
-                    press = code.pressure * code.reference_pressure  # Atm
+                    press = code.pressure * code.ref_pressure  # Atm
                     logger.info(f"{code.step} {Epot:.2f} {press:.2f}")
 
 .. label:: end_logger_class
@@ -125,14 +125,12 @@ a LAMMPS dump format, and can be read by molecular dynamics softwares like VMD.
         if code.dumping_period is not None:
             if code.step % code.dumping_period == 0:
                 # Convert units to the *real* dimensions
-                box_boundaries = code.box_boundaries\
-                    * code.reference_distance # Angstrom
-                atoms_positions = code.atoms_positions\
-                    * code.reference_distance # Angstrom
+                box_boundaries = code.box_boundaries*code.ref_length # Angstrom
+                atoms_positions = code.atoms_positions*code.ref_length # Angstrom
                 atoms_types = code.atoms_type
                 if velocity:
                     atoms_velocities = code.atoms_velocities \
-                        * code.reference_distance/code.reference_time # Angstrom/femtosecond
+                        * code.ref_length/code.ref_time # Angstrom/femtosecond
                 # Start writting the file
                 if code.step == 0: # Create new file
                     f = open(code.data_folder + filename, "w")
@@ -141,18 +139,18 @@ a LAMMPS dump format, and can be read by molecular dynamics softwares like VMD.
                 f.write("ITEM: TIMESTEP\n")
                 f.write(str(code.step) + "\n")
                 f.write("ITEM: NUMBER OF ATOMS\n")
-                f.write(str(code.total_number_atoms) + "\n")
+                f.write(str(np.sum(code.number_atoms)) + "\n")
                 f.write("ITEM: BOX BOUNDS pp pp pp\n")
                 for dim in np.arange(3):
-                    f.write(str(box_boundaries[dim][0]) + " "
-                            + str(box_boundaries[dim][1]) + "\n")
+                    f.write(str(box_boundaries[dim][0].magnitude) + " "
+                            + str(box_boundaries[dim][1].magnitude) + "\n")
                 cpt = 1
                 if velocity:
                     f.write("ITEM: ATOMS id type x y z vx vy vz\n")
                     characters = "%d %d %.3f %.3f %.3f %.3f %.3f %.3f %s"
                     for type, xyz, vxyz in zip(atoms_types,
-                                            atoms_positions,
-                                            atoms_velocities):
+                                            atoms_positions.magnitude,
+                                            atoms_velocities.magnitude):
                         v = [cpt, type, xyz[0], xyz[1], xyz[2],
                                 vxyz[0], vxyz[1], vxyz[2]]
                         f.write(characters % (v[0], v[1], v[2], v[3], v[4],
@@ -162,7 +160,7 @@ a LAMMPS dump format, and can be read by molecular dynamics softwares like VMD.
                     f.write("ITEM: ATOMS id type x y z\n")
                     characters = "%d %d %.3f %.3f %.3f %s"
                     for type, xyz in zip(atoms_types,
-                                        atoms_positions):
+                                        atoms_positions.magnitude):
                         v = [cpt, type, xyz[0], xyz[1], xyz[2]]
                         f.write(characters % (v[0], v[1], v[2],
                                             v[3], v[4], '\n'))
@@ -197,6 +195,35 @@ Add the same lines at the top of the *MinimizeEnergy.py* file:
 
 .. label:: end_MinimizeEnergy_class
 
+Finally, let us make sure that *thermo_period*, *dumping_period*, and *thermo_outputs*
+parameters are passed the InitializeSimulation method:
+
+.. label:: start_InitializeSimulation_class
+
+.. code-block:: python
+
+    def __init__(self,
+            (...)
+                neighbor=1, # Integer
+                thermo_period = None,
+                dumping_period = None,
+                thermo_outputs = None,
+
+.. label:: end_InitializeSimulation_class
+
+.. label:: start_InitializeSimulation_class
+
+.. code-block:: python
+
+    def __init__(self,
+        (...)
+        self.initial_positions = initial_positions
+        self.thermo_period = thermo_period
+        self.dumping_period = dumping_period
+        self.thermo_outputs = thermo_outputs
+
+.. label:: end_InitializeSimulation_class
+
 Test the code
 -------------
 
@@ -208,21 +235,38 @@ files were indeed created without the *Outputs/* folder:
 
 .. code-block:: python
 
-    import os
     from MinimizeEnergy import MinimizeEnergy
+    from pint import UnitRegistry
+    ureg = UnitRegistry()
+    import os
 
-    # Initialize the MinimizeEnergy object and run the minimization
+    # 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
+    # Define box size
+    L = 20*ureg.angstrom
+    # Define a cut off
+    rc = 2.5*sig_1
+
+    # Initialize the prepare object
     minimizer = MinimizeEnergy(
+        ureg = ureg,
         maximum_steps=100,
         thermo_period=25,
         dumping_period=25,
-        thermo_outputs="Epot-MaxF",
-        number_atoms=[2, 3],
-        epsilon=[0.1, 1.0], # kcal/mol
-        sigma=[3, 4], # A
-        atom_mass=[10, 20], # g/mol
-        box_dimensions=[20, 20, 20], # A
+        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
+        box_dimensions=[L, L, L], # A
+        cut_off=rc,
         data_folder="Outputs/",
+        thermo_outputs="Epot-MaxF",
     )
     minimizer.run()
 
diff --git a/tests/build-documentation.py b/tests/build-documentation.py
index 79cc619..adb8619 100644
--- a/tests/build-documentation.py
+++ b/tests/build-documentation.py
@@ -18,7 +18,7 @@
     os.mkdir("generated-codes/")
 
 # loop on the different chapter
-for chapter_id in [1, 2, 3]:
+for chapter_id in [1, 2, 3, 4, 5]:
     # for each chapter, create the corresponding code
     RST_EXISTS, created_tests, folder = sphinx_to_python(path_to_docs, chapter_id)
     if RST_EXISTS:

From 199a099b07b3d89bb832b03e38cd137aeb7796a3 Mon Sep 17 00:00:00 2001
From: Simon Gravelle <simon.gravelle@live.fr>
Date: Sat, 31 Aug 2024 22:35:54 +0200
Subject: [PATCH 04/13] continue simplify parameter import

---
 docs/source/chapters/chapter1.rst |  7 ++-
 docs/source/chapters/chapter2.rst |  2 +-
 docs/source/chapters/chapter4.rst |  4 --
 docs/source/chapters/chapter5.rst |  4 ++
 docs/source/chapters/chapter6.rst | 75 ++++++++++++++-----------------
 tests/build-documentation.py      |  2 +-
 6 files changed, 45 insertions(+), 49 deletions(-)

diff --git a/docs/source/chapters/chapter1.rst b/docs/source/chapters/chapter1.rst
index 3fd1b8a..822dd92 100644
--- a/docs/source/chapters/chapter1.rst
+++ b/docs/source/chapters/chapter1.rst
@@ -167,6 +167,7 @@ Within the *InitializeSimulation.py* file, copy the following lines:
 
 .. code-block:: python
 
+    import os
     import numpy as np
     from Prepare import Prepare
     from Utilities import Utilities
@@ -298,12 +299,14 @@ and copy the following lines into it:
 
     # Make sure that MonteCarlo correctly inherits from Utilities
     def test_montecarlo_inherits_from_utilities():
-        assert issubclass(MonteCarlo, Utilities), "MonteCarlo should inherit from Utilities"
+        assert issubclass(MonteCarlo, Utilities), \
+            "MonteCarlo should inherit from Utilities"
         print("MonteCarlo correctly inherits from Utilities")
 
     # Make sure that Utilities does not inherit from MonteCarlo
     def test_utilities_does_not_inherit_from_montecarlo():
-        assert not issubclass(Utilities, MonteCarlo), "Utilities should not inherit from MonteCarlo"
+        assert not issubclass(Utilities, MonteCarlo), \
+            "Utilities should not inherit from MonteCarlo"
         print("Utilities does not inherit from MonteCarlo, as expected")
 
     # In the script is launched with Python, call Pytest
diff --git a/docs/source/chapters/chapter2.rst b/docs/source/chapters/chapter2.rst
index 468edb4..f6682d4 100644
--- a/docs/source/chapters/chapter2.rst
+++ b/docs/source/chapters/chapter2.rst
@@ -163,7 +163,7 @@ unit system to the *LJ* unit system:
             /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_mass, self.ref_time, self.ref_pressure, self.ref_temperature]
         self.ref_units = [ref.units for ref in self.ref_quantities]
 
 .. label:: end_Prepare_class
diff --git a/docs/source/chapters/chapter4.rst b/docs/source/chapters/chapter4.rst
index e1cf34c..6178904 100644
--- a/docs/source/chapters/chapter4.rst
+++ b/docs/source/chapters/chapter4.rst
@@ -69,14 +69,10 @@ Then, let us fill the *__init__()* method:
         def __init__(self,
                     maximum_steps,
                     thermo_outputs="MaxF",
-                    data_folder="Outputs/",
                     *args,
                     **kwargs):
             self.maximum_steps = maximum_steps
             self.thermo_outputs = thermo_outputs
-            self.data_folder = data_folder
-            if os.path.exists(self.data_folder) is False:
-                os.mkdir(self.data_folder)
             super().__init__(*args, **kwargs)
             
 .. label:: end_MinimizeEnergy_class
diff --git a/docs/source/chapters/chapter5.rst b/docs/source/chapters/chapter5.rst
index b1c7f97..848e90d 100644
--- a/docs/source/chapters/chapter5.rst
+++ b/docs/source/chapters/chapter5.rst
@@ -208,6 +208,7 @@ parameters are passed the InitializeSimulation method:
                 thermo_period = None,
                 dumping_period = None,
                 thermo_outputs = None,
+                data_folder="Outputs/",
 
 .. label:: end_InitializeSimulation_class
 
@@ -221,6 +222,9 @@ parameters are passed the InitializeSimulation method:
         self.thermo_period = thermo_period
         self.dumping_period = dumping_period
         self.thermo_outputs = thermo_outputs
+        self.data_folder = data_folder
+        if os.path.exists(self.data_folder) is False:
+            os.mkdir(self.data_folder)
 
 .. label:: end_InitializeSimulation_class
 
diff --git a/docs/source/chapters/chapter6.rst b/docs/source/chapters/chapter6.rst
index 7daaaf1..149732a 100644
--- a/docs/source/chapters/chapter6.rst
+++ b/docs/source/chapters/chapter6.rst
@@ -96,44 +96,17 @@ and the desired temperature (:math:`T`). Let us add these parameters to the
     class MonteCarlo(Measurements):
         def __init__(self,
                     maximum_steps,
-                    cut_off = 9,
+                    desired_temperature,
                     displace_mc = None,
-                    neighbor = 1,
-                    desired_temperature = 300,
                     thermo_outputs = "press",
-                    data_folder = None,
                     *args,
                     **kwargs):
             self.maximum_steps = maximum_steps
-            self.cut_off = cut_off
             self.displace_mc = displace_mc
-            self.neighbor = neighbor
             self.desired_temperature = desired_temperature
             self.thermo_outputs = thermo_outputs
-            self.data_folder = data_folder
-            if self.data_folder is not None:
-                if os.path.exists(self.data_folder) is False:
-                    os.mkdir(self.data_folder)
             super().__init__(*args, **kwargs)
-            self.nondimensionalize_units_3()
-
-.. label:: end_MonteCarlo_class
-
-Here, we anticipate that some of the parameters have to be nondimensionalized, which
-is done with the *nondimensionalize_units_3* method that must also be added to
-the *MonteCarlo* class:
-
-.. label:: start_MonteCarlo_class
-
-.. code-block:: python
-
-    def nondimensionalize_units_3(self):
-        """Use LJ prefactors to convert units into non-dimensional."""
-        self.cut_off = self.cut_off/self.reference_distance
-        self.desired_temperature = self.desired_temperature \
-            /self.reference_temperature
-        if self.displace_mc is not None:
-            self.displace_mc /= self.reference_distance
+            self.nondimensionalize_units(["desired_temperature", "displace_mc"])
 
 .. label:: end_MonteCarlo_class
 
@@ -201,21 +174,41 @@ One can use a similar test as previously. Let us use a displace distance of
 
 .. code-block:: python
 
-    import os
     from MonteCarlo import MonteCarlo
+    from pint import UnitRegistry
+    ureg = UnitRegistry()
+    import os
 
-    mc = MonteCarlo(maximum_steps=1000,
-        dumping_period=100,
+    # Define atom number of each group
+    nmb_1= 50
+    # Define LJ parameters (sigma)
+    sig_1 = 3*ureg.angstrom
+    # Define LJ parameters (epsilon)
+    eps_1 = 0.1*ureg.kcal/ureg.mol
+    # Define atom mass
+    mss_1 = 10*ureg.gram/ureg.mol
+    # Define box size
+    L = 20*ureg.angstrom
+    # Define a cut off
+    rc = 2.5*sig_1
+    # Pick the desired temperature
+    T = 300*ureg.kelvin
+
+    # Initialize the prepare object
+    mc = MonteCarlo(
+        ureg = ureg,
+        maximum_steps=1000,
         thermo_period=100,
-        thermo_outputs = "Epot",
-        displace_mc = 0.5,
-        number_atoms=[50],
-        epsilon=[0.1], # kcal/mol
-        sigma=[3], # A
-        atom_mass=[10], # g/mol
-        box_dimensions=[20, 20, 20], # A
-        data_folder="Outputs/",
-        )
+        dumping_period=100,
+        number_atoms=[nmb_1],
+        epsilon=[eps_1], # kcal/mol
+        sigma=[sig_1], # A
+        atom_mass=[mss_1], # g/mol
+        box_dimensions=[L, L, L], # A
+        cut_off=rc,
+        thermo_outputs="Epot",
+        desired_temperature=T, # K
+    )
     mc.run()
 
 .. label:: end_test_6a_class
diff --git a/tests/build-documentation.py b/tests/build-documentation.py
index adb8619..ea64d66 100644
--- a/tests/build-documentation.py
+++ b/tests/build-documentation.py
@@ -18,7 +18,7 @@
     os.mkdir("generated-codes/")
 
 # loop on the different chapter
-for chapter_id in [1, 2, 3, 4, 5]:
+for chapter_id in [1, 2, 3, 4, 5, 6]:
     # for each chapter, create the corresponding code
     RST_EXISTS, created_tests, folder = sphinx_to_python(path_to_docs, chapter_id)
     if RST_EXISTS:

From 042742f646c35fd6ea3f214bd6259b142d2a02c7 Mon Sep 17 00:00:00 2001
From: Simon Gravelle <simon.gravelle@live.fr>
Date: Sat, 31 Aug 2024 23:07:12 +0200
Subject: [PATCH 05/13] fixed Monte Carlo units import

---
 docs/source/chapters/chapter4.rst |  2 --
 docs/source/chapters/chapter5.rst |  6 +++---
 docs/source/chapters/chapter6.rst | 22 +++++++++-------------
 3 files changed, 12 insertions(+), 18 deletions(-)

diff --git a/docs/source/chapters/chapter4.rst b/docs/source/chapters/chapter4.rst
index 6178904..0f1a439 100644
--- a/docs/source/chapters/chapter4.rst
+++ b/docs/source/chapters/chapter4.rst
@@ -68,11 +68,9 @@ Then, let us fill the *__init__()* method:
     class MinimizeEnergy(Measurements):
         def __init__(self,
                     maximum_steps,
-                    thermo_outputs="MaxF",
                     *args,
                     **kwargs):
             self.maximum_steps = maximum_steps
-            self.thermo_outputs = thermo_outputs
             super().__init__(*args, **kwargs)
             
 .. label:: end_MinimizeEnergy_class
diff --git a/docs/source/chapters/chapter5.rst b/docs/source/chapters/chapter5.rst
index 848e90d..12f46b7 100644
--- a/docs/source/chapters/chapter5.rst
+++ b/docs/source/chapters/chapter5.rst
@@ -96,13 +96,13 @@ All quantities are re-dimensionalized before getting outputed.
                     elif code.thermo_outputs == "Epot-press":
                         logger.info(f"step Epot press")
                 if code.thermo_outputs == "Epot":
-                    logger.info(f"{code.step} {Epot:.2f}")
+                    logger.info(f"{code.step} {Epot.magnitude:.2f}")
                 elif code.thermo_outputs == "Epot-MaxF":
-                    logger.info(f"{code.step} {Epot:.2f} {code.MaxF:.2f}")
+                    logger.info(f"{code.step} {Epot.magnitude:.2f} {code.MaxF:.2f}")
                 elif code.thermo_outputs == "Epot-press":
                     code.calculate_pressure()
                     press = code.pressure * code.ref_pressure  # Atm
-                    logger.info(f"{code.step} {Epot:.2f} {press:.2f}")
+                    logger.info(f"{code.step} {Epot.magnitude:.2f} {press.magnitude:.2f}")
 
 .. label:: end_logger_class
 
diff --git a/docs/source/chapters/chapter6.rst b/docs/source/chapters/chapter6.rst
index 149732a..475ab44 100644
--- a/docs/source/chapters/chapter6.rst
+++ b/docs/source/chapters/chapter6.rst
@@ -49,22 +49,17 @@ Let us add a method named *monte_carlo_move* to the *MonteCarlo* class:
         if self.displace_mc is not None: # only trigger if displace_mc was provided by the user
             # If needed, recalculate neighbor/coeff lists
             self.update_neighbor_lists()
-            self.identify_atom_properties()
             self.update_cross_coefficients()
-            try: # try using the last saved Epot, if it exists
-                initial_Epot = self.Epot
-            except: # If self.Epot does not exists yet, calculate it
-                initial_Epot = self.compute_potential()
+            if hasattr(self, 'Epot') is False: # If self.Epot does not exists yet, calculate it
+                self.Epot = self.compute_potential()
+            initial_Epot = self.Epot
             # Make a copy of the initial atoms positions
             initial_positions = copy.deepcopy(self.atoms_positions)
             # Pick an atom id randomly
-            atom_id = np.random.randint(self.total_number_atoms)
+            atom_id = np.random.randint(np.sum(self.number_atoms))
             # Move the chosen atom in a random direction
             # The maximum displacement is set by self.displace_mc
-            if self.dimensions == 3:
-                move = (np.random.random(3)-0.5)*self.displace_mc 
-            elif self.dimensions == 2: # the third value will be 0
-                move = np.append((np.random.random(2) - 0.5) * self.displace_mc, 0)
+            move = (np.random.random(3)-0.5)*self.displace_mc 
             self.atoms_positions[atom_id] += move
             # Measure the optential energy of the new configuration
             trial_Epot = self.compute_potential()
@@ -76,7 +71,6 @@ Let us add a method named *monte_carlo_move* to the *MonteCarlo* class:
             if random_number <= acceptation_probability: # Accept new position
                 self.Epot = trial_Epot
             else: # Reject new position
-                self.Epot = initial_Epot
                 self.atoms_positions = initial_positions # Revert to initial positions
 
 .. label:: end_MonteCarlo_class
@@ -98,13 +92,11 @@ and the desired temperature (:math:`T`). Let us add these parameters to the
                     maximum_steps,
                     desired_temperature,
                     displace_mc = None,
-                    thermo_outputs = "press",
                     *args,
                     **kwargs):
             self.maximum_steps = maximum_steps
             self.displace_mc = displace_mc
             self.desired_temperature = desired_temperature
-            self.thermo_outputs = thermo_outputs
             super().__init__(*args, **kwargs)
             self.nondimensionalize_units(["desired_temperature", "displace_mc"])
 
@@ -193,6 +185,8 @@ One can use a similar test as previously. Let us use a displace distance of
     rc = 2.5*sig_1
     # Pick the desired temperature
     T = 300*ureg.kelvin
+    # choose the displace_mc
+    displace_mc = sig_1/4
 
     # Initialize the prepare object
     mc = MonteCarlo(
@@ -208,6 +202,8 @@ One can use a similar test as previously. Let us use a displace distance of
         cut_off=rc,
         thermo_outputs="Epot",
         desired_temperature=T, # K
+        neighbor=20,
+        displace_mc = displace_mc,
     )
     mc.run()
 

From d0a6352aa067bef64ef425c715c8b5fc7f185908 Mon Sep 17 00:00:00 2001
From: Simon Gravelle <simon.gravelle@live.fr>
Date: Sat, 31 Aug 2024 23:22:56 +0200
Subject: [PATCH 06/13] added successful_move and failled move counter

---
 docs/source/chapters/chapter6.rst | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/docs/source/chapters/chapter6.rst b/docs/source/chapters/chapter6.rst
index 475ab44..ff810ec 100644
--- a/docs/source/chapters/chapter6.rst
+++ b/docs/source/chapters/chapter6.rst
@@ -70,8 +70,10 @@ Let us add a method named *monte_carlo_move* to the *MonteCarlo* class:
             acceptation_probability = np.min([1, np.exp(-beta*delta_E)])
             if random_number <= acceptation_probability: # Accept new position
                 self.Epot = trial_Epot
+                self.successful_move += 1
             else: # Reject new position
                 self.atoms_positions = initial_positions # Revert to initial positions
+                self.failed_move += 1
 
 .. label:: end_MonteCarlo_class
 
@@ -99,6 +101,8 @@ and the desired temperature (:math:`T`). Let us add these parameters to the
             self.desired_temperature = desired_temperature
             super().__init__(*args, **kwargs)
             self.nondimensionalize_units(["desired_temperature", "displace_mc"])
+            self.successful_move = 0
+            self.failed_move = 0
 
 .. label:: end_MonteCarlo_class
 

From cc4b43fffdf012d6a15046f3b0ade546627ac1c1 Mon Sep 17 00:00:00 2001
From: Simon Gravelle <simon.gravelle@live.fr>
Date: Sun, 1 Sep 2024 10:20:25 +0200
Subject: [PATCH 07/13] continue fixing units

---
 docs/source/chapters/chapter1.rst |   1 +
 docs/source/chapters/chapter6.rst |  35 ++++----
 docs/source/chapters/chapter7.rst | 145 ++++++++++++------------------
 docs/source/chapters/chapter8.rst | 129 ++++++++++++++++++--------
 tests/build-documentation.py      |   2 +-
 5 files changed, 170 insertions(+), 142 deletions(-)

diff --git a/docs/source/chapters/chapter1.rst b/docs/source/chapters/chapter1.rst
index 822dd92..5540209 100644
--- a/docs/source/chapters/chapter1.rst
+++ b/docs/source/chapters/chapter1.rst
@@ -192,6 +192,7 @@ Within the *Measurements.py* file, copy the following lines:
 
 .. code-block:: python
 
+    import numpy as np
     from InitializeSimulation import InitializeSimulation
 
 
diff --git a/docs/source/chapters/chapter6.rst b/docs/source/chapters/chapter6.rst
index ff810ec..6e4beec 100644
--- a/docs/source/chapters/chapter6.rst
+++ b/docs/source/chapters/chapter6.rst
@@ -145,21 +145,6 @@ Let us add a dump too:
 
 .. label:: end_MonteCarlo_class
 
-To output the density, let us add the following method to the *Utilities* class:
-# TOFIX: not used yet
-
-.. label:: start_Utilities_class
-
-.. code-block:: python
-
-    def calculate_density(self):
-        """Calculate the mass density."""
-        volume = np.prod(self.box_size[:3])  # Unitless
-        total_mass = np.sum(self.atoms_mass)  # Unitless
-        return total_mass/volume  # Unitless
-
-.. label:: end_Utilities_class
-
 Test the code
 -------------
 
@@ -195,9 +180,9 @@ One can use a similar test as previously. Let us use a displace distance of
     # Initialize the prepare object
     mc = MonteCarlo(
         ureg = ureg,
-        maximum_steps=1000,
-        thermo_period=100,
-        dumping_period=100,
+        maximum_steps=100,
+        thermo_period=10,
+        dumping_period=10,
         number_atoms=[nmb_1],
         epsilon=[eps_1], # kcal/mol
         sigma=[sig_1], # A
@@ -209,8 +194,22 @@ One can use a similar test as previously. Let us use a displace distance of
         neighbor=20,
         displace_mc = displace_mc,
     )
+
+    # Run the Monte Carlo simulation
     mc.run()
 
+    # Test function using pytest
+    def test_output_files():
+        assert os.path.exists("Outputs/dump.mc.lammpstrj"), "Test failed: dump file was not created"
+        assert os.path.exists("Outputs/simulation.log"), "Test failed: log file was not created"
+        print("Test passed")
+
+    # If the script is run directly, execute the tests
+    if __name__ == "__main__":
+        import pytest
+        # Run pytest programmatically
+        pytest.main(["-s", __file__])
+
 .. label:: end_test_6a_class
 
 The evolution of the potential energy as a function of the number of steps
diff --git a/docs/source/chapters/chapter7.rst b/docs/source/chapters/chapter7.rst
index eb4d26c..d634e7a 100644
--- a/docs/source/chapters/chapter7.rst
+++ b/docs/source/chapters/chapter7.rst
@@ -30,9 +30,9 @@ vector direction between atoms pairs will have to be written here.
 Implement the Virial equation
 -----------------------------
 
-Let us add the following method to the *Utilities* class.
+Let us add the following method to the *Measurements* class.
 
-.. label:: start_Utilities_class
+.. label:: start_Measurements_class
 
 .. code-block:: python
 
@@ -40,21 +40,24 @@ Let us add the following method to the *Utilities* class.
         """Evaluate p based on the Virial equation (Eq. 4.4.2 in Frenkel-Smit,
         Understanding molecular simulation: from algorithms to applications, 2002)"""
         # Compute the ideal contribution
-        Ndof = self.dimensions*self.total_number_atoms-self.dimensions    
-        volume = np.prod(self.box_size[:self.dimensions])
+        Nat = np.sum(self.number_atoms) # total number of atoms
+        dimension = 3 # 3D is the only possibility here
+        Ndof = dimension*Nat-dimension    
+        volume = np.prod(self.box_size[:3]) # box volume
         try:
             self.calculate_temperature() # this is for later on, when velocities are computed
             temperature = self.temperature
         except:
             temperature = self.desired_temperature # for MC, simply use the desired temperature
-        p_ideal = Ndof*temperature/(volume*self.dimensions)
+        p_ideal = Ndof*temperature/(volume*dimension)
         # Compute the non-ideal contribution
-        distances_forces = np.sum(self.compute_force(return_vector = False)*self.evaluate_rij_matrix())
-        p_nonideal = distances_forces/(volume*self.dimensions)
+        distances_forces = np.sum(self.compute_force(return_vector = False) \
+                                  *self.evaluate_rij_matrix())
+        p_nonideal = distances_forces/(volume*dimension)
         # Final pressure
         self.pressure = p_ideal+p_nonideal
 
-.. label:: end_Utilities_class
+.. label:: end_Measurements_class
 
 To evaluate all the vectors between all the particles, let us also add the
 *evaluate_rij_matrix* method to the *Utilities* class:
@@ -65,95 +68,75 @@ To evaluate all the vectors between all the particles, let us also add the
 
     def evaluate_rij_matrix(self):
         """Matrix of vectors between all particles."""
-        box_size = self.box_size[:3]
-        half_box_size = self.box_size[:3]/2.0
-        rij_matrix = np.zeros((self.total_number_atoms,self.total_number_atoms,3))
-        positions_j = self.atoms_positions
-        for Ni in range(self.total_number_atoms-1):
-            position_i = self.atoms_positions[Ni]
-            rij_xyz = (np.remainder(position_i - positions_j + half_box_size, box_size) - half_box_size)
+        Nat = np.sum(self.number_atoms)
+        Box = self.box_size[:3]
+        rij_matrix = np.zeros((Nat, Nat,3))
+        pos_j = self.atoms_positions
+        for Ni in range(Nat-1):
+            pos_i = self.atoms_positions[Ni]
+            rij_xyz = (np.remainder(pos_i - pos_j + Box/2.0, Box) - Box/2.0)
             rij_matrix[Ni] = rij_xyz
-        # use nan_to_num to avoid "nan" values in 2D
-        return np.nan_to_num(rij_matrix)
+        return rij_matrix
 
 .. label:: end_Utilities_class
 
 Test the code
 -------------
 
-Let us test the outputed pressure. An interesting test is to contront the output from
-our code with some data from the litterature. Let us used the same parameters
-as in Ref. :cite:`woodMonteCarloEquation1957`, where Monte Carlo simulations
-are used to simulate argon bulk phase. All we have to do is to apply our current
-code using their parameters, i.e. :math:`\sigma = 3.405~\text{Å}`, :math:`\epsilon = 0.238~\text{kcal/mol}`,
-and :math:`T = 55~^\circ\text{C}`. More details are given in the first illustration, :ref:`project1-label`.
-
-On the side note, a relatively small cut-off as well as a small number of atoms were
-chosen to make the calculation faster. 
+Let us test the outputed pressure. 
 
 .. label:: start_test_7a_class
 
 .. code-block:: python
 
-    import numpy as np
     from MonteCarlo import MonteCarlo
-    from scipy import constants as cst
     from pint import UnitRegistry
-    from reader import import_data
-
-    # Initialize the unit registry
     ureg = UnitRegistry()
-
-    # Constants
-    kB = cst.Boltzmann * ureg.J / ureg.kelvin  # Boltzmann constant
-    Na = cst.Avogadro / ureg.mole  # Avogadro's number
-    R = kB * Na  # Gas constant
-
-    # Parameters taken from Wood1957
-    tau = 2  # Ratio between volume / reduced volume
-    epsilon = (119.76 * ureg.kelvin * kB * Na).to(ureg.kcal / ureg.mol)  # kcal/mol
-    r_star = 3.822 * ureg.angstrom  # Angstrom
-    sigma = r_star / 2**(1/6)  # Angstrom
-    N_atom = 50  # Number of atoms
-    m_argon = 39.948 * ureg.gram / ureg.mol  # g/mol
-    T = 328.15 * ureg.degK  # 328 K or 55°C
-    volume_star = r_star**3 * Na * 2**(-0.5)
-    cut_off = sigma * 2.5  # Angstrom
-    displace_mc = sigma / 5  # Angstrom
-    volume = N_atom * volume_star * tau / Na  # Angstrom^3
-    box_size = volume**(1/3)  # Angstrom
-
-    # Initialize and run the Monte Carlo simulation
+    import os
+
+    # Define atom number of each group
+    nmb_1= 50
+    # Define LJ parameters (sigma)
+    sig_1 = 3*ureg.angstrom
+    # Define LJ parameters (epsilon)
+    eps_1 = 0.1*ureg.kcal/ureg.mol
+    # Define atom mass
+    mss_1 = 10*ureg.gram/ureg.mol
+    # Define box size
+    L = 20*ureg.angstrom
+    # Define a cut off
+    rc = 2.5*sig_1
+    # Pick the desired temperature
+    T = 300*ureg.kelvin
+    # choose the displace_mc
+    displace_mc = sig_1/4
+
+    # Initialize the prepare object
     mc = MonteCarlo(
-        maximum_steps=15000,
-        dumping_period=1000,
-        thermo_period=1000,
+        ureg = ureg,
+        maximum_steps=100,
+        thermo_period=10,
+        dumping_period=10,
+        number_atoms=[nmb_1],
+        epsilon=[eps_1], # kcal/mol
+        sigma=[sig_1], # A
+        atom_mass=[mss_1], # g/mol
+        box_dimensions=[L, L, L], # A
+        cut_off=rc,
         thermo_outputs="Epot-press",
-        neighbor=50,
-        number_atoms=[N_atom],
-        epsilon=[epsilon.magnitude],
-        sigma=[sigma.magnitude],
-        atom_mass=[m_argon.magnitude],
-        box_dimensions=[box_size.magnitude, box_size.magnitude, box_size.magnitude],
-        displace_mc=displace_mc.magnitude,
-        desired_temperature=T.magnitude,
-        cut_off=cut_off.magnitude,
-        data_folder="Outputs/",
+        desired_temperature=T, # K
+        neighbor=20,
+        displace_mc = displace_mc,
     )
+
+    # Run the Monte Carlo simulation
     mc.run()
 
     # Test function using pytest
-    def test_pV_over_RT():
-        # Import the data and calculate pV / RT
-        pressure_data = np.array(import_data("Outputs/simulation.log")[1:])[0][:,2][10:]  # Skip initial values
-        pressure_atm = np.mean(pressure_data) # atm
-        pressure_Pa = (pressure_atm * ureg.atm).to(ureg.pascal) # Pa
-        volume = (volume_star * tau / Na).to(ureg.meter**3) # m3
-        pV_over_RT = (pressure_Pa * volume / (R * T) * Na).magnitude
-        
-        # Assert that pV_over_RT is close to 1.5
-        assert np.isclose(pV_over_RT, 1.5, atol=1.0), f"Test failed: pV/Rt = {pV_over_RT}, expected close to 1.5"
-        print(f"pV/RT = {pV_over_RT:.2f} --- (The expected value from Wood1957 is 1.5)")
+    def test_output_files():
+        assert os.path.exists("Outputs/dump.mc.lammpstrj"), "Test failed: dump file was not created"
+        assert os.path.exists("Outputs/simulation.log"), "Test failed: log file was not created"
+        print("Test passed")
 
     # If the script is run directly, execute the tests
     if __name__ == "__main__":
@@ -162,13 +145,3 @@ chosen to make the calculation faster.
         pytest.main(["-s", __file__])
 
 .. label:: end_test_7a_class
-
-Which should return a value for :math:`p V / R T` that is close to the expected value
-of 1.5 by Wood and Parker for :math:`\tau = V/V^* = 2` (see Fig. 4 in Ref. :cite:`woodMonteCarloEquation1957`):
-
-.. code-block:: bw
-
-    (...)
-    p v / R T = 1.56  --- (The expected value from Wood1957 is 1.5)
-
-The exact value will vary from one simulation to the other due to noise.
\ No newline at end of file
diff --git a/docs/source/chapters/chapter8.rst b/docs/source/chapters/chapter8.rst
index 9982fc4..e1bf6cc 100644
--- a/docs/source/chapters/chapter8.rst
+++ b/docs/source/chapters/chapter8.rst
@@ -45,18 +45,29 @@ Let us add the following method called *monte_carlo_exchange* to the *MonteCarlo
         initial_sigma_lists = copy.deepcopy(self.sigma_ij_list)
         initial_epsilon_lists = copy.deepcopy(self.epsilon_ij_list)
         # Apply a 50-50 probability to insert or delete
-        if np.random.random() < 0.5:
+        insert_or_delete = np.random.random()
+        if np.random.random() < insert_or_delete:
             self.monte_carlo_insert()
         else:
             self.monte_carlo_delete()
-        # If np.random.random() < self.acceptation_probability, do nothing
-        if np.random.random() > self.acceptation_probability:
+        if np.random.random() < self.acceptation_probability: # accepted move
+            # Update the success counters
+            if np.random.random() < insert_or_delete:
+                self.successful_insert += 1
+            else:
+                self.successful_delete += 1
+        else:
             # Reject the new position, revert to inital position
             self.neighbor_lists = initial_neighbor_lists
             self.sigma_ij_list = initial_sigma_lists
             self.epsilon_ij_list = initial_epsilon_lists
             self.atoms_positions = initial_positions
             self.number_atoms = initial_number_atoms
+            # Update the failed counters
+            if np.random.random() < insert_or_delete:
+                self.failed_insert += 1
+            else:
+                self.failed_delete += 1
 
 .. label:: end_MonteCarlo_class
 
@@ -143,8 +154,6 @@ Complete the *__init__* method as follows:
                     displace_mc = None,
                     desired_mu = None,
                     inserted_type = 0,
-                    desired_temperature = 300,
-                    (...)
 
 .. label:: end_MonteCarlo_class
 
@@ -160,23 +169,26 @@ and
             self.displace_mc = displace_mc
             self.desired_mu = desired_mu
             self.inserted_type = inserted_type
-            self.desired_temperature = desired_temperature
-            (...)
 
 .. label:: end_MonteCarlo_class
 
-Let us non-dimentionalize desired_mu by adding:
+Let us also normalised the "desired_mu":
 
 .. label:: start_MonteCarlo_class
 
 .. code-block:: python
 
-    def nondimensionalize_units_3(self):
-        (...)
-            self.displace_mc /= self.reference_distance
-        if self.desired_mu is not None:
-            self.desired_mu /= self.reference_energy
-        (...)
+    class MonteCarlo(Outputs):
+        def __init__(self,
+            (...)
+            self.nondimensionalize_units(["desired_temperature", "displace_mc"])
+            self.nondimensionalize_units(["desired_mu"])
+            self.successful_move = 0
+            self.failed_move = 0
+            self.successful_insert = 0
+            self.failed_insert = 0
+            self.successful_delete = 0
+            self.failed_delete = 0
 
 .. label:: end_MonteCarlo_class
 
@@ -203,47 +215,90 @@ We need to calculate Lambda:
 
     def calculate_Lambda(self, mass):
         """Estimate the de Broglie wavelength."""
-        # Is it normal that mass is unitless ???
         m = mass/cst.Avogadro*cst.milli  # kg
-        kB = cst.Boltzmann  # J/K
-        Na = cst.Avogadro
-        kB *= Na/cst.calorie/cst.kilo #  kCal/mol/K
         T = self.desired_temperature  # N
-        T_K = T*self.reference_energy/kB  # K
-        Lambda = cst.h/np.sqrt(2*np.pi*kB*m*T_K)  # m
-        Lambda /= cst.angstrom  # Angstrom
-        return Lambda / self.reference_distance # dimensionless
+        return 1/np.sqrt(2*np.pi*mass*T)
 
 .. label:: end_MonteCarlo_class
 
+To output the density, let us add the following method to the *Measurements* class:
+
+.. label:: start_Measurements_class
+
+.. code-block:: python
+
+    def calculate_density(self):
+        """Calculate the mass density."""
+        volume = np.prod(self.box_size[:3])  # Unitless
+        total_mass = np.sum(self.atoms_mass)  # Unitless
+        return total_mass/volume  # Unitless
+
+.. label:: end_Measurements_class
+
 Test the code
 -------------
 
 One can use a similar test as previously. Let us use a displace distance of
 0.5 Angstrom, and make 1000 steps.
 
-.. label:: start_test_MonteCarlo_class
+.. label:: start_test_8a_class
 
 .. code-block:: python
 
-    import os
     from MonteCarlo import MonteCarlo
+    from pint import UnitRegistry
+    ureg = UnitRegistry()
+    import os
 
-    mc = MonteCarlo(maximum_steps=1000,
-        dumping_period=100,
-        thermo_period=100,
-        desired_mu = -3,
-        inserted_type = 0,
-        number_atoms=[50],
-        epsilon=[0.1], # kcal/mol
-        sigma=[3], # A
-        atom_mass=[1], # g/mol
-        box_dimensions=[20, 20, 20], # A
-        data_folder = "outputs/",
-        )
+    # Define atom number of each group
+    nmb_1= 50
+    # Define LJ parameters (sigma)
+    sig_1 = 3*ureg.angstrom
+    # Define LJ parameters (epsilon)
+    eps_1 = 0.1*ureg.kcal/ureg.mol
+    # Define atom mass
+    mss_1 = 10*ureg.gram/ureg.mol
+    # Define box size
+    L = 20*ureg.angstrom
+    # Define a cut off
+    rc = 2.5*sig_1
+    # Pick the desired temperature
+    T = 300*ureg.kelvin
+    # choose the desired_mu
+    desired_mu = -3*ureg.kcal/ureg.mol
+
+    # Initialize the prepare object
+    mc = MonteCarlo(
+        ureg = ureg,
+        maximum_steps=100,
+        thermo_period=10,
+        dumping_period=10,
+        number_atoms=[nmb_1],
+        epsilon=[eps_1], # kcal/mol
+        sigma=[sig_1], # A
+        atom_mass=[mss_1], # g/mol
+        box_dimensions=[L, L, L], # A
+        cut_off=rc,
+        thermo_outputs="Epot-press",
+        desired_temperature=T, # K
+        neighbor=1,
+        desired_mu = desired_mu,
+    )
     mc.run()
 
-.. label:: end_test_MonteCarlo_class
+    # Test function using pytest
+    def test_output_files():
+        assert os.path.exists("Outputs/dump.mc.lammpstrj"), "Test failed: dump file was not created"
+        assert os.path.exists("Outputs/simulation.log"), "Test failed: log file was not created"
+        print("Test passed")
+
+    # If the script is run directly, execute the tests
+    if __name__ == "__main__":
+        import pytest
+        # Run pytest programmatically
+        pytest.main(["-s", __file__])
+
+.. label:: end_test_8a_class
 
 The evolution of the potential energy as a function of the
 number of steps are written in the *Outputs/Epot.dat* file
diff --git a/tests/build-documentation.py b/tests/build-documentation.py
index ea64d66..b6c5395 100644
--- a/tests/build-documentation.py
+++ b/tests/build-documentation.py
@@ -18,7 +18,7 @@
     os.mkdir("generated-codes/")
 
 # loop on the different chapter
-for chapter_id in [1, 2, 3, 4, 5, 6]:
+for chapter_id in [1, 2, 3, 4, 5, 6, 7, 8]:
     # for each chapter, create the corresponding code
     RST_EXISTS, created_tests, folder = sphinx_to_python(path_to_docs, chapter_id)
     if RST_EXISTS:

From 694644bd950c52d1b604624a4641026e0e377357 Mon Sep 17 00:00:00 2001
From: Simon Gravelle <simon.gravelle@live.fr>
Date: Sun, 1 Sep 2024 11:47:00 +0200
Subject: [PATCH 08/13] imporved text and grammar

---
 docs/source/chapters/chapter1.rst | 101 +++++++++++++++---------------
 docs/source/chapters/chapter2.rst |  28 ++++-----
 docs/source/chapters/chapter4.rst |  20 +-----
 docs/source/chapters/chapter8.rst |   1 -
 docs/source/journal-article.bib   |  25 +++++++-
 5 files changed, 91 insertions(+), 84 deletions(-)

diff --git a/docs/source/chapters/chapter1.rst b/docs/source/chapters/chapter1.rst
index 5540209..632d616 100644
--- a/docs/source/chapters/chapter1.rst
+++ b/docs/source/chapters/chapter1.rst
@@ -39,32 +39,32 @@ containing either Python functions or classes:
 
    * - File Name 
      - Content
-   * - *Prepare.py* 
+   * - Prepare.py 
      - *Prepare* class: Methods for preparing the non-dimensionalization of the
        units
-   * - *Utilities.py* 
+   * - Utilities.py
      - *Utilities* class: General-purpose methods, inherited by all other
        classes
-   * - *InitializeSimulation.py*
+   * - InitializeSimulation.py
      - *InitializeSimulation* class: Methods necessary to set up the system and
        prepare the simulation, inherited by all the classes below
-   * - *MinimizeEnergy.py* 
+   * - MinimizeEnergy.py
      - *MinimizeEnergy* class: Methods for performing energy minimization
-   * - *MonteCarlo.py*
+   * - MonteCarlo.py
      - *MonteCarlo* class: Methods for performing Monte Carlo simulations in
        different ensembles (e.g., Grand Canonical, Canonical)
-   * - *MolecularDynamics.py*
+   * - MolecularDynamics.py
      - *MolecularDynamics* class: Methods for performing molecular dynamics in
        different ensembles (NVE, NPT, NVT)
-   * - *measurements.py* 
-     - Functions for performing specific measurements on the system
+   * - Measurements.py
+     - *Measurements* class: Methods for for performing specific measurements on the system
    * - *potentials.py* 
      - Functions for calculating the potentials and forces between atoms
-   * - *logger.py*
+   * - logger.py
      - Functions for outputting data into text files
-   * - *dumper.py*
+   * - dumper.py
      - Functions for outputting data into trajectory files for visualization
-   * - *reader.py*
+   * - reader.py
      - Functions for importing data from text files
 
 Some of these files are created in this chapter; others will be created later
@@ -73,14 +73,13 @@ on. All of these files must be created within the same folder.
 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
+In molecular simulations, potential functions are used to model the interaction
+between atoms. Although more complex options exist, potentials are usually
 defined as functions of the distance :math:`r` between two atoms.
 
-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.
+Create a file named *potentials.py*. This file will contain a function called
+*potentials*. For now, the only potential that can be returned by this function
+is the Lennard-Jones (LJ) potential, but this may change in the future.
 
 Copy the following lines into *potentials.py*:
 
@@ -88,16 +87,11 @@ Copy the following lines into *potentials.py*:
 
 .. code-block:: python
 
-    import numpy as np
-
-    def potentials(potential_type, epsilon, sigma, r, derivative=False):
-        if potential_type == "Lennard-Jones":
-            if derivative:
-                return 48 * epsilon * ((sigma / r) ** 12 - 0.5 * (sigma / r) ** 6) / r
-            else:
-                return 4 * epsilon * ((sigma / r) ** 12 - (sigma / r) ** 6)
+    def potentials(epsilon, sigma, r, derivative=False):
+        if derivative:
+            return 48 * epsilon * ((sigma / r) ** 12 - 0.5 * (sigma / r) ** 6) / r
         else:
-            raise ValueError(f"Unknown potential type: {potential_type}")
+            return 4 * epsilon * ((sigma / r) ** 12 - (sigma / r) ** 6)
 
 .. label:: end_potentials_class
 
@@ -108,14 +102,19 @@ 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],
+    - \frac{1}{2} \left( \frac{\sigma}{r} \right)^6 \right], ~ \text{for} ~ r < r_\text{c},
 
 or the potential energy:
 
 .. math::
 
     U_\text{LJ} = 4 \epsilon \left[ \left( \frac{\sigma}{r} \right)^{12}
-    - \left( \frac{\sigma}{r} \right)^6 \right].
+    - \left( \frac{\sigma}{r} \right)^6 \right], ~ \text{for} ~ r < r_\text{c}.
+
+Here, :math:`\sigma` is the distance at which the potential :math:`U_\text{LJ}`
+is zero, :math:`\epsilon` is the depth of the potential well, and
+:math:`r_\text{c}` is a cutoff distance. For :math:`r > r_\text{c}`,
+:math:`U_\text{LJ} = 0` and :math:`F_\text{LJ} = 0`.
 
 Create the Classes
 ------------------
@@ -183,8 +182,14 @@ Within the *InitializeSimulation.py* file, copy the following lines:
 .. label:: end_InitializeSimulation_class
 
 The *InitializeSimulation* class inherits from the previously created
-*Prepare* and Utilities classes. Additionally, we anticipate that *NumPy* will
-be required.
+*Prepare* and *Utilities* classes. Additionally, we anticipate that |NumPy|
+will be required :cite:`harris2020array`. We also anticipate that the *os*
+module, which provides a way to interact with the operating system, will
+be required :cite:`Rossum2009Python3`.
+
+.. |NumPy| raw:: html
+
+   <a href="https://numpy.org/" target="_blank">NumPy</a>
 
 Within the *Measurements.py* file, copy the following lines:
 
@@ -204,13 +209,13 @@ Within the *Measurements.py* file, copy the following lines:
           
 .. label:: end_Measurements_class
 
-The *Measurements* class inherits both the *InitializeSimulation* and
-*Utilities* classes. 
+The *Measurements* class inherits from *InitializeSimulation* (and thus
+also inherits from the *Prepare* and *Utilities* classes).
 
-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*.
+Finally, let us create the three remaining classes: *MinimizeEnergy*,
+*MonteCarlo*, and *MolecularDynamics*. Each of these classes inherits
+from the *Measurements* class (and thus also from the *Prepare*, *Utilities*,
+and *InitializeSimulation* classes).
 
 Within the *MinimizeEnergy.py* file, copy the following lines:
 
@@ -219,6 +224,8 @@ Within the *MinimizeEnergy.py* file, copy the following lines:
 .. code-block:: python
 
     from Measurements import Measurements
+    import numpy as np
+    import copy
     import os
 
 
@@ -230,8 +237,8 @@ Within the *MinimizeEnergy.py* file, copy the following lines:
 
 .. label:: end_MinimizeEnergy_class
 
-We anticipate that the *os* module, which provides a way to interact with the
-operating system, will be required :cite:`Rossum2009Python3`.
+The *copy* library, which provides functions to create shallow or deep copies of
+objects, is imported, along with *NumPy* and *os*.
 
 Within the *MonteCarlo.py* file, copy the following lines:
 
@@ -239,10 +246,8 @@ Within the *MonteCarlo.py* file, copy the following lines:
 
 .. code-block:: python
 
-    from scipy import constants as cst
     import numpy as np
     import copy
-    import os
     from Measurements import Measurements
 
     import warnings
@@ -257,12 +262,10 @@ Within the *MonteCarlo.py* file, copy the following lines:
 
 .. label:: end_MonteCarlo_class
 
-Several libraries were imported, namely *Constants* from *SciPy*, *NumPy*, *copy*
-and *os*.
-
-The *warnings* was placed to avoid the anoying message "*RuntimeWarning: overflow
-encountered in exp*" that is sometimes triggered by the exponential of the
-*acceptation_probability* (see :ref:`chapter6-label`).
+The *ignore warnings* commands are optional; they were added to avoid the
+annoying message "*RuntimeWarning: overflow encountered in exp*" that is sometimes
+triggered by the exponential function of *acceptation_probability* (see the
+:ref:`chapter6-label` chapter).
 
 Finally, within the *MolecularDynamics.py* file, copy the following lines:
 
@@ -331,9 +334,9 @@ Alternatively, this test can also be launched using *Pytest* by typing in a term
 
     pytest .
 
-We can also test that calling the *__init__*
-method of the *MonteCarlo* class does not return any error. In new Python file
-called *test_1b.py*, copy the following lines:
+We can also test that calling the *__init__* method of the *MonteCarlo* class
+does not return any error. In new Python file called *test_1b.py*, copy the
+following lines:
 
 .. label:: start_test_1b_class
 
diff --git a/docs/source/chapters/chapter2.rst b/docs/source/chapters/chapter2.rst
index f6682d4..2714e19 100644
--- a/docs/source/chapters/chapter2.rst
+++ b/docs/source/chapters/chapter2.rst
@@ -29,7 +29,7 @@ The *real* unit system follows the conventions outlined in the |lammps-unit-syst
 - Forces are in (kcal/mol)/Ångström,
 - Temperature is in Kelvin,
 - Pressure is in atmospheres,
-- Density is in g/cm\ :sup:`3` (in 3D).
+- Density is in g/cm\ :sup:`3`.
 
 .. |lammps-unit-systems| raw:: html
 
@@ -56,12 +56,8 @@ where :math:`k_\text{B}` is the Boltzmann constant.
 Start coding
 ------------
 
-Let's fill in the previously created class named Prepare. To facilitate unit
-conversion, we will import |NumPy| and the constants module from |SciPy|.
-
-.. |NumPy| raw:: html
-
-   <a href="https://numpy.org/" target="_blank">NumPy</a>
+Let's fill in the previously created class named *Prepare*. To facilitate unit
+conversion, we will import NumPy and the constants module from |SciPy|.
 
 .. |SciPy| raw:: html
 
@@ -78,7 +74,7 @@ In the file named *Prepare.py*, add the following lines:
 
 .. label:: end_Prepare_class
 
-Four parameters are provided to the *Prepare* class:
+Four atom parameters are provided to the *Prepare* class:
 
 - the atom masses :math:`m`,
 - the LJ parameters :math:`\sigma` and :math:`\epsilon`,
@@ -100,7 +96,6 @@ Modify the *Prepare* class as follows:
                     epsilon, # List - Kcal/mol
                     sigma, # List - Angstrom
                     atom_mass,  # List - g/mol
-                    potential_type="Lennard-Jones",
                     *args,
                     **kwargs):
             self.ureg = ureg
@@ -108,18 +103,19 @@ Modify the *Prepare* class as follows:
             self.epsilon = epsilon
             self.sigma = sigma
             self.atom_mass = atom_mass
-            self.potential_type = potential_type
             super().__init__(*args, **kwargs)
 
 .. label:: end_Prepare_class
 
-Here, the four lists *number_atoms* :math:`N`, *epsilon* :math:`\epsilon`,
-*sigma* :math:`\sigma`, and *atom_mass* :math:`m` are given default values of
-:math:`10`, :math:`0.1~\text{[Kcal/mol]}`, :math:`3~\text{[Å]}`, and
-:math:`10~\text{[g/mol]}`, respectively.
+Here, *number_atoms* :math:`N`, *epsilon* :math:`\epsilon`,
+*sigma* :math:`\sigma`, and *atom_mass* :math:`m` must be provided as lists
+where the elements have no units, kcal/mol, angstrom, and g/mol units,
+respectively. The units will be enforced with the |Pint| unit registry, *ureg*,
+which must also be provided as a parameter.
+
+.. |Pint| raw:: html
 
-The type of potential is also specified, with Lennard-Jones being chosen as
-the default option.
+   <a href="https://pint.readthedocs.io" target="_blank">Pint</a>
 
 All the parameters are assigned to *self*, allowing other methods to access
 them. The *args* and *kwargs* parameters are used to accept an arbitrary number
diff --git a/docs/source/chapters/chapter4.rst b/docs/source/chapters/chapter4.rst
index 0f1a439..0746a69 100644
--- a/docs/source/chapters/chapter4.rst
+++ b/docs/source/chapters/chapter4.rst
@@ -47,19 +47,7 @@ minimized energy state.
 Prepare the minimization
 ------------------------
 
-Let us start by importing NumPy and the copy libraries. Add the following
-to the beginning of the *MinimizeEnergy.py* file:
-
-.. label:: start_MinimizeEnergy_class
-
-.. code-block:: python
-
-    import numpy as np
-    import copy
-
-.. label:: end_MinimizeEnergy_class
-
-Then, let us fill the *__init__()* method:
+Let us fill the *__init__()* method:
 
 .. label:: start_MinimizeEnergy_class
 
@@ -154,8 +142,7 @@ class:
             # Measure potential using information about cross coefficients
             sigma_ij = self.sigma_ij_list[Ni]
             epsilon_ij = self.epsilon_ij_list[Ni]
-            energy_potential += np.sum(potentials(self.potential_type,
-                                                  epsilon_ij, sigma_ij, rij))
+            energy_potential += np.sum(potentials(epsilon_ij, sigma_ij, rij))
         return energy_potential
     
 .. label:: end_Utilities_class
@@ -207,8 +194,7 @@ let us create a new method that is dedicated solely to measuring forces:
             # Measure force using information about cross coefficients
             sigma_ij = self.sigma_ij_list[Ni]
             epsilon_ij = self.epsilon_ij_list[Ni]       
-            fij_xyz = potentials(self.potential_type, epsilon_ij,
-                                 sigma_ij, rij, derivative = True)
+            fij_xyz = potentials(epsilon_ij, sigma_ij, rij, derivative = True)
             if return_vector:
                 # Add the contribution to both Ni and its neighbors
                 force_vector[Ni] += np.sum((fij_xyz*rij_xyz.T/rij).T, axis=0)
diff --git a/docs/source/chapters/chapter8.rst b/docs/source/chapters/chapter8.rst
index e1bf6cc..e4c9b8c 100644
--- a/docs/source/chapters/chapter8.rst
+++ b/docs/source/chapters/chapter8.rst
@@ -215,7 +215,6 @@ We need to calculate Lambda:
 
     def calculate_Lambda(self, mass):
         """Estimate the de Broglie wavelength."""
-        m = mass/cst.Avogadro*cst.milli  # kg
         T = self.desired_temperature  # N
         return 1/np.sqrt(2*np.pi*mass*T)
 
diff --git a/docs/source/journal-article.bib b/docs/source/journal-article.bib
index 7148bbc..8869996 100644
--- a/docs/source/journal-article.bib
+++ b/docs/source/journal-article.bib
@@ -61,4 +61,27 @@ @book{Rossum2009Python3
  isbn = {1441412697},
  publisher = {CreateSpace},
  address = {Scotts Valley, CA}
-} 
\ No newline at end of file
+}
+
+@article{         harris2020array,
+ title         = {Array programming with {NumPy}},
+ author        = {Charles R. Harris and K. Jarrod Millman and St{\'{e}}fan J.
+                 van der Walt and Ralf Gommers and Pauli Virtanen and David
+                 Cournapeau and Eric Wieser and Julian Taylor and Sebastian
+                 Berg and Nathaniel J. Smith and Robert Kern and Matti Picus
+                 and Stephan Hoyer and Marten H. van Kerkwijk and Matthew
+                 Brett and Allan Haldane and Jaime Fern{\'{a}}ndez del
+                 R{\'{i}}o and Mark Wiebe and Pearu Peterson and Pierre
+                 G{\'{e}}rard-Marchant and Kevin Sheppard and Tyler Reddy and
+                 Warren Weckesser and Hameer Abbasi and Christoph Gohlke and
+                 Travis E. Oliphant},
+ year          = {2020},
+ month         = sep,
+ journal       = {Nature},
+ volume        = {585},
+ number        = {7825},
+ pages         = {357--362},
+ doi           = {10.1038/s41586-020-2649-2},
+ publisher     = {Springer Science and Business Media {LLC}},
+ url           = {https://doi.org/10.1038/s41586-020-2649-2}
+}
\ No newline at end of file

From 277678f4088bb65f796c64a040b6b6e597dc4a97 Mon Sep 17 00:00:00 2001
From: Simon Gravelle <simon.gravelle@live.fr>
Date: Sun, 1 Sep 2024 12:08:24 +0200
Subject: [PATCH 09/13] more typo correction

---
 docs/source/chapters/chapter2.rst | 66 +++++++++++++++++--------------
 1 file changed, 37 insertions(+), 29 deletions(-)

diff --git a/docs/source/chapters/chapter2.rst b/docs/source/chapters/chapter2.rst
index 2714e19..3b1ccdd 100644
--- a/docs/source/chapters/chapter2.rst
+++ b/docs/source/chapters/chapter2.rst
@@ -111,7 +111,7 @@ Here, *number_atoms* :math:`N`, *epsilon* :math:`\epsilon`,
 *sigma* :math:`\sigma`, and *atom_mass* :math:`m` must be provided as lists
 where the elements have no units, kcal/mol, angstrom, and g/mol units,
 respectively. The units will be enforced with the |Pint| unit registry, *ureg*,
-which must also be provided as a parameter.
+which must also be provided as a parameter (more details later).
 
 .. |Pint| raw:: html
 
@@ -124,9 +124,9 @@ of positional and keyword arguments, respectively.
 Calculate LJ units prefactors
 -----------------------------
 
-Within the *Prepare* class, let us create a method called *calculate_LJunits_prefactors*
-that will be used to calculate the prefactors necessary to convert units from the *real*
-unit system to the *LJ* unit system:
+Within the *Prepare* class, create a method called *calculate_LJunits_prefactors*
+that will be used to calculate the prefactors necessary to convert units from the
+*real* unit system to the *LJ* unit system:
 
 .. label:: start_Prepare_class
 
@@ -134,7 +134,7 @@ unit system to the *LJ* unit system:
 
     def calculate_LJunits_prefactors(self):
         """Calculate the Lennard-Jones units prefactors."""
-        # First define constants
+        # First define Boltzmann and Avogadro constants
         kB = cst.Boltzmann*cst.Avogadro/cst.calorie/cst.kilo  # kcal/mol/K
         kB *= self.ureg.kcal/self.ureg.mol/self.ureg.kelvin
         Na = cst.Avogadro/self.ureg.mol
@@ -157,22 +157,24 @@ unit system to the *LJ* unit system:
         # 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
+        # Group all the reference quantities into a list for practicality
         self.ref_quantities = [self.ref_length, self.ref_energy,
             self.ref_mass, self.ref_time, self.ref_pressure, self.ref_temperature]
         self.ref_units = [ref.units for ref in self.ref_quantities]
 
 .. label:: end_Prepare_class
 
-This method defines the reference distance as the first element in the
-*sigma* list, i.e., :math:`\sigma_{11}`. Therefore, atoms of type one will
-always be used for the normalization. Similarly, the first element
-in the *epsilon* list (:math:`\epsilon_{11}`) is used as the reference energy,
-and the first element in the *atom_mass* list (:math:`m_1`) is used as the
-reference mass. Then, the reference_time in femtoseconds is calculated
-as :math:`\sqrt{m_1 \sigma_{11}^2 / \epsilon_{11}}`, the reference temperature
-in Kelvin as :math:`\epsilon_{11} / k_\text{B}`, and the reference_pressure
-in atmospheres is calculated as :math:`\epsilon_{11}/\sigma_{11}^3`.
+This method defines the reference length as the first element in the
+*sigma* list, i.e., :math:`\sigma_{11}`. Therefore, atoms of type 1 will
+always be used for normalization. Similarly, the first element in the
+*epsilon* list (:math:`\epsilon_{11}`) is used as the reference energy, and
+the first element in the *atom_mass* list (:math:`m_1`) is used as the
+reference mass.
+
+The reference time in femtoseconds is then calculated as
+:math:`\sqrt{m_1 \sigma_{11}^2 / \epsilon_{11}}`, the reference
+temperature in Kelvin as :math:`\epsilon_{11} / k_\text{B}`, and the
+reference pressure in atmospheres as :math:`\epsilon_{11} / \sigma_{11}^3`.
 
 Finally, let us ensure that the *calculate_LJunits_prefactors* method is
 called systematically by adding the following line to the *__init__()* method:
@@ -233,11 +235,11 @@ class:
 
 .. label:: end_Prepare_class
 
-The index *0* is used to differentiate this method from other methods that
-will be used to nondimensionalize units in future classes. We anticipate that
-*epsilon*, *sigma*, and *atom_mass* may contain more than one element, so
-each element is normalized with the corresponding reference value. The
-*zip()* function allows us to loop over all three lists at once.
+When a *quantities_to_normalise* list containing parameter names is provided
+to the *nondimensionalize_units* method, a loop is performed over all the
+quantities. The value of each quantity is extracted using *getattr*. The units
+of the quantity of interest are then detected and normalized by the appropriate
+reference quantities defined by *calculate_LJunits_prefactors*.
 
 Let us also call the *nondimensionalize_units* from the *__init__()* method
 of the *Prepare* class:
@@ -253,22 +255,25 @@ of the *Prepare* class:
 
 .. label:: end_Prepare_class
 
-Identify atom properties
+Here, the *epsilon*, *sigma*, and *atom_mass* parameters will be
+nondimensionalized.
+
+Identify Atom Properties
 ------------------------
 
 Anticipating the future use of multiple atom types, where each type will be
 associated with its own :math:`\sigma`, :math:`\epsilon`, and :math:`m`, let
 us create arrays containing the properties of each atom in the simulation. For
-instance, in the case of a simulation with two atoms of type 1 and three atoms
-of type 2, the corresponding *atoms_sigma* array will be:
+instance, in a simulation with two atoms of type 1 and three atoms of type 2,
+the corresponding *atoms_sigma* array will be:
 
 .. math::
 
     \text{atoms_sigma} = [\sigma_{11}, \sigma_{11}, \sigma_{22}, \sigma_{22}, \sigma_{22}]
 
 where :math:`\sigma_{11}` and :math:`\sigma_{22}` are the sigma values for
-atoms of type 1 and 2, respectively. The *atoms_sigma* array will allow for
-future calculations of force.
+atoms of type 1 and 2, respectively. The *atoms_sigma* array will facilitate
+future calculations of forces.
 
 Create a new method called *identify_atom_properties*, and place it
 within the *Prepare* class:
@@ -318,7 +323,8 @@ Test the code
 
 Let's test the *Prepare* class to make sure that it does what is expected.
 Here, a system containing 2 atoms of type 1, and 3 atoms of type 2 is
-prepared. LJs parameters and masses for each groups are also defined.
+prepared. LJs parameters and masses for each groups are also defined, and given
+physical units thanks to the *UnitRegistry* of Pint.
 
 .. label:: start_test_2a_class
 
@@ -351,7 +357,8 @@ prepared. LJs parameters and masses for each groups are also defined.
     def test_atoms_epsilon():
         expected = np.array([1., 1., 2., 2., 2.])
         result = prep.atoms_epsilon
-        assert np.array_equal(result, expected), f"Test failed: {result} != {expected}"
+        assert np.array_equal(result, expected), \
+        f"Test failed: {result} != {expected}"
         print("Test passed")
 
     # In the script is launched with Python, call Pytest
@@ -361,5 +368,6 @@ prepared. LJs parameters and masses for each groups are also defined.
 
 .. label:: end_test_2a_class
 
-This test assert that the generated *atoms_epsilon* array is consistent with
-its expected value (see the previous paragraphs).
+This test asserts that the generated *atoms_epsilon* array is consistent
+with its expected value (see the previous paragraphs). Note that the expected
+values are in LJ units, i.e., they have been nondimensionalized by :math:`\epsilon_1`.

From dbae49fc45adb999b6537b4ae6d60fb377a0e7da Mon Sep 17 00:00:00 2001
From: Simon Gravelle <simon.gravelle@live.fr>
Date: Sun, 1 Sep 2024 13:11:35 +0200
Subject: [PATCH 10/13] improved descripition chapter 3

---
 docs/source/chapters/chapter3.rst | 117 +++++++++++++++---------------
 docs/source/chapters/chapter4.rst |   7 +-
 docs/source/journal-article.bib   |  13 +++-
 3 files changed, 73 insertions(+), 64 deletions(-)

diff --git a/docs/source/chapters/chapter3.rst b/docs/source/chapters/chapter3.rst
index 8f40ffd..72a6bf6 100644
--- a/docs/source/chapters/chapter3.rst
+++ b/docs/source/chapters/chapter3.rst
@@ -15,9 +15,9 @@ Initialize the simulation
     :align: right
     :class: only-light
 
-Here, the *InitializeSimulation* class is created. This class is used to
-prepare the simulation box and populate the box randomly with atoms.
-Let us improve the previously created *InitializeSimulation* class:
+Here, the *InitializeSimulation* class is completed. This class is used to
+prepare the simulation box and populate it randomly with atoms. Improve the
+previously created *InitializeSimulation* class as follows:
 
 .. label:: start_InitializeSimulation_class
 
@@ -43,30 +43,20 @@ Let us improve the previously created *InitializeSimulation* class:
 
 Several parameters are provided to the *InitializeSimulation* class:
 
-- The first parameter is the box dimensions, which is a list with a length
-  corresponding to the dimension of the system. Each element of the list
-  corresponds to a dimension of the box in Ångström in the x, y, and z
-  directions, respectively.
-- Optionally, a seed can be provided as an integer. If the seed is given
-  by the user, the random number generator will always produce the same
-  displacements.
+- The box dimensions, which is a list with 3 elements. Each element of
+  the list corresponds to a dimension of the box in Ångström in the :math:`x`,
+  :math:`y`, and :math:`z`` directions, respectively.
+- The cutoff :math:`r_\text{c}` for the potential interaction.
 - Optionally, initial positions for the atoms can be provided as an array
   of length corresponding to the number of atoms. If *initial_positions* 
-  is left equal to *None*, positions will be randomly assigned to the
-  atoms (see below).
-- A cut off value with a default of 10 Ångströms is defined.
-- The *neighbor* parameter determines the interval between recalculations of
-  the neighbor lists. By default, a value of 1 is used, meaning that neighbor
-  list will be reconstructed every steps. In certain cases, the number will be
-  increased to recude the computation time.
-
-..
-    - Optionally, a thermo period and a dumping_period can be provided to  # TOFIX to be added later
-    control the outputs from the simulation (it will be implemented  # TOFIX to be added later
-    in :ref:`chapter5-label`).  # TOFIX to be added later
-
-Finally, the *dimensions* of the system are calculated as the length of
-*box_dimensions*.
+  is set to *None*, positions will be randomly assigned to the atoms (see
+  below).
+- The *neighbor* parameter determines the interval between recalculations
+  of the neighbor lists. By default, a value of 1 is used, meaning that the
+  neighbor list will be reconstructed every step. In certain cases, this
+  number may be increased to reduce computation time.
+
+Note that the simulation step is initialized to 0.
 
 Nondimensionalize units
 -----------------------
@@ -89,8 +79,8 @@ parameters, here the *box_dimensions*, the *cut_off*, as well as the *initial_po
 Define the box
 --------------
 
-Let us define the simulation box using the values from *box_dimensions*. Add the following
-method to the *InitializeSimulation* class:
+Let us define the simulation box using the values from *box_dimensions*. Add the
+following method to the *InitializeSimulation* class:
 
 .. label:: start_InitializeSimulation_class
 
@@ -99,10 +89,8 @@ method to the *InitializeSimulation* class:
     def define_box(self):
         """Define the simulation box. Only 3D boxes are supported."""
         box_boundaries = np.zeros((3, 2))
-        dim = 0
-        for L in self.box_dimensions:
+        for dim, L in enumerate(self.box_dimensions):
             box_boundaries[dim] = -L/2, L/2
-            dim += 1
         self.box_boundaries = box_boundaries
         box_size = np.diff(self.box_boundaries).reshape(3)
         box_geometry = np.array([90, 90, 90])
@@ -113,12 +101,16 @@ method to the *InitializeSimulation* class:
 The *box_boundaries* are calculated from the *box_dimensions*. They
 represent the lowest and highest coordinates in all directions. By symmetry,
 the box is centered at 0 along all axes. A *box_size* is also defined,
-following the MDAnalysis conventions: Lx, Ly, Lz, 90, 90, 90, where the
-last three numbers are angles in degrees. Values different from *90* for
-the angles would define a triclinic (non-orthogonal) box, which is not
-currently supported by the existing code.
+following the MDAnalysis |mdanalysis_box|: Lx, Ly, Lz, 90, 90, 90, where
+the last three numbers are angles in degrees :cite:`michaud2011mdanalysis`.
+Values different from *90* for the angles would define a triclinic (non-orthogonal)
+box, which is not currently supported by the existing code.
+
+.. |mdanalysis_box| raw:: html
+
+   <a href="https://docs.mdanalysis.org/stable/documentation_pages/transformations/boxdimensions.html" target="_blank">conventions</a>
 
-Let us call the *define_box* method from the *__init__* class:
+Let us call *define_box()* from the *__init__()* method:
 
 .. label:: start_InitializeSimulation_class
 
@@ -132,13 +124,13 @@ Let us call the *define_box* method from the *__init__* class:
 
 .. label:: end_InitializeSimulation_class
 
-Populate the box
+Populate the Box
 ----------------
 
 Here, the atoms are placed within the simulation box. If initial
 positions were not provided (i.e., *initial_positions = None*), atoms
-are placed randomly within the box. If *initial_positions* was provided
-as an array, the provided positions are used instead. Note that, in this
+are placed randomly within the box. If *initial_positions* is provided
+as an array, the provided positions are used instead. Note that in this
 case, the array must be of size 'number of atoms' times 'number of dimensions'.
 
 .. label:: start_InitializeSimulation_class
@@ -146,11 +138,12 @@ case, the array must be of size 'number of atoms' times 'number of dimensions'.
 .. code-block:: python
 
     def populate_box(self):
+        Nat = np.sum(self.number_atoms) # total number of atoms
         if self.initial_positions is None:
-            atoms_positions = np.zeros((np.sum(self.number_atoms), 3))
+            atoms_positions = np.zeros((Nat, 3))
             for dim in np.arange(3):
                 diff_box = np.diff(self.box_boundaries[dim])
-                random_pos = np.random.random(np.sum(self.number_atoms))
+                random_pos = np.random.random(Nat)
                 atoms_positions[:, dim] = random_pos*diff_box-diff_box/2
             self.atoms_positions = atoms_positions
         else:
@@ -158,14 +151,14 @@ case, the array must be of size 'number of atoms' times 'number of dimensions'.
 
 .. label:: end_InitializeSimulation_class
 
-In case *initial_positions is None*, and array is first created. Then, random
-positions that are constrained within the box boundaries are defined using the
-random function of NumPy. Note that, here, the newly added atoms are added
-randomly within the box, without taking care of avoiding overlaps with
-existing atoms. Overlaps will be dealt with using energy minimization,
-see :ref:`chapter4-label`.
+In case *initial_positions* is None, an array is first created. Then,
+random positions constrained within the box boundaries are defined using
+the random function from NumPy. Note that newly added atoms are placed
+randomly within the box, without considering overlaps with existing
+atoms. Overlaps will be addressed using energy minimization (see
+the :ref:`chapter4-label` chapter).
 
-Let us call the *populate_box* method from the *__init__* class:
+Let us call *populate_box* from the *__init__* method:
 
 .. label:: start_InitializeSimulation_class
 
@@ -178,7 +171,7 @@ Let us call the *populate_box* method from the *__init__* class:
 
 .. label:: end_InitializeSimulation_class
 
-Build neighbor lists
+Build Neighbor Lists
 --------------------
 
 In molecular simulations, it is common practice to identify neighboring atoms
@@ -210,12 +203,16 @@ atom, ensuring that only relevant interactions are considered in the
 calculations. These lists will be recalculated at intervals specified by
 the *neighbor* input parameter.
 
-Update cross coefficients
+For efficiency, the *contact_matrix* function from MDAnalysis is
+used :cite:`michaud2011mdanalysis`. The *contact_matrix* function returns
+information about atoms located at a distance less than the cutoff from one another.
+
+Update Cross Coefficients
 -------------------------
 
-At the same time as the neighbor lists are getting build up, let us also
-pre-calculate the cross coefficients. This will make the force calculation
-more practical (see below).
+As the neighbor lists are being built, let us also pre-calculate the cross
+coefficients. This will make the force calculation more efficient (see
+below).
 
 .. label:: start_Utilities_class
 
@@ -250,7 +247,7 @@ Here, the values of the cross coefficients between atom of type 1 and 2,
     \sigma_{12} = (\sigma_{11}+\sigma_{22})/2 \\
     \epsilon_{12} = (\epsilon_{11}+\epsilon_{22})/2
 
-Finally, import the following library in the *Utilities.py* file:
+Finally, import the following libraries in the *Utilities.py* file:
 
 .. label:: start_Utilities_class
 
@@ -261,8 +258,8 @@ Finally, import the following library in the *Utilities.py* file:
 
 .. label:: end_Utilities_class
 
-Let us call the *update_neighbor_lists* and *update_cross_coefficients*
-methods from the *__init__* class:
+Let us call *update_neighbor_lists* and *update_cross_coefficients*
+from the *__init__* method:
 
 .. label:: start_InitializeSimulation_class
 
@@ -276,12 +273,11 @@ methods from the *__init__* class:
 
 .. label:: end_InitializeSimulation_class
         
-Test the code
+Test the Code
 -------------
 
-Let us test the *InitializeSimulation* class to make sure that it does what
-is expected, i.e. that it does create atoms that are located within the box
-boundaries along all 3 dimensions of space:
+Let us test the *InitializeSimulation* class to ensure that it behaves as
+expected, i.e., that it creates atoms located within the box boundaries:
 
 .. label:: start_test_3a_class
 
@@ -333,3 +329,6 @@ boundaries along all 3 dimensions of space:
         pytest.main(["-s", __file__])
 
 .. label:: end_test_3a_class
+
+The value of the cutoff, chosen as :math:`2.5 \sigma_{11}`, is relatively
+common for Lennard-Jones fluids and will often be the default choice for us.
diff --git a/docs/source/chapters/chapter4.rst b/docs/source/chapters/chapter4.rst
index 0746a69..fd3f35b 100644
--- a/docs/source/chapters/chapter4.rst
+++ b/docs/source/chapters/chapter4.rst
@@ -138,7 +138,7 @@ class:
             # Measure distance
             rij = self.compute_distance(self.atoms_positions[Ni],
                                         self.atoms_positions[neighbor_of_i],
-                                        self.box_size[:3])
+                                        self.box_size)
             # Measure potential using information about cross coefficients
             sigma_ij = self.sigma_ij_list[Ni]
             epsilon_ij = self.epsilon_ij_list[Ni]
@@ -158,11 +158,10 @@ class as well:
     def compute_distance(self,position_i, positions_j, box_size, only_norm = True):
         """
         Measure the distances between two particles.
-        The nan_to_num is crutial in 2D to avoid nan value along third dimension.
         # TOFIX: Move as function instead of a method?
         """
         rij_xyz = np.nan_to_num(np.remainder(position_i - positions_j
-                                + box_size[:3]/2.0, box_size) - box_size[:3]/2.0)
+                  + box_size[:3]/2.0, box_size[:3]) - box_size[:3]/2.0)
         if only_norm:
             return np.linalg.norm(rij_xyz, axis=1)
         else:
@@ -190,7 +189,7 @@ let us create a new method that is dedicated solely to measuring forces:
             # Measure distance
             rij, rij_xyz = self.compute_distance(self.atoms_positions[Ni],
                                         self.atoms_positions[neighbor_of_i],
-                                        self.box_size[:3], only_norm = False)
+                                        self.box_size, only_norm = False)
             # Measure force using information about cross coefficients
             sigma_ij = self.sigma_ij_list[Ni]
             epsilon_ij = self.epsilon_ij_list[Ni]       
diff --git a/docs/source/journal-article.bib b/docs/source/journal-article.bib
index 8869996..9461751 100644
--- a/docs/source/journal-article.bib
+++ b/docs/source/journal-article.bib
@@ -84,4 +84,15 @@ @article{         harris2020array
  doi           = {10.1038/s41586-020-2649-2},
  publisher     = {Springer Science and Business Media {LLC}},
  url           = {https://doi.org/10.1038/s41586-020-2649-2}
-}
\ No newline at end of file
+}
+
+@article{michaud2011mdanalysis,
+  title={MDAnalysis: a toolkit for the analysis of molecular dynamics simulations},
+  author={Michaud-Agrawal, Naveen and Denning, Elizabeth J and Woolf, Thomas B and Beckstein, Oliver},
+  journal={Journal of computational chemistry},
+  volume={32},
+  number={10},
+  pages={2319--2327},
+  year={2011},
+  publisher={Wiley Online Library}
+}

From 9762e9906ca15b380d2a122e6499a3e764aa9d39 Mon Sep 17 00:00:00 2001
From: Simon Gravelle <simon.gravelle@live.fr>
Date: Sun, 1 Sep 2024 14:07:50 +0200
Subject: [PATCH 11/13] fix grammar issue, chapter 4

---
 docs/source/chapters/chapter4.rst | 81 ++++++++++++++++---------------
 1 file changed, 42 insertions(+), 39 deletions(-)

diff --git a/docs/source/chapters/chapter4.rst b/docs/source/chapters/chapter4.rst
index fd3f35b..8118713 100644
--- a/docs/source/chapters/chapter4.rst
+++ b/docs/source/chapters/chapter4.rst
@@ -32,8 +32,10 @@ The steepest descent method is used for energy minimization, following these ste
 - 3) Move the atoms in the opposite direction of the maximum
      force to minimize the potential energy by a displacement step.
      The size of the step is adjusted iteratively based on the reduction in energy.
-- 4) Compute the new potential energy after the displacement, :math:`E_\text{pot}^\text{trial}`.
-- 5) Evaluate the change in energy: :math:`\Delta E = E_\text{pot}^\text{trial} - E_\text{pot}^\text{initial}`.
+- 4) Compute the new potential energy after the displacement,
+     :math:`E_\text{pot}^\text{trial}`.
+- 5) Evaluate the change in energy:
+     :math:`\Delta E = E_\text{pot}^\text{trial} - E_\text{pot}^\text{initial}`.
   
   - If :math:`\Delta E < 0`, the new configuration is accepted as it results in
     lower energy, and the step size is increased.
@@ -63,13 +65,8 @@ Let us fill the *__init__()* method:
             
 .. label:: end_MinimizeEnergy_class
 
-An important parameter is *maximum_steps*, which sets the maximum number
-of steps for the energy minimization process. The *displacement*
-parameter, with a default value of 0.01 Ångström, sets the initial atom
-displacement value.
-
-The *thermo_outputs* and *data_folder* parameters are used for printing data
-to files. These two parameters will be useful in the next chapter, :ref:`chapter5-label`.
+The parameter *maximum_steps* sets the maximum number
+of steps for the energy minimization process.
 
 Energy minimizer
 ----------------
@@ -89,8 +86,13 @@ following *run()* method to the *MinimizeEnergy* class:
             self.update_neighbor_lists() # Rebuild neighbor list, if necessary
             self.update_cross_coefficients() # Recalculate the cross coefficients, if necessary
             # Compute Epot/MaxF/force
-            init_Epot = self.compute_potential()
-            forces, init_MaxF = self.compute_force()
+            if hasattr(self, 'Epot') is False: # If self.Epot does not exists yet, calculate it
+                self.Epot = self.compute_potential()
+            if hasattr(self, 'MaxF') is False: # If self.MaxF does not exists yet, calculate it
+                forces = self.compute_force()
+                self.MaxF = np.max(np.abs(forces))
+            init_Epot = self.Epot
+            init_MaxF = self.MaxF
             # Save the current atom positions
             init_positions = copy.deepcopy(self.atoms_positions)
             # Move the atoms in the opposite direction of the maximum force
@@ -101,10 +103,9 @@ following *run()* method to the *MinimizeEnergy* class:
             # Keep the more favorable energy
             if trial_Epot < init_Epot: # accept new position
                 self.Epot = trial_Epot
-                # calculate the new max force and save it
-                forces, init_MaxF = self.compute_force()
+                forces = self.compute_force() # calculate the new max force and save it
                 self.MaxF = np.max(np.abs(forces))
-                self.wrap_in_box()  # Wrap atoms in the box, if necessary
+                self.wrap_in_box()  # Wrap atoms in the box
                 self.displacement *= 1.2 # Multiply the displacement by a factor 1.2
             else: # reject new position
                 self.Epot = init_Epot # Revert to old energy
@@ -113,17 +114,19 @@ following *run()* method to the *MinimizeEnergy* class:
 
 .. label:: end_MinimizeEnergy_class
 
-The displacement, which has an initial value of 0.01, is adjusted through energy
-minimization. When the trial is successful, its value is multiplied by 1.2. When
-the trial is rejected, its value is multiplied by 0.2.
+The displacement, which has an initial value of 0.01, is adjusted through
+energy minimization. When the trial is successful, its value is multiplied
+by 1.2. When the trial is rejected, its value is multiplied by 0.2. Thus,
+as the minimization progresses, the displacements that the particles
+perform increase.
 
-Compute_potential
------------------
+Compute the Potential
+---------------------
 
-Computing the potential energy of the system is central to the energy minimizer,
-as the value of the potential is used to decide if the trial is accepted or
-rejected. Add the following method called *compute_potential()*  to the *Utilities*
-class:
+Computing the potential energy of the system is central to the energy
+minimizer, as the value of the potential is used to decide if the trial is
+accepted or rejected. Add the following method called *compute_potential()*
+to the *Utilities* class:
 
 .. label:: start_Utilities_class
 
@@ -139,7 +142,7 @@ class:
             rij = self.compute_distance(self.atoms_positions[Ni],
                                         self.atoms_positions[neighbor_of_i],
                                         self.box_size)
-            # Measure potential using information about cross coefficients
+            # Measure potential using pre-calculated cross coefficients
             sigma_ij = self.sigma_ij_list[Ni]
             epsilon_ij = self.epsilon_ij_list[Ni]
             energy_potential += np.sum(potentials(epsilon_ij, sigma_ij, rij))
@@ -147,8 +150,8 @@ class:
     
 .. label:: end_Utilities_class
 
-Measuring the distance is an important step of computing the potential. Let us
-do it using a dedicated method. Add the following method to the *Utilities*
+Measuring the distance is a central step in computing the potential. Let us
+do this using a dedicated method. Add the following method to the *Utilities*
 class as well:
 
 .. label:: start_Utilities_class
@@ -158,7 +161,7 @@ class as well:
     def compute_distance(self,position_i, positions_j, box_size, only_norm = True):
         """
         Measure the distances between two particles.
-        # TOFIX: Move as function instead of a method?
+        # TOFIX: Move as a function instead of a method?
         """
         rij_xyz = np.nan_to_num(np.remainder(position_i - positions_j
                   + box_size[:3]/2.0, box_size[:3]) - box_size[:3]/2.0)
@@ -171,7 +174,9 @@ class as well:
 
 Finally, the energy minimization requires the computation of the minimum
 force in the system. Although not very different from the potential measurement,
-let us create a new method that is dedicated solely to measuring forces:
+let us create a new method that is dedicated solely to measuring forces.
+The force can be returned as a vector (default), or as a matrix (which will be
+useful for pressure calculation, see the :ref:`chapter7-label` chapter):
 
 .. label:: start_Utilities_class
 
@@ -202,23 +207,20 @@ let us create a new method that is dedicated solely to measuring forces:
                 # Add the contribution to the matrix
                 force_matrix[Ni][neighbor_of_i] += (fij_xyz*rij_xyz.T/rij).T
         if return_vector:
-            max_force = np.max(np.abs(force_vector))
-            return force_vector, max_force
+            return force_vector
         else:
             return force_matrix
     
 .. label:: end_Utilities_class
 
-Here, two types of outputs can
-be requested by the user: *force-vector*, and *force-matrix*. 
-The *force-matrix* option will be useful for pressure calculation, see
-:ref:`chapter7-label`.
+The force can be returned as a vector (default) or as a matrix. The latter
+will be useful for pressure calculation; see the :ref:`chapter7-label` chapter.
 
-Wrap in box
+Wrap in Box
 -----------
 
-Every time atoms are being displaced, one has to ensure that they remain in
-the box. This is done by the *wrap_in_box()* method that must be placed
+Every time atoms are displaced, one has to ensure that they remain within
+the box. This is done by the *wrap_in_box()* method, which must be placed
 within the *Utilities* class:
 
 .. label:: start_Utilities_class
@@ -296,5 +298,6 @@ typically negative.
 
 .. label:: end_test_4a_class
 
-For such as low density in particle, we can reasonably expect the energy to be always
-negative after 100 steps.
\ No newline at end of file
+For such a low particle density, we can reasonably expect that the potential
+energy will always be negative after 100 steps, and that the maximum force
+will be smaller than 10 (unitless).
\ No newline at end of file

From 5b3d37e2d8ea25299593a758d9c156ba148e38a8 Mon Sep 17 00:00:00 2001
From: Simon Gravelle <simon.gravelle@live.fr>
Date: Sun, 1 Sep 2024 15:07:34 +0200
Subject: [PATCH 12/13] improved grammar in chapters 5-7

---
 docs/source/chapters/chapter4.rst |   2 +-
 docs/source/chapters/chapter5.rst | 126 ++++++++++++++++++------------
 docs/source/chapters/chapter6.rst |  55 +++++++------
 docs/source/chapters/chapter7.rst |  54 ++++++++-----
 docs/source/journal-article.bib   |  22 ++++++
 5 files changed, 164 insertions(+), 95 deletions(-)

diff --git a/docs/source/chapters/chapter4.rst b/docs/source/chapters/chapter4.rst
index 8118713..e471092 100644
--- a/docs/source/chapters/chapter4.rst
+++ b/docs/source/chapters/chapter4.rst
@@ -300,4 +300,4 @@ typically negative.
 
 For such a low particle density, we can reasonably expect that the potential
 energy will always be negative after 100 steps, and that the maximum force
-will be smaller than 10 (unitless).
\ No newline at end of file
+will be smaller than 10 (unitless).
diff --git a/docs/source/chapters/chapter5.rst b/docs/source/chapters/chapter5.rst
index 12f46b7..3c04c34 100644
--- a/docs/source/chapters/chapter5.rst
+++ b/docs/source/chapters/chapter5.rst
@@ -4,39 +4,19 @@ Output the simulation state
 ===========================
 
 Here, the code is modified to allow us to follow the evolution of the system
-during the simulation. To do so, the *Output* class is modified.
-
-Update the MinimizeEnergy class
--------------------------------
-
-Let us start by calling two additional methods within the for loop of the
-*MinimizeEnergy* class, within the *run()* method.
-
-.. label:: start_MinimizeEnergy_class
-
-.. code-block:: python
-
-    def run(self):
-        (...)
-        for self.step in range(0, self.maximum_steps+1):
-            (...)
-                self.displacement *= 0.2 # Multiply the displacement by a factor 0.2
-            log_simulation_data(self)
-            update_dump_file(self, "dump.min.lammpstrj")
-
-.. label:: end_MinimizeEnergy_class
-
-The two methods named *update_log_minimize()* and *update_dump_file()*, are used
-to print the information in the terminal and in a LAMMPS-type data file, respectively.
-These two methods will be written in the following.
+during the simulation. To do so, two new files will be created: one dedicated
+to logging information, and the other to printing the atom positions to a
+file, thus allowing for their visualization with software like
+VMD :cite:`humphrey1996vmd` or Ovito :cite:`stukowski2009visualization`.
 
 Create logger
 -------------
 
-Let us create functions named *log_simulation_data* to a file named *logger.py*.
+Let us create a function named *log_simulation_data()* to a file named *logger.py*.
 With the logger, some output are being printed in a file, as well as in the terminal.
-The frequency of printing is set by *thermo_period*, see :ref:`chapter3-label`.
-All quantities are re-dimensionalized before getting outputed.
+The period of printing, which is a multiple of the simulation steps, will be set
+by a new parameter named *thermo_period* (integer). All quantities are 
+converted into *real* units before getting outputed (see :ref:`chapter2-label`).
 
 .. label:: start_logger_class
 
@@ -76,6 +56,7 @@ All quantities are re-dimensionalized before getting outputed.
         return logger
 
     def log_simulation_data(code):
+        # TOFIX: ceurrently, MaxF is returned dimensionless
 
         # Setup the logger with the folder name, overwriting the log if code.step is 0
         logger = setup_logger(code.data_folder, overwrite=(code.step == 0))
@@ -106,14 +87,22 @@ All quantities are re-dimensionalized before getting outputed.
 
 .. label:: end_logger_class
 
-Create dumper
+For simplicify, the *logger* uses the |logging| library, which provides a
+flexible framework for emitting log messages from Python programs. Depending on
+the value of *thermo_outputs*, different informations are printed, such as
+*step*, *Epot*, or/and *MaxF*.
+
+.. |logging| raw:: html
+
+   <a href="https://docs.python.org/3/library/logging.html" target="_blank">logging</a>
+
+Create Dumper
 -------------
 
-Let us create a function named *update_dump_file* to a file named 
-*dumper.py*. The dumper will print a *.lammpstrj file*, which contains the box
-dimensions and atom positions at every chosen frame (set by *dumping_period*,
-see :ref:`chapter3-label`). All quantities are dimensionalized before getting outputed, and the file follows
-a LAMMPS dump format, and can be read by molecular dynamics softwares like VMD.
+Let us create a function named *update_dump_file()* in a file named
+*dumper.py*. The dumper will print a *.lammpstrj* file, which contains
+the box dimensions and atom positions at every chosen frame (set by
+*dumping_period*). All quantities are dimensionalized before being output.
 
 .. label:: start_dumper_class
 
@@ -228,11 +217,35 @@ parameters are passed the InitializeSimulation method:
 
 .. label:: end_InitializeSimulation_class
 
+Update the MinimizeEnergy class
+-------------------------------
+
+As a final step, let us call two functions *log_simulation_data* and
+*update_dump_file* within the *for* loop of the
+*MinimizeEnergy* class, within the *run()* method:
+
+.. label:: start_MinimizeEnergy_class
+
+.. code-block:: python
+
+    def run(self):
+        (...)
+        for self.step in range(0, self.maximum_steps+1):
+            (...)
+                self.displacement *= 0.2 # Multiply the displacement by a factor 0.2
+            log_simulation_data(self)
+            update_dump_file(self, "dump.min.lammpstrj")
+
+.. label:: end_MinimizeEnergy_class
+
+The name *dump.min.lammpstrj* will be used when printing the dump file
+during energy minimization.
+
 Test the code
 -------------
 
-One can use a test similar as the previous ones. Let us ask out code to print
-information in the dump and the log files, and then let us assert the
+One can use a test similar as the previous ones. Let us ask our code to print
+information in the *dump* and the *log* files, and then let us assert that the
 files were indeed created without the *Outputs/* folder:
 
 .. label:: start_test_5a_class
@@ -276,8 +289,10 @@ files were indeed created without the *Outputs/* folder:
 
     # Test function using pytest
     def test_output_files():
-        assert os.path.exists("Outputs/dump.min.lammpstrj"), "Test failed: dump file was not created"
-        assert os.path.exists("Outputs/simulation.log"), "Test failed: log file was not created"
+        assert os.path.exists("Outputs/dump.min.lammpstrj"), \
+        "Test failed: the dump file was not created"
+        assert os.path.exists("Outputs/simulation.log"), \
+        "Test failed: the log file was not created"
         print("Test passed")
 
     # If the script is run directly, execute the tests
@@ -288,8 +303,9 @@ files were indeed created without the *Outputs/* folder:
 
 .. label:: end_test_5a_class
 
-I addition to the files getting created, information must be printed in the terminal
-during the simulation:
+In addition to making sure that both files were created, let us verify
+that the expected information was printed to the terminal during the
+simulation. The content of the *simulation.log* file should resemble:
 
 .. code-block:: bw
 
@@ -302,8 +318,8 @@ during the simulation:
 
 The data from the *simulation.log* can be used to generate plots using softwares
 line XmGrace, GnuPlot, or Python/Pyplot. For the later, one can use a simple data
-reader to import the data from *Outputs/simulation.log* into Python. Copy the
-following lines in a file named *reader.py*:
+reader to import the data from *simulation.log* into Python. Copy the
+following lines in a new file named *reader.py*:
 
 .. label:: start_reader_class
 
@@ -342,14 +358,28 @@ The *import_data* function from *reader.py* can simply be used as follows:
 
 .. label:: start_test_5b_class
 
+.. code-block:: python
+
     from reader import import_data
 
-    file_path = "Outputs/simulation.log"
-    header, data = import_data(file_path)
+    def test_file_not_empty():
+        # Import data from the file
+        file_path = "Outputs/simulation.log"
+        header, data = import_data(file_path)
+        # Check if the header and data are not empty
+        assert header, "Error, no header in simulation.log"
+        assert data, "Error, no data in simulation.log"
+        assert len(data) > 0, "Data should contain at least one row"
 
-    print(header)
-    for row in data:
-        print(row)
+        print(header)
+        for row in data:
+            print(row)
+
+    # If the script is run directly, execute the tests
+    if __name__ == "__main__":
+        import pytest
+        # Run pytest programmatically
+        pytest.main(["-s", __file__])
 
 .. label:: end_test_5b_class
 
@@ -362,4 +392,4 @@ Which must return:
     [25.0, -2.12, 1.22]
     [50.0, -2.19, 2.85]
     [75.0, -2.64, 0.99]
-    [100.0, -2.64, 0.51]
\ No newline at end of file
+    [100.0, -2.64, 0.51]
diff --git a/docs/source/chapters/chapter6.rst b/docs/source/chapters/chapter6.rst
index 6e4beec..6b7ea4d 100644
--- a/docs/source/chapters/chapter6.rst
+++ b/docs/source/chapters/chapter6.rst
@@ -47,13 +47,15 @@ Let us add a method named *monte_carlo_move* to the *MonteCarlo* class:
     def monte_carlo_move(self):
         """Monte Carlo move trial."""
         if self.displace_mc is not None: # only trigger if displace_mc was provided by the user
-            # If needed, recalculate neighbor/coeff lists
+            # When needed, recalculate neighbor/coeff lists
             self.update_neighbor_lists()
             self.update_cross_coefficients()
-            if hasattr(self, 'Epot') is False: # If self.Epot does not exists yet, calculate it
+            # If self.Epot does not exists yet, calculate it
+            # It should only be necessary when step = 0
+            if hasattr(self, 'Epot') is False:
                 self.Epot = self.compute_potential()
+            # Make a copy of the initial atoms positions and initial energy
             initial_Epot = self.Epot
-            # Make a copy of the initial atoms positions
             initial_positions = copy.deepcopy(self.atoms_positions)
             # Pick an atom id randomly
             atom_id = np.random.randint(np.sum(self.number_atoms))
@@ -61,7 +63,7 @@ Let us add a method named *monte_carlo_move* to the *MonteCarlo* class:
             # The maximum displacement is set by self.displace_mc
             move = (np.random.random(3)-0.5)*self.displace_mc 
             self.atoms_positions[atom_id] += move
-            # Measure the optential energy of the new configuration
+            # Measure the potential energy of the new configuration
             trial_Epot = self.compute_potential()
             # Evaluate whether the new configuration should be kept or not
             beta =  1/self.desired_temperature
@@ -77,13 +79,16 @@ Let us add a method named *monte_carlo_move* to the *MonteCarlo* class:
 
 .. label:: end_MonteCarlo_class
 
+The counters *successful_move* and *failed_move* are incremented with each
+successful and failed attempt, respectively.
+
 Parameters
 ----------
 
-The *monte_carlo_move* method requires a few parameters to be selected by the
-users, such as *displace_mc* (:math:`d_\text{mc}`), the maximum number of steps,
-and the desired temperature (:math:`T`). Let us add these parameters to the
-*__init__* method:
+The *monte_carlo_move* method requires a few parameters to be selected by
+the users, such as *displace_mc* (:math:`d_\text{mc}`, in Ångströms), the
+maximum number of steps, and the desired temperature (:math:`T`, in Kelvin).
+Let us add these parameters to the *__init__()* method:
 
 .. label:: start_MonteCarlo_class
 
@@ -106,11 +111,11 @@ and the desired temperature (:math:`T`). Let us add these parameters to the
 
 .. label:: end_MonteCarlo_class
 
-Run method
+Run Method
 ----------
 
-Finally, let us add a *run* method to the *MonteCarlo* class, that is used to
-perform a loop over the desired number of steps *maximum_steps*:
+Finally, let us add a *run* method to the *MonteCarlo* class that performs
+a loop over the desired number of steps, *maximum_steps*:
 
 .. label:: start_MonteCarlo_class
 
@@ -124,12 +129,10 @@ perform a loop over the desired number of steps *maximum_steps*:
 
 .. label:: end_MonteCarlo_class
 
-At each step, the *monte_carlo_move* method is called. The previously defined
-mthe *wrap_in_box* method is also called to ensure that
-the atoms remain inside the box, respectively.
-
-Let us call *update_log_md_mc* from the run method of the MonteCarlo class.
-Let us add a dump too:
+At each step, the *monte_carlo_move()* method is called. The previously
+defined *wrap_in_box* method is also called to ensure that the atoms remain
+inside the box. Additionally, let us call *log_simulation_data()* and
+*update_dump_file()*:
 
 .. label:: start_MonteCarlo_class
 
@@ -145,11 +148,11 @@ Let us add a dump too:
 
 .. label:: end_MonteCarlo_class
 
-Test the code
+Test the Code
 -------------
 
-One can use a similar test as previously. Let us use a displace distance of
-0.5 Angstrom, and make 1000 steps.
+Let us use a similar test as before. Set a displacement distance corresponding
+to a quarter of sigma, and perform a very small number of steps:
 
 .. label:: start_test_6a_class
 
@@ -194,14 +197,14 @@ One can use a similar test as previously. Let us use a displace distance of
         neighbor=20,
         displace_mc = displace_mc,
     )
-
-    # Run the Monte Carlo simulation
     mc.run()
 
     # Test function using pytest
     def test_output_files():
-        assert os.path.exists("Outputs/dump.mc.lammpstrj"), "Test failed: dump file was not created"
-        assert os.path.exists("Outputs/simulation.log"), "Test failed: log file was not created"
+        assert os.path.exists("Outputs/dump.mc.lammpstrj"), \
+        "Test failed: dump file was not created"
+        assert os.path.exists("Outputs/simulation.log"), \
+        "Test failed: log file was not created"
         print("Test passed")
 
     # If the script is run directly, execute the tests
@@ -213,5 +216,5 @@ One can use a similar test as previously. Let us use a displace distance of
 .. label:: end_test_6a_class
 
 The evolution of the potential energy as a function of the number of steps
-are written in the *simulation.log* file. The data can be used to plot
-the evolution of the system with time.
+is recorded in the *simulation.log* file. The data in *simulation.log* can
+be used to plot the evolution of the system over time.
diff --git a/docs/source/chapters/chapter7.rst b/docs/source/chapters/chapter7.rst
index d634e7a..0d0a8e8 100644
--- a/docs/source/chapters/chapter7.rst
+++ b/docs/source/chapters/chapter7.rst
@@ -1,36 +1,38 @@
 .. _chapter7-label:
 
-Pressure measurement
+Pressure Measurement
 ====================
 
-In order to extract the equation of state in our simulation, we need to measure the
-pressure of the system, :math:`p`. The pressure in a molecular simulation can be
-calculated from the interactions between particles. The pressure can be measured as
-the sum of the ideal contribution, :math:`p_\text{ideal} = N_\text{DOF} k_\text{B} T / V d`,
-which comes from the ideal gas law, and a Virial term which accounts for the
-pressure contribution from the forces between particles,
+To extract the equation of state in our simulation, we need to measure the
+pressure of the system, :math:`p`. The pressure in a molecular simulation can
+be calculated from the interactions between particles. The pressure can be
+measured as the sum of the ideal contribution,
+:math:`p_\text{ideal} = N_\text{DOF} k_\text{B} T / V d`, 
+which comes from the ideal gas law, and a Virial term that
+accounts for the pressure contribution from the forces between particles,
 :math:`p_\text{non_ideal} = \left< \sum_i r_i \cdot F_i \right> / V d`. The final
-expression reads:
+expression reads :cite:`frenkel2023understanding`:
 
-.. math:: 
+.. math::
 
     p = \dfrac{1}{V d} \left[ N_\text{DOF} k_\text{B} T +  \left< \sum_i r_i \cdot F_i \right> \right]
 
 :math:`N_\text{DOF}` is the number of degrees-of-freedom, which can be calculated
-from the number of particles, :math:`N`, and the dimension of the system, :math:`d`, as 
-:math:`N_\text{DOF} = d N - d` :cite:`frenkel2023understanding`.
+from the number of particles, :math:`N`, and the dimension of the system,
+:math:`d = 3`, as :math:`N_\text{DOF} = d N - d` :cite:`frenkel2023understanding`.
 
-The calculation of :math:`p_\text{ideal}` is straighforward. For Monte Carlo simulation,
-as atoms do not have temperature, the *imposed* temperature will be used instead.
-The calculation of :math:`p_\text{non_ideal}` requires the measurement of all the
-force and distance between the atoms. The calculation of the forces was already
-implemented in a previous chapter, but a new function that returns all the
-vector direction between atoms pairs will have to be written here.
+The calculation of :math:`p_\text{ideal}` is straightforward. For Monte Carlo
+simulation, as atoms do not have a temperature, the *imposed* temperature will
+be used as ":math:`T`". The calculation of :math:`p_\text{non_ideal}` requires the
+measurement of all the forces and distances between atoms. The calculation of
+the forces was already implemented in a previous chapter (see *compute_force*
+the :ref:`chapter4-label` chapter), but a new function that returns all the
+vector directions between atom pairs will need to be written here.
 
 Implement the Virial equation
 -----------------------------
 
-Let us add the following method to the *Measurements* class.
+Let us add the following method to the *Measurements* class:
 
 .. label:: start_Measurements_class
 
@@ -134,8 +136,10 @@ Let us test the outputed pressure.
 
     # Test function using pytest
     def test_output_files():
-        assert os.path.exists("Outputs/dump.mc.lammpstrj"), "Test failed: dump file was not created"
-        assert os.path.exists("Outputs/simulation.log"), "Test failed: log file was not created"
+        assert os.path.exists("Outputs/dump.mc.lammpstrj"), \
+        "Test failed: dump file was not created"
+        assert os.path.exists("Outputs/simulation.log"), \
+        "Test failed: log file was not created"
         print("Test passed")
 
     # If the script is run directly, execute the tests
@@ -145,3 +149,13 @@ Let us test the outputed pressure.
         pytest.main(["-s", __file__])
 
 .. label:: end_test_7a_class
+
+The pressure should be returned alongside the potential energy within *simulation.log*:
+
+.. code-block:: bw
+
+    step Epot press
+    0 134248.72 4608379.27
+    10 124905.76 4287242.75
+    20 124858.91 4285584.40
+    (...)
\ No newline at end of file
diff --git a/docs/source/journal-article.bib b/docs/source/journal-article.bib
index 9461751..6616f3b 100644
--- a/docs/source/journal-article.bib
+++ b/docs/source/journal-article.bib
@@ -96,3 +96,25 @@ @article{michaud2011mdanalysis
   year={2011},
   publisher={Wiley Online Library}
 }
+
+@article{humphrey1996vmd,
+  title={VMD: visual molecular dynamics},
+  author={Humphrey, William and Dalke, Andrew and Schulten, Klaus},
+  journal={Journal of molecular graphics},
+  volume={14},
+  number={1},
+  pages={33--38},
+  year={1996},
+  publisher={Elsevier}
+}
+
+@article{stukowski2009visualization,
+  title={Visualization and analysis of atomistic simulation data with OVITO--the Open Visualization Tool},
+  author={Stukowski, Alexander},
+  journal={Modelling and simulation in materials science and engineering},
+  volume={18},
+  number={1},
+  pages={015012},
+  year={2009},
+  publisher={IOP Publishing}
+}

From 8951d8bae00e65d23bcb573e6b2a541d50973323 Mon Sep 17 00:00:00 2001
From: Simon Gravelle <simon.gravelle@live.fr>
Date: Sun, 1 Sep 2024 15:20:37 +0200
Subject: [PATCH 13/13] fix grammar in chapter 8

---
 docs/source/chapters/chapter8.rst | 48 +++++++++++++++----------------
 1 file changed, 23 insertions(+), 25 deletions(-)

diff --git a/docs/source/chapters/chapter8.rst b/docs/source/chapters/chapter8.rst
index e4c9b8c..920c899 100644
--- a/docs/source/chapters/chapter8.rst
+++ b/docs/source/chapters/chapter8.rst
@@ -16,12 +16,6 @@ simulation resemble the Monte Carlo move from :ref:`chapter6-label`:
 - 4) We then decide to keep or reject the move by calculating
   the difference in energy between the trial and the initial configurations:
   :math:`\Delta E = E_\text{pot}^\text{trial} - E_\text{pot}^\text{initial}`.
-  
-  - If :math:`\Delta E < 0`, then the move is automatically accepted. 
-  - If :math:`\Delta E > 0`, then the move is accepted with a probability given
-    by the Boltzmann factor :math:`\exp{- \beta \Delta E}`, where
-    :math:`\beta = 1 / k_\text{B} T` and :math:`T` is the imposed temperature.
-
 - 5) Steps 1-4 are repeated a large number of times, generating a broad range of
      possible configurations.
 
@@ -38,7 +32,7 @@ Let us add the following method called *monte_carlo_exchange* to the *MonteCarlo
         # The first step is to make a copy of the previous state
         # Since atoms numbers are evolving, its also important to store the
         # neighbor, sigma, and epsilon lists
-        self.Epot = self.compute_potential()
+        self.Epot = self.compute_potential() # TOFIX: not necessary every time
         initial_positions = copy.deepcopy(self.atoms_positions)
         initial_number_atoms = copy.deepcopy(self.number_atoms)
         initial_neighbor_lists = copy.deepcopy(self.neighbor_lists)
@@ -72,15 +66,15 @@ Let us add the following method called *monte_carlo_exchange* to the *MonteCarlo
 .. label:: end_MonteCarlo_class
 
 First, the potential energy is calculated, and the current state of the
-simulations, such as positions and atom numbers, are saved. Then, an insertion
-try or a deletion trial is made, each with a probability of 0.5. The
-*monte_carlo_insert* and *monte_carlo_delete*, that are implemented here below,
-are both calculting the *acceptation_probability*. A random number is again selected,
-and if that number is larger than the acceptation probability, then the system
-revert to its previous position. If the random number is lower than the acceptation
-probability, nothing appends, which means that the last trial is accepted.
+simulation, such as positions and atom numbers, is saved. Then, an insertion
+trial or a deletion trial is made, each with a probability of 0.5. The
+*monte_carlo_insert* and *monte_carlo_delete* methods, which are implemented
+below, both calculate the *acceptance_probability*. A random number is selected,
+and if that number is larger than the acceptance probability, the system reverts
+to its previous position. If the random number is lower than the acceptance
+probability, nothing happens, meaning that the last trial is accepted.
 
-Then, let us write the *monte_carlo_insert* method:
+Then, let us write the *monte_carlo_insert()* method:
 
 .. label:: start_MonteCarlo_class
 
@@ -96,7 +90,6 @@ Then, let us write the *monte_carlo_insert* method:
         self.atoms_positions = np.vstack([self.atoms_positions[:shift_id],
                                         new_atom_position,
                                         self.atoms_positions[shift_id:]])
-        self.total_number_atoms = np.sum(self.number_atoms)
         self.update_neighbor_lists()
         self.identify_atom_properties()
         self.update_cross_coefficients()
@@ -111,7 +104,10 @@ Then, let us write the *monte_carlo_insert* method:
 
 .. label:: end_MonteCarlo_class
 
-as well as the *monte_carlo_delete* method:
+After trying to insert a new particle, neighbor lists and cross coefficients
+must be re-evaluated. Then, the acceptance probability is calculated.
+
+Let us add the very similar *monte_carlo_delete()* method:
 
 .. label:: start_MonteCarlo_class
 
@@ -192,7 +188,7 @@ Let us also normalised the "desired_mu":
 
 .. label:: end_MonteCarlo_class
 
-Finally, the *monte_carlo_insert_delete()* method must be included in the run:
+Finally, the *monte_carlo_exchange()* method must be included in the run:
 
 .. label:: start_MonteCarlo_class
 
@@ -207,7 +203,7 @@ Finally, the *monte_carlo_insert_delete()* method must be included in the run:
 
 .. label:: end_MonteCarlo_class
 
-We need to calculate Lambda:
+We need to calculate :math:`\Lambda`:
 
 .. label:: start_MonteCarlo_class
 
@@ -228,6 +224,7 @@ To output the density, let us add the following method to the *Measurements* cla
 
     def calculate_density(self):
         """Calculate the mass density."""
+        # TOFIX: not used yet
         volume = np.prod(self.box_size[:3])  # Unitless
         total_mass = np.sum(self.atoms_mass)  # Unitless
         return total_mass/volume  # Unitless
@@ -237,8 +234,8 @@ To output the density, let us add the following method to the *Measurements* cla
 Test the code
 -------------
 
-One can use a similar test as previously. Let us use a displace distance of
-0.5 Angstrom, and make 1000 steps.
+One can use a similar test as previously, but with an imposed chemical
+potential *desired_mu*:
 
 .. label:: start_test_8a_class
 
@@ -287,8 +284,10 @@ One can use a similar test as previously. Let us use a displace distance of
 
     # Test function using pytest
     def test_output_files():
-        assert os.path.exists("Outputs/dump.mc.lammpstrj"), "Test failed: dump file was not created"
-        assert os.path.exists("Outputs/simulation.log"), "Test failed: log file was not created"
+        assert os.path.exists("Outputs/dump.mc.lammpstrj"), \
+        "Test failed: dump file was not created"
+        assert os.path.exists("Outputs/simulation.log"), \
+        "Test failed: log file was not created"
         print("Test passed")
 
     # If the script is run directly, execute the tests
@@ -300,5 +299,4 @@ One can use a similar test as previously. Let us use a displace distance of
 .. label:: end_test_8a_class
 
 The evolution of the potential energy as a function of the
-number of steps are written in the *Outputs/Epot.dat* file
-and can be plotted.
+number of steps are written in the *Outputs/Epot.dat*.