started to simplify potential functions

simongravelle · simongravelle · commit 23b110ee95bb · 2024-08-29T23:32:15.000+02:00
diff --git a/docs/source/chapters/chapter2.rst b/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
 -------------
 
diff --git a/docs/source/chapters/chapter4.rst b/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_vector_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_vector_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
@@ -274,42 +270,92 @@ 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
+
+.. 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
+
+.. label:: start_Utilities_class
+
+.. code-block:: python
+
+    def compute_vector_force(self):
+        force_vector = np.zeros((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)
+            # 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 
+        max_force = np.max(np.abs(force_vector))
+        return force_vector, max_force
+    
+.. label:: end_Utilities_class
+
+.. label:: start_Utilities_class
+
+.. code-block:: python
+
+    def compute_vector_matrix(self):
+        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)
+            # Add the contribution to the force matrix
+            force_matrix[Ni][neighbor_of_i] += (fij_xyz*rij_xyz.T/rij).T
+        return force_matrix
 
 .. label:: end_Utilities_class
 
diff --git 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, 3, 4]: #, 5, 6, 7]:
     # for each chapter, create the corresponding code
     RST_EXISTS, created_tests, folder = sphinx_to_python(path_to_docs, chapter_id)
     if RST_EXISTS:
