How to compute forces in multi-particle MD

Question

Suppose we have a system of $N$ particles that interact via the Lennard-Jones potential $$V(r)=V_0\left[\left(\frac{r_0}{r}\right)^{12}-2\ \left(\frac{r_0}{r}\right)^{6}\right].$$ No other forces exist. The objective is to perform molecular dynamics simulations with this system.

I'm stuck on an important part. I do not know exactly how to calculate the forces on each atom. I know that to calculate the force on one atom, I must consider the interaction of it with all other atoms. However, I'm stuck as to how I put together these forces. Particularly, I don't know how to get the direction (sign) of the forces right.

Currently, my approach is to consider the $x$, $y$, and $z$ components of a given atom separately in calculating the forces/accelerations. I visit the interaction of one atom with all other atoms and compute the force in $x$, $y$, $z$. For each component, I add them all up, assuming this approach leaves me with a net acceleration in each direction at the end. I think I can use the final values of atm.ax, atm.ay, and atm.az to update the position and velocity of the atom.

void computeForce(Atom atm){
   for(int i=0;i<N_particles;i++){
   if(i==atm.index)
      continue;
   else{
      dx = PBC1D(atm.x - allAtm[i].x);
      dy = PBC1D(atm.y - allAtm[i].y);
      dz = PBC1D(atm.z - allAtm[i].z);

      atm.ax += force(dx) / m;
      atm.ay += force(dy) / m;
      atm.az += force(dz) / m;
   }
}

(Atom is a user-defined class; PBC1D allows for periodic boundary conditions---not important to my question here.)

Although my current approach makes sense to me, I have a feeling that it is incomplete and perhaps very inefficient. Can anyone help me understand if what I'm doing is correct or not? And, perhaps teach me better ways?

Answer 1

I will only be developing the computation of forces on a given atom (part 1 of BalazsToth answer), for this point is actually a bit tricky and I think it deserves to be a bit more explained. You can find useful remarks in any book of molecular simulation my personal reference being "Computational Simulation of Liquids" from Allen and Tildesley and the "Molecular Dynamics" Lecture Notes on Physics from W. Hoover.

The force of a pair of particles interacting is directed along the vector defined by the difference of the coordinates of the particles. Let's call this vector $\mathbf{r}_{ij}$ and the unitary vector with the same direction: $$ \mathbf{r}_{ij} = \mathbf{r}_i - \mathbf{r}_j\ ;\ \mathbf{n}_{ij}=\frac{\mathbf{r}_{ij}}{|\mathbf{r}_{ij}|}=\frac{\mathbf{r}_{ij}}{r}$$

With $|\mathbf{r}_{ij}|=\sqrt{r_{ij,x}^2+r_{ij,y}^2+r_{ij,z}^2}=r$

If you compute the force $f$ with regard to $r$ : $$f_{ij} = -\frac{dV}{dr}$$ You will note that :

1. This value is positive when $r<r_0$ and negative when $r>r_0$.
2. The force is pushing the atoms away from one another in the first case and getting them closer in the second.

Therefore the force exerced on atom $i$ is opposed in direction to vector $\mathbf{n}_{ij}$ as it is defined here. Some authors define it with opposite value, or for example, BalazsToth used $\mathbf{n}_{ji}$.

You can then add the force to the total force on the atom $i$ : $$ \mathbf{F}_{i} = \sum_{j=1\ ;\ j \ne i}^N -f_{ij}\mathbf{n}_{ij}$$

One remark though concerning your algorithm. Given that $f_{ij}=-f_{ji}$ due to Newton's third law, you can make your loop outside of your function computing the forces. You can actually make your piece of code more like :

void ComputeForces{
for(int i=0;i<N_particles;i++){
  for(int j=i+1;j<N_particles,j++){
    Computeforce(atm[i], atm[j]);
  }
}

In the function ComputeForce, taking two arguments allow you to shorten the loop as you go on and on as you can update the forces on atom $j$ as well.

Answer 2

I have the following remarks:

1. Force calculation

The equation of motion of the $i^\text{th}$ atom with mass $m_i$ reads as

$\textbf{a}_i=\frac{d\textbf{v}_i}{dt}=\frac{1}{m_i}\cdot f(\vert\textbf{r}_j-\textbf{r}_i\vert)\textbf{n}_{ji}$,

where the interparticle force $f(\vert\textbf{r}_j-\textbf{r}_i\vert)$ depends only on the distance between particles $i$ and $j$. In your implementation you calculate the force using the distances in each direction, which is incorrect. Instead, you need to compute the force for each pair of particles only once (using the scalar distance between them) and multiply it with the normalised direction vector $\textbf{n}_{ji}=\frac{\textbf{r}_j-\textbf{r}_i}{\vert\textbf{r}_j-\textbf{r}_i\vert}$ to obtain the components of the force-vector.

2. Neighours

It's not so efficient in case of large number of atoms because of the $\mathcal{O}(n^2)$ computational complexity.

Since the interaction between the atoms becomes smaller and smaller with increasing distance, it is beneficial to apply a certain threshold or cutoff distance to consider only the neighbouring ones of each atom. Based on the threshold you can divide your domain into threshold-sized cells and perform a neighbour search. It is a technique to find potential neighbours without the actual calculation of any of the distances between any pair of atoms. Based on the list of potential neighbours to be considered you can calculate the distances and check if they are really inside the influence radius or not. There are several methods to find potential neighbours, but in general $\mathcal{O}(n)$ complexity can be achieved. Read more about the neighbour search algorithms here.

Of course, ignoring the particles outside of the influence radius may cause errors in your simulation if your interaction law has infinite radius. Therefore, application of a threshold is a trade-off between accuracy and requirement of computational performance.

3. Container of particles

Your current implementation stores the particles as an AoS (array of structures). It is a convenient approach from programming point of view, but less efficient than the SoA (structure of arrays), especially for parallel computations. In the latter you would create a large structure with arrays inside storing all the particle properties. It ends up in a well-aligned data in the memory that is preferable for this type of calculations.