计算科学 - 可能是最简单的 MD 模拟的困难 - 吾爱随笔录

我将继续询问可能出现的调试问题。但是，我正在寻找一些我可能错过的概念，这个问题与我之前的问题有些相关：惰性气体的 NVE MD 模拟：保持平衡的问题，尽管现在是一个不同的问题。如果版主决定关闭它，我理解。

我正在使用以下方法对 100x100 2D 盒子中的 49 个氩粒子进行恒定能量 (NVE) 模拟，Lennard Jones 势在 240 KI 温度下：

1) 单位系统：自然单位系统，质量，LJ 参数，，，玻尔兹曼常数，使得 120 K 的温度为 1。 $m=1$ $\sigma=1$ $\epsilon=1$ $k=1$

2）初始条件：方格位置，麦克斯韦-玻尔兹曼速度

3) 积分：具有周期性边界条件的 Velocity-verlet。时间步长 0.001（实际时间为 2fs）

4）我正在检查非零平均速度（飞冰块）的发生：它没有发生。

我无法在超过 2000 个时间步长（实际总时间为 4ps）内将能量守恒保持在 1% 以内。大多数研究、讲义、在线教程等似乎都在使用这些相同的时间步长、参数等值。虽然这些从未说明它们能够保存多长时间（以及在什么精度范围内），但我发现研究论文确实如此平均10ns。我在以下 3 个订单中失败了。至少对于如此简单的事情来说，Velocity verlet 似乎是一个被广泛接受的算法！

这似乎不是一个基本的编码（即非概念性）错误，因为代码在达到大约 4ps 之前运行良好。有些运行我可以达到 6 或 7ps，我猜这只是幸运的初始条件。关于我可能出错的地方有什么有用的建议吗？谢谢。

更新： 1）错误似乎与时间步长无关。大多数论文似乎使用 =0.001，我尝试了从1e-2到1e-5的范围。 $\Delta t$ $\Delta t$

2) 我的积分器代码如下所示： { 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++;
}

\力计算器是：（无截止）

 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;
}