# Difficulty with possibly the simplest MD simulation

I am going ahead with asking what might appear as a debugging question. But,I am looking for something conceptual that I might have missed, and this question is somewhat related to my previous question: NVE MD simulation of inert gas: Problem maintaining equilibrium, although a different issue now. If the moderators decide to close this, I understand.

I am doing a constant energy (NVE) simulation of 49 Argon particles in a 100x100 2D box with Lennard Jones potential at temperature 240 K I am using the following:

1) Unit system: Natural unit system such that mass $m=1$, LJ parameters, $\sigma=1$, $\epsilon=1$, Boltzmann constant $k=1$ such that temperatureof 120 K is 1.

2) Initial Conditions: Square lattice position, Maxwell-Boltzmann velocities

3) Integration: Velocity-verlet with periodic boundary conditions. Time step of 0.001 (2fs in actual time)

4) I am checking for occurrence of a non-zero mean velocity (flying ice cubes): It does not occur.

I am unable to keep energy conservation within 1% for more than 2000 time steps (total actual time of 4ps). Most studies, lecture notes, online tutorials,etc., seem to be using these same values for time step, parameters, etc. While these never say how long they are able to conserve energy (and within what accuracy) I find research papers do 10ns on an average. I am failing 3 orders below. Velocity verlet seems like a fairly-accepted algorithm at least for something so simple!

It doesn't seem like a basic coding (i.e.,non-conceptual) mistake since the code does well until it gets to around 4ps. Some runs I can go upto 6 or 7ps, just lucky initial conditions I guess. Any helpful suggestions on where I could be going wrong? Thanks.

Update: 1) The error doesn't seem to scale with time-step $\Delta t$. Most papers seem to use $\Delta t$ =0.001, I have tried the range from 1e-2 to 1e-5.

2) My integrator code looks like this: { while(t

  //Energy Calculation
Ven(t)=ven;
KE(U,V,ken,N);
Ken(t)=ken;
if(t==0)
{
Hen1=Ken(0)+Ven(0);
cout<<"Total energy at t=0  is: "<<Hen1<<endl;
}

//Error Checking
ErrMeanVel=(U.sum())/N;
ErrEConserve=(abs(Ken(t)+Ven(t)-Hen1)*100)/Hen1;
if(ErrMeanVel>1e-6)
{
cout<<"Flying ice cube at t="<<t*Deltat<<"ps. The mean velocity is"<<ErrMeanVel<<endl;
break;
}
if(ErrEConserve>ETol)
{
cout<<"Energy conservation violated at t="<<t*Deltat<<"ps. The energy at this time is "<<(Ken(t)+Ven(t))<<" and the percentage error is "<<((abs(Ken(t)+Ven(t)-Hen1)*100)/Hen1)<<endl;
break;
}

//Velocity Verlet
for (i=0;i<N;i++)
{
X(0,i)=X(0,i)+(U(0,i)*Deltat)+0.5*Ax(i)*Deltat*(Deltat/mass);
Y(0,i)=Y(0,i)+(V(0,i)*Deltat)+0.5*Ay(i)*Deltat*(Deltat/mass);
U(0,i)=U(0,i)+(0.5*Ax(i)*Deltat/mass);
V(0,i)=V(0,i)+(0.5*Ay(i)*Deltat/mass);
}
LJ(X,Y,N,ven,Fx,Fy,ctr);
Ax=Fx.rowwise().sum();
Ay=Fy.rowwise().sum();
for (i=0;i<N;i++)
{
U(0,i)=U(0,i)+(0.5*Ax(i)*Deltat/mass);
V(0,i)=V(0,i)+(0.5*Ay(i)*Deltat/mass);
}
//Boundary Conditions
PeriodicBoundary(X,Y,L,N);
t++;
}


\Force calculator is: (no cutoff)

 int i=0, j=0;
double dx=0, dy=0,r2=0,F6=0,F12=0,F=0,Vcut=0;
fx=MatrixXd::Zero(n,n);
fy=MatrixXd::Zero(n,n);
ven=0;
MatrixXd ft(n,n);
for (j=0;j<n;j++)
{
for(i=j+1;i<n;i++)
{
dx=(a(0,j)-a(0,i));
dy=(b(0,j)-b(0,i));
r2=dx*dx+dy*dy;
if(r2<sigma)
{
ctr++;
}
F6=pow((sigma*sigma/r2),3);
F12=pow((sigma*sigma/r2),6);
F=(24*epsilon/r2)*(F6-2*F12);
ven=ven+4*epsilon*(F12-F6);
fx(i,j)=F*dx;
fy(i,j)=F*dy;
}
}
ft=fx.transpose();
fx=fx-ft;
ft=fy.transpose();
fy=fy-ft;
}

• Can you post the code for your integrator? Commented Feb 22, 2013 at 3:31
• Have you done an error analysis to estimate what energy drift you should expect given machine precision? What happens when you reduce $\Delta t$ to 0.0001? Commented Feb 22, 2013 at 19:52
• @Juan and Deathbreath: I have updated my question based on your comments. Commented Feb 23, 2013 at 19:01
• Glad to see you're progressing with your simulations! Could you try plotting both the potential and kinetic energy for each time step? Both curves should appear smooth. If, however, they have jumps, they will at least tell you where and when your code is going wrong. Commented Feb 23, 2013 at 19:28
• Ok, forget my last comment, the bug is in F=(24*epsilon/r2)*(F6-2*F12). You should be dividing by $r$, not $r^2$. The mismatch between the energy function and its derivative is most probably what's causing you problems with the total energy. Commented Feb 23, 2013 at 20:22

This is a re-formulation of my comments above as an answer.

You need to enforce the boundary conditions when computing dx and dy to account for particles that interact over the periodic boundary. Imagine two particles at opposite ends of the domain: In your current setup they won't interact, but as soon as one of them moves slightly and crosses the boundary, they will all of a sudden see each other.

Most simulations just use what's referred to as the minimum image convention.

• Interestingly I can conserve energy below 0.01% for say 49 molecules or 100 molecules, but can't for fewer molecules. I can conserve better energy for higher temperatures, i.e., error is less significant in comparison with kinetic energies. That seems about right. What do you think? Commented Feb 25, 2013 at 3:33
• @Sankaran: I'd recommend you test your code with just two atoms interacting and/or moving across the periodic boundaries, keeping an eye on both the potential and kinetic energies. If you see a jump in either anywhere, you can then trace it back to what actually happened, e.g. a particle moving across a boundary, within or outside of a cutoff distance, etc... Commented Feb 25, 2013 at 10:49

Two thoughts

A) A certain degree of energy dissipation in MD simulations is considered normal and can be counterbalanced with so called thermostat algorithms e.g. this. In your particular case I can not say wheter the ammount of dissipation you have is to be expected in your setup, or result of some bug.

B) The Lennard Jones potential is very steep when particles are near each other and diminishes fast when they are apart. That means that upon collision you would need to have very small timesteps to get reasonable results. Have you tried running your simulation with smaller timesteps just to eliminate the possiblilty that it's simply the temporal resolution?