Minimize the energy¶

Now that the code for placing the atoms within the box has been written, let us proceed to write the code for performing energy minimization of the system. This step helps ensure that there is no unreasonable overlapping between the atoms.

The steepest descent method is used for energy minimization, following these steps:

1. Start with an initial configuration and measure the initial potential energy, \(E_\text{pot}^\text{initial}\).
1. Calculate the gradient of the potential energy with respect to atomic positions to determine the direction of the steepest ascent in energy. The magnitude of this gradient is used to compute the maximum force acting on the atoms. This maximum force indicates the direction in which the energy decreases most rapidly.
1. Move the atoms in the opposite direction of the maximum force to minimize the potential energy by a displacement step. The size of the step is adjusted iteratively based on the reduction in energy.
1. Compute the new potential energy after the displacement, \(E_\text{pot}^\text{trial}\).
1. Evaluate the change in energy: \(\Delta E = E_\text{pot}^\text{trial} - E_\text{pot}^\text{initial}\).
- If \(\Delta E < 0\), the new configuration is accepted as it results in lower energy, and the step size is increased.
- If \(\Delta E \geq 0\), the new configuration is rejected, and the step size is decreased.

The process is repeated until the maximum number of steps is reached. The goal is to iteratively reduce the potential energy and reach a stable, minimized energy state.

Prepare the minimization¶

Let us fill the __init__() method:

class MinimizeEnergy(Measurements):
    def __init__(self,
                maximum_steps,
                *args,
                **kwargs):
        self.maximum_steps = maximum_steps
        super().__init__(*args, **kwargs)

The parameter maximum_steps sets the maximum number of steps for the energy minimization process.

Energy minimizer¶

Let us implement the energy minimized described at the top of this page. Add the following run() method to the MinimizeEnergy class:

def run(self):
    self.displacement = 0.01 # pick a random initial displacement (dimentionless)
    # *step* loops for 0 to *maximum_steps*+1
    for self.step in range(0, self.maximum_steps+1):
        # First, evaluate 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
        # Compute Epot/MaxF/force
        if hasattr(self, 'Epot') is False: # If self.Epot does not exist yet, calculate it
            self.Epot = self.compute_potential()
        if hasattr(self, 'MaxF') is False: # If self.MaxF does not exist yet, calculate it
            forces = self.compute_force()
            self.MaxF = np.max(np.abs(forces))
        init_Epot = self.Epot
        init_MaxF = self.MaxF
        # 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()
        # Keep the more favorable energy
        if trial_Epot < init_Epot: # accept new position
            self.Epot = trial_Epot
            forces = self.compute_force() # calculate the new max force and save it
            self.MaxF = np.max(np.abs(forces))
            self.wrap_in_box()  # Wrap atoms in the box
            self.displacement *= 1.2 # Multiply the displacement by a factor 1.2
        else: # reject new position
            self.Epot = init_Epot # Revert to old energy
            self.atoms_positions = init_positions # Revert to old positions
            self.displacement *= 0.2 # Multiply the displacement by a factor 0.2

The displacement, which has an initial value of 0.01, is adjusted through energy minimization. When the trial is successful, its value is multiplied by 1.2. When the trial is rejected, its value is multiplied by 0.2. Thus, as the minimization progresses, the displacements that the particles perform increase.

Compute the 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:

def compute_potential(self):
    """Compute the potential energy by summing up all pair contributions."""
    energy_potential = 0
    for Ni in np.arange(np.sum(self.number_atoms)-1):
        # Read neighbor list
        neighbor_of_i = self.neighbor_lists[Ni]
        # Measure distance
        rij = self.compute_distance(self.atoms_positions[Ni],
                                    self.atoms_positions[neighbor_of_i],
                                    self.box_size)
        # Measure potential using pre-calculated cross coefficients
        sigma_ij = self.sigma_ij_list[Ni]
        epsilon_ij = self.epsilon_ij_list[Ni]
        energy_potential += np.sum(potentials(epsilon_ij, sigma_ij, rij))
    return energy_potential

Measuring the distance is a central step in computing the potential. Let us do this using a dedicated method. Add the following method to the Utilities class as well:

def compute_distance(self,position_i, positions_j, box_size, only_norm = True):
    """
    Measure the distances between two particles.
    # TOFIX: Move as a function instead of a method?
    """
    rij_xyz = np.nan_to_num(np.remainder(position_i - positions_j
              + box_size[:3]/2.0, box_size[:3]) - 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

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. The force can be returned as a vector (default), or as a matrix (which will be useful for pressure calculation, see the Pressure Measurement chapter):

def compute_force(self, return_vector = True):
    if return_vector: # return a N-size vector
        force_vector = np.zeros((np.sum(self.number_atoms),3))
    else: # return a N x N matrix
        force_matrix = np.zeros((np.sum(self.number_atoms),
                                np.sum(self.number_atoms),3))
    for Ni in np.arange(np.sum(self.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, 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(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:
        return force_vector
    else:
        return force_matrix

The force can be returned as a vector (default) or as a matrix. The latter will be useful for pressure calculation; see the Pressure Measurement chapter.

Wrap in Box¶

Every time atoms are displaced, one has to ensure that they remain within the box. This is done by the wrap_in_box() method, which must be placed within the Utilities class:

def wrap_in_box(self):
    for dim in np.arange(3):
        out_ids = self.atoms_positions[:, dim] \
            > self.box_boundaries[dim][1]
        self.atoms_positions[:, dim][out_ids] \
            -= np.diff(self.box_boundaries[dim])[0]
        out_ids = self.atoms_positions[:, dim] \
            < self.box_boundaries[dim][0]
        self.atoms_positions[:, dim][out_ids] \
            += np.diff(self.box_boundaries[dim])[0]

Test the code¶

Let us test the MinimizeEnergy class to make sure that it does what is expected, i.e. that it leads to a potential energy that is small, and typically negative.

from MinimizeEnergy import MinimizeEnergy
from pint import UnitRegistry
ureg = UnitRegistry()

# Define atom number of each group
nmb_1, nmb_2= [2, 3]
# Define LJ parameters (sigma)
sig_1, sig_2 = [3, 4]*ureg.angstrom
# Define LJ parameters (epsilon)
eps_1, eps_2 = [0.2, 0.4]*ureg.kcal/ureg.mol
# Define atom mass
mss_1, mss_2 = [10, 20]*ureg.gram/ureg.mol
# Define box size
L = 20*ureg.angstrom
# Define a cut off
rc = 2.5*sig_1

# Initialize the prepare object
minimizer = MinimizeEnergy(
    ureg = ureg,
    maximum_steps=100,
    number_atoms=[nmb_1, nmb_2],
    epsilon=[eps_1, eps_2], # kcal/mol
    sigma=[sig_1, sig_2], # A
    atom_mass=[mss_1, mss_2], # g/mol
    box_dimensions=[L, L, L], # A
    cut_off=rc,
)
minimizer.run()

# Test function using pytest
def test_energy_and_force():
    Final_Epot = minimizer.Epot
    Final_MaxF = minimizer.MaxF
    assert Final_Epot < 0, f"Test failed: Final energy too large: {Final_Epot}"
    assert Final_MaxF < 10, f"Test failed: Final max force too large: {Final_MaxF}"
    print("Test passed")

# If the script is run directly, execute the tests
if __name__ == "__main__":
    import pytest
    # Run pytest programmatically
    pytest.main(["-s", __file__])

For such a low particle density, we can reasonably expect that the potential energy will always be negative after 100 steps, and that the maximum force will be smaller than 10 (unitless).