# Kinetic energy goes gradually to zero in my velocity Verlet MD simulation

I am trying to figure out what might have caused such a problem: I am simulating a binary Lennard-Jones mixture (two types of atoms interacting via LJ potential) in NVE ensemble with velocity Verlet algorithm. The problem that I have is that, my kinetic energy goes gradually (no sharp change at any single step) to zero, and the potential energy drops gradually (no sharp change at any single step) to a negative value, which is obviously not correct. I printed out the positions and velocities of the atoms, and I find that the velocities go to vanishing values and the positions almost do not change if running very long. The boundary condition is good, so the volume is still fixed. It basically means that in my simulation, the atoms eventually come to a stop.

It is easy to verify that my velocity Verlet algorithm is correct, but what do you think might be the cause of this kind of weird results? It looks like that the energy is dissipating in my simulation. The code for calculating the potential energy and the force is here, with a cutoff at distance rc, the parameters being the Kob Andersen parameters. Nsmall is the number of smaller atoms, while N is the total number of atoms. For each array (r, f), the first (N-Nsmall) elements are the values for the larger atoms, with the rest Nsmall elements for the smaller atoms.

double total_e ( double * rx, double * ry, double * rz,
double * fx, double * fy, double * fz,
int N, int Nsmall, double L) {
int i,j;
double dx, dy, dz, r2, rc2, r6i, rc6i, epsilon, sigma;
double e = 0.0, hL=L/2.0,f;

/* Zero the forces */
for (i=0;i<N;i++) {
fx[i]=fy[i]=fz[i]=0.0;
}

for (i=0;i<(N-1);i++) {
for (j=i+1;j<N;j++) {
/* Kob Andersen parameters here */
if (i < N - Nsmall && j < N - Nsmall){
epsilon = 1.0;
sigma = 1.0;
} else if (i >= N - Nsmall && j >= N - Nsmall){
epsilon = 0.5;
sigma = 0.88;
} else {
epsilon = 1.5;
sigma = 0.8;
}

dx  = (rx[i]-rx[j]);
dy  = (ry[i]-ry[j]);
dz  = (rz[i]-rz[j]);
/* Periodic boundary conditions: Apply the minimum image
convention; note that this is *not* used to truncate the
potential as long as there an explicit cutoff. */
if (dx>hL)       dx-=L;
else if (dx<-hL) dx+=L;
if (dy>hL)       dy-=L;
else if (dy<-hL) dy+=L;
if (dz>hL)       dz-=L;
else if (dz<-hL) dz+=L;

r2 = dx*dx + dy*dy + dz*dz;
rc2 = 6.25*sigma*sigma; /* Kob Andersen parameter for rcut*/
if (r2<rc2) {
r6i   = 1.0/(r2*r2*r2);
rc6i = 1.0/(rc2*rc2*rc2); // cutoff correction to the potential and to the force
e    += 4*epsilon*(pow(sigma,12)*r6i*r6i - pow(sigma,6)*r6i) - 4*epsilon*(pow(sigma,12)*rc6i*rc6i - pow(sigma,6)*rc6i) + 48*epsilon*(pow(sigma,12)*rc6i*rc6i/(sqrt(rc2))-0.5*pow(sigma,6)*rc6i/(sqrt(rc2)))*(sqrt(r2) - sqrt(rc2));
f     = 48*epsilon*(pow(sigma,12)*r6i*r6i-0.5*pow(sigma,6)*r6i) - 48*epsilon*(pow(sigma,12)*rc6i*rc6i-0.5*pow(sigma,6)*rc6i);
/* Newton's third law*/
fx[i] += dx*f/r2;
fx[j] -= dx*f/r2;
fy[i] += dy*f/r2;
fy[j] -= dy*f/r2;
fz[i] += dz*f/r2;
fz[j] -= dz*f/r2;
}
}
}
return e;


}

By the way, another question is, should the total energy be strictly constant in a velocity Verlet numerical simulation?

• Would Computational Science be a better home for this question? – Qmechanic Oct 8 '16 at 20:49
• The debugging of code is certainly not on-topic here (do consider Computational Science for that), through the general topic is. – dmckee Oct 8 '16 at 21:23
• Aha, I don't even know there exists such a site. Thanks! – Damien Oct 8 '16 at 21:40

Energy will not necessarily be conserved exactly in any dynamics simulation, because the numerical algorithm is "sampling" the value of some quantity (in your code, apparently the forces fx, fy, fz) during each time step and assuming it remains constant over the whole of the finite length step. In the real world situation, the force is changing continuously during the time step.
"Conservation of energy" may not be the whole story in any case. For example, there are well known numerical integration methods which conserve energy "exactly" for simple problems like a linear harmonic oscillator, but the constant value of the energy is a function of the time step size, and the incorrect amount of energy in the numerical system results in an error in the oscillation frequency compared with the original differential equation - i.e. the frequency of the numerical solution is not exactly $\omega = \sqrt{k/m}$ but rather $\omega = \sqrt{k/m} + f(h)$ where $h$ is the numerical time step.