# Error in MD simulation with Lennard-Jones potential in C, resulting in a lack of conservation of energy

I'm working on a MD simulation using the Lennard-Jones potential and the Velocity Verlet algorithm. The code compiles and runs, but the energy for each time step is not conserved; it decreases constantly when run. The magnitude of the potential energy is much larger than that of the kinetic energy, but I haven't been able to detect where the error is.

#include <stdio.H>
#include <math.H>
#include <stdlib.H>
struct particle
{
double x[3],v[3],a[3];
};

double sigma, epsilon, m, dt;
#define N 1000
struct particle particles[N];

double rfinder(int id_a, int id_b)
{
int i;
double rtot,dist,dist2;
dist2 = 0;
for(i=0;i<3;i++)
{
dist = particles[id_a].x[i] - particles[id_b].x[i];
dist2 += dist * dist;
}
rtot = sqrt(dist2);
return rtot;
}

double calcacceleration()
{
int b,c,i;
double dist,r,r7,r13,sigma6,sigma12, r12, sor6, sor12, pe;
pe = 0;
for(b=0;b<N;b++)
{
for(i=0;i<3;i++)
{
particles[b].a[i] = 0;
}
}

for(b=0;b<N;b++)
{
for(c=0;c<N;c++)
{
if(b!=c)
{
double dist2;
dist2 = 0;
for(i=0;i<3;i++)
{
dist = particles[b].x[i] - particles[c].x[i];
dist2 += dist * dist;
}
if (dist2<8*8)
{
r = sqrt(dist2);
//printf("%d %d %f\n",b,c,r);
double rinv = 1.0/r;
double r2 = rinv*rinv;
double r4 = r2*r2;
double r6 = r4 * r2;
r7 = r4*r2*rinv;
r13 = r7*r4*r2;
r12 = r6 * r6;
sigma6 = sigma*sigma*sigma*sigma*sigma*sigma;
sigma12 = sigma6*sigma6;
sor6 = sigma6*r6;
sor12 = sigma12*r12;
for(i=0;i<3;i++)
{
dist = particles[b].x[i] - particles[c].x[i];
particles[b].a[i] += -(0.5 * (4 * epsilon *((-12 * sigma12 * r13) +(6 * sigma6 * r7))) * ((dist/r))/m);
particles[c].a[i] +=  (0.5 * (4 * epsilon *((-12 * sigma12 * r13) +(6 * sigma6 * r7))) * ((dist/r))/m);
pe += 0.5*(4*epsilon*(sor12 - sor6));
}
}
}
}
}
return pe;
}

void position_mover()
{
double x, x0, v;
int p, i;
for(p=0;p<N;p++)
{
for(i=0;i<3;i++)
{
particles[p].x[i] = particles[p].x[i] + (particles[p].v[i] * dt);
}
}
return;
}

double velocity_mover()
{
double ke;
int p, i;
ke = 0;
for(p=0;p<N;p++)
{
for(i=0;i<3;i++)
{
particles[p].v[i] = particles[p].v[i] + (0.5 * particles[p].a[i] * dt);
ke += 0.5 * m * particles[p].v[i] * particles[p].v[i];
}
}
return ke;
}

void output_to_file()
{
int p;
FILE * fp;
fp = fopen("data.xyz", "a");
fprintf(fp,"%d \n\n",N);
for(p=0;p<N;p++)
{
fprintf(fp,"H %f %f %f\n",particles[p].x[0],particles[p].x[1],particles[p].x[2]);
}
fclose(fp);
}

int main()
{
int p, i, j, k;
double ft,t,pe, ke,te;
sigma = 1;
epsilon = 1;
FILE * fp;
fp = fopen("data.xyz", "w");
fclose(fp);
printf("What timestep would you like?\n");
scanf("%lf", &dt);
printf("In seconds, what is the final time?\n");
scanf("%lf", &ft);
printf("What mass would you like?\n");
scanf("%lf", &m);
p = 0;
//sets initial positions of atoms to be evenly spaced in a 10 * 10 * 10 cube
for(i=0;i<10;i++)
{
for(j=0;j<10;j++)
{
for(k=0;k<10;k++)
{
particles[p].x[0] = i*2;
particles[p].x[1] = j*2;
particles[p].x[2] = k*2;
p++;
}
}
}

/*sets initial velocities and accelerations of atoms to 0 to get rid of any
garbage in memory */
for(i=0;i<3;i++)
{
for(p=0;p<N;p++)
{
particles[p].v[i] = 0;
particles[p].a[i] = 0;
}
}

for(t=0;t<ft;t+=dt)
{
output_to_file();
calcacceleration();
velocity_mover();
position_mover();
pe = calcacceleration();
ke = velocity_mover();
te = ke + pe;
printf("Time: %f\n PE: %f\n KE: %f\n TE: %f\n",t,pe, ke, te);
}
printf("Done!\n");
return 0;
}

• Are you using a time-step that is small enough? – nicoguaro Jun 30 '16 at 21:58
• @nicoguaro I've been using .005. – L L Jun 30 '16 at 22:04
• What are the units? You are not showing the system, so I can't tell that 0.005 is small or big. How do you determine that value? What happens when you increase it or decrease it? – nicoguaro Jun 30 '16 at 22:13
• @nicoguaro I'm using reduced units. A timestep of 1 results in the PE going to 0 after 18 seconds while the kinetic energy blows up. A tilmestep of of .000005 results in the same constant decrease. – L L Jun 30 '16 at 22:16
• There are infinitely many reduce units systems. I was trying to get the information for people to help you. But if you don't want to provide it, it's your choice. You are the one asking the question... – nicoguaro Jun 30 '16 at 22:18