Merge pull request #8 from mdcourse/simplify_compute_potential

Simplify potential functions

Author: simongravelle

docs/source/chapters/chapter2.rst

@@ -293,86 +293,6 @@ Let us call the *identify_atom_properties* from the *__init__()* method:
 
 .. label:: end_Prepare_class
 
-..
-    Calculate cross coefficients
-    ----------------------------
-
-    Let us calculate all cross coefficients that are required to calculate the interactions
-    between atom :math:`i` and atom :math:`j`. From the example described previously,
-    where:
-
-    .. math::
-
-        \text{atoms_sigma} = [\sigma_{11}, \sigma_{11}, \sigma_{22}, \sigma_{22}, \sigma_{22}]
-
-    one expects all direct and cross coefficients to be:
-
-    .. math::
-
-        \text{array_sigma_ij} = \begin{bmatrix}
-                \sigma_{11} & \sigma_{11} & \sigma_{12} & \sigma_{12} & \sigma_{12} \\
-                \sigma_{11} & \sigma_{11} & \sigma_{12} & \sigma_{12} & \sigma_{12} \\
-                \sigma_{12} & \sigma_{12} & \sigma_{22} & \sigma_{22} & \sigma_{22} \\
-                \sigma_{12} & \sigma_{12} & \sigma_{22} & \sigma_{22} & \sigma_{22} \\
-                \sigma_{12} & \sigma_{12} & \sigma_{22} & \sigma_{22} & \sigma_{22} \\
-            \end{bmatrix}
-
-
-    The matrix is symmetric, so the coefficients in the bottom left corner are 
-    the same as the coefficient in the top right corner. The first value in the top left corner of the matrix,
-    :math:`\sigma_{11}`, indicates that the :math:`\sigma` value for the interaction
-    between the atom 1 and itself is is :math:`\sigma_{11}`. A similar matrix can
-    be written for epsilon_sigma_ij.
-
-    Here, the values of the cross coefficients :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
-
-    Create the following method called *calculate_cross_coefficients* within the 
-    *Prepare* class:
-
-    . . label:: start_Prepare_class
-
-    .. code-block:: python
-
-        def calculate_cross_coefficients(self):
-            """Calculate all the cross-coefficients for the LJ interations."""
-            self.identify_atom_properties() # TOFIX: this was left because of GCMC. Remove? Move?
-            matrix_epsilon_ij = []
-            matrix_sigma_ij = []
-            for i in range(self.total_number_atoms):
-                matrix_epsilon_ij.append([])
-                matrix_sigma_ij.append([])
-                for j in range(self.total_number_atoms):
-                    matrix_epsilon_ij[-1].append((self.atoms_epsilon[i]+self.atoms_epsilon[j])/2)
-                    matrix_sigma_ij[-1].append((self.atoms_sigma[i]+self.atoms_sigma[j])/2)
-            self.matrix_sigma_ij = matrix_sigma_ij
-            self.matrix_epsilon_ij = matrix_epsilon_ij
-
-    . . label:: end_Prepare_class
-
-    After calling for the *identify_atom_properties()* method, a double loop
-    is performed over all direct coefficients, and the cross coefficients
-    are stored within *matrix_sigma_ij* and *matrix_epsilon_ij*, two matrices.
-
-    Finally, let us call the *calculate_cross_coefficients* method from the
-    *__init__()* method.
-
-    . . label:: start_Prepare_class
-
-    .. code-block:: python
-
-        def __init__(self,
-            (...)
-            self.identify_atom_properties()
-            self.calculate_cross_coefficients()
-
-    . . label:: end_Prepare_class
-
 Test the code
 -------------
 

docs/source/chapters/chapter4.rst

@@ -146,25 +146,21 @@ following *run()* method to the *MinimizeEnergy* class:
             # First, meevaluate the initial energy and max force
             self.update_neighbor_lists() # Rebuild neighbor list, if necessary
             self.update_cross_coefficients() # Recalculate the cross coefficients, if necessary
-            try: # try using the last saved Epot and MaxF, if it exists
-                init_Epot = self.Epot
-                init_MaxF = self.MaxF
-            except: # If Epot/MaxF do not exists yet, calculate them both
-                init_Epot = self.compute_potential(output="potential")
-                forces = self.compute_potential(output="force-vector")
-                init_MaxF = np.max(np.abs(forces))
+            if self.step == 0: # At the first step, Epot/MaxF do not exists yet, calculate them both
+                init_Epot = self.compute_potential()
+                forces, init_MaxF = self.compute_force()
             # Save the current atom positions
             init_positions = copy.deepcopy(self.atoms_positions)
             # Move the atoms in the opposite direction of the maximum force
             self.atoms_positions = self.atoms_positions \
                 + forces/init_MaxF*self.displacement
             # Recalculate the energy
-            trial_Epot = self.compute_potential(output="potential")
+            trial_Epot = self.compute_potential()
             # 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 = self.compute_potential(output="force-vector")
+                forces, init_MaxF = self.compute_force()
                 self.MaxF = np.max(np.abs(forces))
                 self.wrap_in_box()  # Wrap atoms in the box, if necessary
                 self.displacement *= 1.2 # Multiply the displacement by a factor 1.2
@@ -268,56 +264,99 @@ Compute_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.
+class:
 
 .. label:: start_Utilities_class
 
 .. code-block:: python
 
-    def compute_potential(self, output):
-        if output == "force-vector":
-            forces = np.zeros((self.total_number_atoms,3))
-        elif output == "force-matrix":
-            forces = np.zeros((self.total_number_atoms,self.total_number_atoms,3))
+    def compute_potential(self):
+        """Compute the potential energy by summing up all pair contributions."""
         energy_potential = 0
-        box_size = self.box_size[:3]
-        half_box_size = self.box_size[:3]/2.0
         for Ni in np.arange(self.total_number_atoms-1):
-            # Read information about atom i
-            position_i = self.atoms_positions[Ni]
+            # Read neighbor list
             neighbor_of_i = self.neighbor_lists[Ni]
-            # Read information about neighbors j and cross coefficient
-            positions_j = self.atoms_positions[neighbor_of_i]
+            # Measure distance
+            rij = self.compute_distance(self.atoms_positions[Ni],
+                                        self.atoms_positions[neighbor_of_i],
+                                        self.box_size[:3])
+            # Measure potential using information about cross coefficients
             sigma_ij = self.sigma_ij_list[Ni]
             epsilon_ij = self.epsilon_ij_list[Ni]
-            # Measure distances
-            # 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(potentials(self.potential_type, epsilon_ij, sigma_ij, rij))
-            else:
-                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 
-                elif output == "force-matrix":
-                    forces[Ni][neighbor_of_i] += (derivative_potential*rij_xyz.T/rij).T
-        if output=="potential":
-            return energy_potential
-        elif (output == "force-vector") | (output == "force-matrix"):
-            return forces
+            energy_potential += np.sum(potentials(self.potential_type,
+                                                  epsilon_ij, sigma_ij, rij))
+        return energy_potential
+    
+.. 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*
+class as well:
+
+.. label:: start_Utilities_class
+
+.. code-block:: python
+
+    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)
+        if only_norm:
+            return np.linalg.norm(rij_xyz, axis=1)
+        else:
+            return np.linalg.norm(rij_xyz, axis=1), rij_xyz
+
+.. label:: end_Utilities_class
 
+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:
+
+.. label:: start_Utilities_class
+
+.. code-block:: python
+
+    def compute_force(self, return_vector = True):
+        if return_vector: # return a N-size vector
+            force_vector = np.zeros((self.total_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):
+            # Read neighbor list
+            neighbor_of_i = self.neighbor_lists[Ni]
+            # 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)
+            # 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)
+            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)
+                force_vector[neighbor_of_i] -= (fij_xyz*rij_xyz.T/rij).T 
+            else:
+                # 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
+        else:
+            return force_matrix
+    
 .. label:: end_Utilities_class
 
-Here, the method is a little bit complicated, because three types of outputs can
-be requested by the user: *force-vector*, *force-matrix*, and *potential*. The last
-one, *potential*, simply returns the value of the potential energy for the entire system.
-If *force-vector* or *force-matrix* are selected instead, then the individual forces
-between atoms are returned.
+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`.
 
 Wrap in box
 -----------

docs/source/chapters/chapter5.rst

@@ -84,7 +84,7 @@ All quantities are re-dimensionalized before getting outputed.
         if code.thermo_period is not None:
             if code.step % code.thermo_period == 0:
                 if code.step == 0:
-                    Epot = code.compute_potential(output="potential") \
+                    Epot = code.compute_potential() \
                         * code.reference_energy  # kcal/mol
                 else:
                     Epot = code.Epot * code.reference_energy  # kcal/mol

docs/source/chapters/chapter6.rst

@@ -50,7 +50,7 @@ Let us add a method named *monte_carlo_move* to the *MonteCarlo* class:
             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(output="potential")
+                initial_Epot = self.compute_potential()
             # Make a copy of the initial atoms positions
             initial_positions = copy.deepcopy(self.atoms_positions)
             # Pick an atom id randomly
@@ -63,7 +63,7 @@ Let us add a method named *monte_carlo_move* to the *MonteCarlo* class:
                 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")
+            trial_Epot = self.compute_potential()
             # Evaluate whether the new configuration should be kept or not
             beta =  1/self.desired_temperature
             delta_E = trial_Epot-initial_Epot

docs/source/chapters/chapter7.rst

@@ -49,7 +49,7 @@ Let us add the following method to the *Utilities* class.
             temperature = self.desired_temperature # for MC, simply use the desired temperature
         p_ideal = Ndof*temperature/(volume*self.dimensions)
         # Compute the non-ideal contribution
-        distances_forces = np.sum(self.compute_potential(output="force-matrix")*self.evaluate_rij_matrix())
+        distances_forces = np.sum(self.compute_force(return_vector = False)*self.evaluate_rij_matrix())
         p_nonideal = distances_forces/(volume*self.dimensions)
         # Final pressure
         self.pressure = p_ideal+p_nonideal