The Lennard-Jones Potential

The Model Fluid

The topic of the fundamental molecular dynamics of non-bonded interactions is studied using a soft-disk fluid simulation, which is an implementation of the Equation of Motion in which a spherical particle (atoms) interact with one another. Interactions between pairs of atoms are calculated for providing the two principal features of interatomic forces:

  • Repulsive forces: resistance to compression.

  • Attractive forces: mutual attraction of pairs of particles over a range of separations (dipole-dipole, dipole-induced dipole, van der Waals forces encompassing London dispersion forces).

    The genral equation that calculates these features is the Lennard-Jones potential (LJP).

    • This potential governs the strength of the interaction and repulsion for a pair of atoms $i$ and $j$ located at $r_i$ and $r_j$ positions. It calculates their accelerations and forces based on their distance of separation $r$.

    $ u(r_{ij}) = 4ϵ \bigg[ \bigg( \frac{\sigma}{r_{ij}}\bigg)^{12} - \bigg(\frac{\sigma}{r_{ij}}\bigg)^{6} \bigg]$

    The LJP model is composed of two components: the repulsive $(σ/r)¹²$ and the attractive term $(σ/r)⁶$, which respectively denote repulsive and attractive forces. The parameter $r$ is the distance between the two atoms (in $Å$ units);

    $σ$ is a length scale representing the distance at which the intermolecular potential between the two atoms is $= 0$ (in $Å$ units), and $ε$ governs the strength of the interaction (in $eV$ units). In essence, it’s a measure of how strongly two atoms attract each other.

    Finally, $u$ is the intermolecular potential between the two particles. The interaction repels at close range, then attracts, and is cut off at some limiting separation $r_c$: as the parameter $r$ increases towards $r_c$, the force drops to $0$.

    The LJP Curve

    For computational simplicity, Rapaport has modified the LJP in the example reported in Chapter Two, simplifying the interaction by ignoring the attractive tail represented by the van der Waals forces. Moreover, $ε$ and $σ$ are set to $1$.

    The Limiting separation is:

    $ r_{ij} < r_{c} = 2^{1/6} σ$

    If we construct the model fluid with this kind of potential, the simulation will look like a little more than a collection of $\textbf{colliding softballs}$

    Including Newtonian Mechanics

    Considering Newton's second law of motion, where $F$ is the force for each pair of particles at a position $r$. The Equation becomes:

    $ F_{i}(r_{ij}) = 48 ϵ \bigg[ \bigg( \frac{\sigma^{12}}{r_{ij}^{13}}\bigg) - 24 ϵ \bigg(\frac{\sigma^6}{r_{ij}^7}\bigg) \bigg]$

    Also, Newton’s third law implies that $F_{ji}$ is equal to -$F_{ij}$, so the force on particle $i$ from the pairwise interaction $u(r_{ij})$ has the opposite direction of the force on particle $j$. The calculation of forces will need to be numerically integrated, and it will allow us to derives coordinates, velocities, and accelerations of each atom within the simulation.

    Soft-Disk Fluid Algorithm for a 2-Dimension

    1. The Main coordinates the flux, firstly calling SetParams and SetupJob that are two functions for the program initialization. SetParams serves to set many global parameters, while SetupJob embodies InitCoords, InitVels, and InitAccels, which are functions for the initialization of the coordinates, the velocities, and the accelerations of all the atoms, respectively.
    1. Everything Main has to do will be calling SingleStep, which represents the function that handles the whole process. SingleStep will call LeapFrogStep, which performs the integration of Equation of Motion using a simple numerical technique: the leapfrog method. This method has excellent energy conservation properties, integrating the coordinates and velocities of the particles. The LeapFrogStep appears twice in the listing of SingleStep, with the argument 1 or 2 that determines which portion of the two-step leapfrog process is to be performed. Finally, SingleStep encompasses ComputeForces, which implements the LJP and the forces and updates atom accelerations, and two functions for the properties measurements: EvalProps and AccumProps

    Introducing one Python class for the molecule definition (class Mol) and one for the properties (class Prop).

    Replacing all the functions for the vector operations with linear algebra functions from NumPy. Also, written all the necessary functions for the randomness and for updating the coordinates in periodic boundaries

    The initialization functions for coordinates, velocities, and accelerations, are included in the following code

    The two functions SetParams and SetupJob are presented

    Introducing the attractive tail represented by the van der Waals forces

    Below code includes the leapfrog method for the Integration of the Equation of Motion

    Below code provides all the functions for the properties measurements (Temperature, Energy, and Pressure).

    Function for plotting the trajectories of the atoms (plotMolCoo), which makes all the plots step by step, and a second function that creates an mp4 video from all the coordinates plot tiles. The two functions have to work together. Also,including GraphOutput, a function for printing properties in a pandas dataframe shape.

    Here, the SingleStep function: the real gear of the whole algorithm.

    And the Main Loop

    A Soft-Disk fluid simulation based on the Lennard-Jones Potential represents a microscopic model of a fluid or gas. It is based on spherical particles simulating atoms that interact among them. The interactions occur between pairs of particles exploiting repulsive and attractive forces. As seen in the video, after 500 timesteps, the two colored atoms, which initially were close, tend to move in different directions.