Merge pull request #7 from mdcourse/make-code-2D-compatible

start adapting the code to 2D

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/chapter3.rst

@@ -35,7 +35,8 @@ Let us improve the previously created *InitializeSimulation* class:
                     ):
             super().__init__(*args, **kwargs)
             self.box_dimensions = box_dimensions
-            self.dimensions = len(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.initial_positions = initial_positions
             self.thermo_period = thermo_period
@@ -137,9 +138,14 @@ method to the *InitializeSimulation* class:
 .. code-block:: python
 
     def define_box(self):
-        box_boundaries = np.zeros((self.dimensions, 2))
-        for dim, L in zip(range(self.dimensions), self.box_dimensions):
+        """Define the simulation box.
+        For 2D simulations, the third dimensions only contains 0.
+        """
+        box_boundaries = np.zeros((3, 2))
+        dim = 0
+        for L in 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])
@@ -183,9 +189,8 @@ 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,
-                                        self.dimensions))
-            for dim in np.arange(self.dimensions):
+            atoms_positions = np.zeros((self.total_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)
                 atoms_positions[:, dim] = random_pos*diff_box-diff_box/2

docs/source/chapters/chapter4.rst

@@ -291,13 +291,16 @@ class.
             sigma_ij = self.sigma_ij_list[Ni]
             epsilon_ij = self.epsilon_ij_list[Ni]
             # Measure distances
-            rij_xyz = (np.remainder(position_i - positions_j + half_box_size, box_size) - half_box_size)
+            # Measure distances
+            # The nan_to_num is crutial in 2D to avoid nan value along third dimension
+            rij_xyz = np.nan_to_num(np.remainder(position_i - positions_j
+                                                 + half_box_size, box_size) - half_box_size)
             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 

docs/source/chapters/chapter6.rst

@@ -57,7 +57,10 @@ Let us add a method named *monte_carlo_move* to the *MonteCarlo* class:
             atom_id = np.random.randint(self.total_number_atoms)
             # Move the chosen atom in a random direction
             # The maximum displacement is set by self.displace_mc
-            move = (np.random.random(self.dimensions)-0.5)*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)
             self.atoms_positions[atom_id] += move
             # Measure the optential energy of the new configuration
             trial_Epot = self.compute_potential(output="potential")

docs/source/chapters/chapter7.rst

@@ -39,16 +39,16 @@ Let us add the following method to the *Utilities* class.
     def calculate_pressure(self):
         """Evaluate p based on the Virial equation (Eq. 4.4.2 in Frenkel-Smit,
         Understanding molecular simulation: from algorithms to applications, 2002)"""
-        # Ideal contribution
+        # Compute the ideal contribution
         Ndof = self.dimensions*self.total_number_atoms-self.dimensions    
-        volume = np.prod(self.box_size[:3])
+        volume = np.prod(self.box_size[:self.dimensions])
         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)
-        # Non-ideal contribution
+        # Compute the non-ideal contribution
         distances_forces = np.sum(self.compute_potential(output="force-matrix")*self.evaluate_rij_matrix())
         p_nonideal = distances_forces/(volume*self.dimensions)
         # Final pressure
@@ -73,7 +73,8 @@ To evaluate all the vectors between all the particles, let us also add the
             position_i = self.atoms_positions[Ni]
             rij_xyz = (np.remainder(position_i - positions_j + half_box_size, box_size) - half_box_size)
             rij_matrix[Ni] = rij_xyz
-        return rij_matrix
+        # use nan_to_num to avoid "nan" values in 2D
+        return np.nan_to_num(rij_matrix)
 
 .. label:: end_Utilities_class