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;
}
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:22dxanddy... – Pedro Feb 23 '13 at 20:46dxanddyto 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. Or how are you enforcing the boundary conditions? – Pedro Feb 23 '13 at 21:21