start adapting the code to 3D

Author: simongravelle

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,11 +138,16 @@ 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_size = np.diff(self.box_boundaries).reshape(self.dimension)
         box_geometry = np.array([90, 90, 90])
         self.box_size = np.array(box_size.tolist()+box_geometry.tolist())
 
@@ -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,7 +291,10 @@ 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":

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")