Take the 2-minute tour ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

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;
}
share|improve this question
    
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
show 7 more comments

1 Answer

up vote 3 down vote accepted

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.

share|improve this answer
    
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? –  Sankaran Feb 25 '13 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... –  Pedro Feb 25 '13 at 10:49
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.