# 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? –  Juan M. Bello Rivas Feb 22 '13 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? –  Deathbreath Feb 22 '13 at 19:52
@Juan and Deathbreath: I have updated my question based on your comments. –  Sankaran Feb 23 '13 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. –  Pedro Feb 23 '13 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. –  Pedro Feb 23 '13 at 20:22

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.