# How to compute forces in multi-particle MD

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?

• Your interaction should be repulsive if $r<r_0$ and attractive if $r>r_0$. From that you can infer the sign of your forces. – nicoguaro Dec 20 '17 at 17:54
• So you would check if $\sqrt{x^2+y^2+z^2}$ is greater than or less than $r_0$? – Ptheguy Dec 21 '17 at 2:59

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.

• Although this reflects mostly what BalazsToth has said, I think it is a more well-written and coherent answer. Additionally, my code looks very much like yours now, as I updated it since this post.Thank you! – Ptheguy Dec 28 '17 at 12:04

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.

• Thank you for the insightful remarks. You're very correct with regards to the cutoff approach. In fact, I am familiar with it but for this research project we like to actually consider all interactions. I'm simulating 50-100 atoms anyways. With respect to the AoS vs. SoA thing, are you suggesting that I use an array to store the data? So for example, double atom? The 9 indices will be three for position, three for velocity, and three for acceleration (accounting for each component). – Ptheguy Dec 21 '17 at 3:07
• An Array of Structures vs Structure of Arrays is never a hard fact. It really depends on the use of it. If you have complicated modelling which you must do on a (big) part of the structure for each structure, the cache efficiency of the Array of Structure may be more important. In other words: It depends on the use case. – André Dec 21 '17 at 11:47
• @Ptheguy Neighbour searching is then senseless, but in case of 50-100 atoms it does not seem to be an issue. The same stands for the AoS vs SoA problem. In case of such a low number of particles, i would not make serious efforts in optimisation or even parallelisation. – BalazsToth Dec 21 '17 at 11:55
• @Appleguy ok thank you. May I ask if my force calculation is correct. Particularly, is it correct to consider the three coordinates and add them all up to obtain the net effect? I think it is as it makes sense :) – Ptheguy Dec 21 '17 at 13:17
• @Ptheguy I think that is wrong in your code. You can calculate the forces separately in each directions, but you have to use the scalar distance d=sqrt(dx^2+dy^2+dz^2) instead of the distances in each direction. So you should calculate the force once for each pair or particles (f=force(d)) and then use it for the force-vector components: fx = dx/df, fy = dy/df, fz = dz/d*f. I'll edit it in my answer. – BalazsToth Dec 21 '17 at 13:36