implemented hard sphere potential

Author: simongravelle

docs/source/chapters/chapter1.rst

@@ -66,33 +66,51 @@ 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 named *LJ_potential* for the Lennard-Jones
-potential (LJ). Copy the following into *potentials.py*:
+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.
+
+Copy the following lines into *potentials.py*:
 
 .. label:: start_potentials_class
 
 .. code-block:: python
 
-    def LJ_potential(epsilon, sigma, r, derivative = False):
-        if derivative:
-            return 48*epsilon*((sigma/r)**12-0.5*(sigma/r)**6)/r
+    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)
+        elif potential_type == "Hard-Sphere":
+            if derivative:
+                raise ValueError("Derivative is not defined for Hard-Sphere potential.")
+            else:
+                return 1000 if r < sigma else 0
         else:
-            return 4*epsilon*((sigma/r)**12-(sigma/r)**6)
+            raise ValueError(f"Unknown potential type: {potential_type}")
 
 .. label:: end_potentials_class
 
-Depending on the value of the optional argument *derivative*, which can be
-either *False* or *True*, the *LJ_potential* function will return the force:
+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:
 
 .. 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
 ------------------
@@ -124,7 +142,7 @@ copy the following lines:
 
 .. code-block:: python
 
-    from potentials import LJ_potential
+    from potentials import potentials
 
 
     class Utilities:

docs/source/chapters/chapter2.rst

@@ -99,12 +99,14 @@ Modify the *Prepare* class as follows:
                     epsilon=[0.1],  # List - Kcal/mol
                     sigma=[3],  # List - Angstrom
                     atom_mass=[10],  # List - g/mol
+                    potential_type="Lennard-Jones",
                     *args,
                     **kwargs):
             self.number_atoms = number_atoms
             self.epsilon = epsilon
             self.sigma = sigma
             self.atom_mass = atom_mass
+            self.potential_type = potential_type
             super().__init__(*args, **kwargs)
 
 .. label:: end_Prepare_class
@@ -114,7 +116,10 @@ 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.
 
-All four parameters are assigned to *self*, allowing other methods to access
+The type of potential is also specified, with Lennard-Jones being closen as
+the default option.
+
+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
 of positional and keyword arguments, respectively.
 
@@ -130,7 +135,7 @@ unit system to the *LJ* unit system:
 .. code-block:: python
 
     def calculate_LJunits_prefactors(self):
-        """Calculate the Lennard-Jones units prefacors."""
+        """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

docs/source/chapters/chapter4.rst

@@ -298,9 +298,9 @@ class.
             rij = np.linalg.norm(rij_xyz, axis=1)
             # Measure potential
             if output == "potential":
-                energy_potential += np.sum(LJ_potential(epsilon_ij, sigma_ij, rij))
+                energy_potential += np.sum(potentials(self.potential_type, epsilon_ij, sigma_ij, rij))
             else:
-                derivative_potential = LJ_potential(epsilon_ij, sigma_ij, rij, derivative = True)
+                derivative_potential = potentials(self.potential_type, epsilon_ij, sigma_ij, rij, derivative = True)
                 if output == "force-vector":
                     forces[Ni] += np.sum((derivative_potential*rij_xyz.T/rij).T, axis=0)
                     forces[neighbor_of_i] -= (derivative_potential*rij_xyz.T/rij).T