如何使velocity verlet算法更稳定？

计算科学 计算物理学 分子动力学
2021-12-24 03:17:06

这个问题的答案意味着减少时间步长会使它更稳定。但是我尝试减少时间步，但系统仍然不稳定（总能量增加到一个非常大的值）。我的理解是，如果两个球非常接近并且在某个步骤结束时，在下一步它们将具有非常大的速度，从而使系统不稳定。如何选择时间步以使系统保持稳定？如果有人可以查看我提供的代码以确定算法的实现是否存在错误，我也将不胜感激？

program mol_dyn_ref
  implicit none
  double precision,allocatable,dimension(:) :: posx,posy,posz,velx,vely,velz,ax,ay,az,tempx,tempy,tempz
  integer,allocatable,dimension(:) :: seed
  double precision :: x,y,z,m,rad,eps,kin,pot,k,p
  double precision :: xinit,yinit,zinit,xlen,ylen,zlen,mindist
  integer :: i,j,st,seedsize,n,iter,t
  double precision :: tot_t,dt,dist,r,Fx,Fy,Fz

  interface
     subroutine acc(n,posx,posy,posz,m,eps,tempx,tempy,tempz)
       implicit none
       integer :: i,j,n
       double precision,dimension(n) :: posx,posy,posz,ax,ay,az,tempx,tempy,tempz
       double precision :: r,dist,Fx,Fy,Fz,eps,m
     end subroutine acc

     subroutine pos_upd(n,posx,posy,posz,velx,vely,velz,ax,ay,az,dt,xlen,ylen,zlen,rad)
       implicit none
       double precision,dimension(n) :: posx,posy,posz,velx,vely,velz,ax,ay,az
       integer :: n
       double precision :: tot_t,dt,xlen,ylen,zlen,rad
     end subroutine pos_upd

     subroutine vel_upd(n,velx,vely,velz,ax,ay,az,dt,tempx,tempy,tempz)
       implicit none
       integer n
       double precision,dimension(n) :: velx,vely,velz,ax,ay,az,tempx,tempy,tempz
       double precision :: dt
     end subroutine vel_upd
  end interface



  n=200          !Number of particles
  eps=1.d0
  m=1.d0
  xlen=20.d0
  ylen=20.d0
  zlen=20.d0
  xinit=0.d0
  yinit=0.d0
  zinit=0.d0
  mindist=1.5d0
  rad=mindist/2.d0
  tot_t=10.d0
  dt=0.001d0

  iter=int(tot_t/dt)

  allocate(posx(n),posy(n),posz(n),velx(n),vely(n),velz(n),ax(n),ay(n),az(n),tempx(n),tempy(n),tempz(n))

  open(100,file="pos.dat",status="replace")
  open(200,file="vel.dat",status="replace")
  open(300,file="acc.dat",status="replace")
  open(400,file="energy.dat",status="replace")

  call random_seed(size=seedsize)
  allocate(seed(seedsize))
  do i=1,seedsize
     call system_clock(st)
     seed(i)=st
  enddo
  call random_seed(put=seed)

  !Assigning initial position to first particle
10 call random_number(x)
  posx(1)=xinit+x*xlen
  call random_number(y)
  posy(1)=yinit+y*ylen
  call random_number(z)
  posz(1)=zinit+z*zlen

  if(posx(1)<rad .OR. posx(1)>xlen-rad) goto 10
  if(posy(1)<rad .OR. posy(1)>ylen-rad) goto 10
  if(posz(1)<rad .OR. posz(1)>zlen-rad) goto 10

  !Assigning initial position
  do i=2,n
20   call random_number(x)
     posx(i)=xinit+x*xlen
     call random_number(y)
     posy(i)=yinit+y*ylen
     call random_number(z)
     posz(i)=zinit+z*zlen
     if(posx(i)<rad .OR. posx(i)>xlen-rad) goto 20
     if(posy(i)<rad .OR. posy(i)>ylen-rad) goto 20
     if(posz(i)<rad .OR. posz(i)>zlen-rad) goto 20
     do j=1,i-1
        if (dist(posx(i),posy(i),posz(i),posx(j),posy(j),posz(j))<mindist) goto 20
     enddo
  enddo
  print*, "Position initialisation finished"

  !Assigning initial velocities
  do i=1,n
     call random_number(x)
     velx(i)=02.1d0*(2.d0*x-1.d0)
     call random_number(y)
     vely(i)=02.1d0*(2.d0*y-1.d0)
     call random_number(z)
     velz(i)=02.1d0*(2.d0*z-1.d0)
   enddo
  print*, "Velocity initialisation finished"

  !Assigning initial acceleration
 do i=1,n
     ax(i)=0
     ay(i)=0
     az(i)=0
     do j=1,n
        if (j==i) cycle
        r=dist(posx(i),posy(i),posz(i),posx(j),posy(j),posz(j))
        ax(i)=ax(i)+ Fx(r,posx(i)-posx(j),posy(i)-posy(j),posz(i)-posz(j),eps)/m
        ay(i)=ay(i)+ Fy(r,posx(i)-posx(j),posy(i)-posy(j),posz(i)-posz(j),eps)/m
        az(i)=az(i)+ Fz(r,posx(i)-posx(j),posy(i)-posy(j),posz(i)-posz(j),eps)/m
enddo
  enddo
  print*, "Acceleration initialisation finished."

!Molecular Dynamics Simulation
  do t=1,iter
     print*, t
     do i=1,n
     write(100,*) posx(i),posy(i),posz(i)
     write(200,*) velx(i),vely(i),velz(i)
     write(300,*) ax(i),ay(i),az(i)
     enddo
     k=kin(velx,vely,velz,m,n)
     p=pot(posx,posy,posz,eps,n)
     write(400,*) t,k,p,k+p

     call pos_upd(n,posx,posy,posz,velx,vely,velz,ax,ay,az,dt,xlen,ylen,zlen,rad)
     call acc(n,posx,posy,posz,m,eps,tempx,tempy,tempz)
     call vel_upd(n,velx,vely,velz,ax,ay,az,dt,tempx,tempy,tempz)
  enddo
  call system("gnuplot --persist plot.gp")
endprogram mol_dyn_ref


!Updating Acceleration
subroutine acc(n,posx,posy,posz,m,eps,tempx,tempy,tempz)
  implicit none

  integer :: i,j,n
  double precision,dimension(n) :: posx,posy,posz,ax,ay,az,tempx,tempy,tempz
  double precision :: r,dist,Fx,Fy,Fz,eps,m
  do i=1,n
     tempx(i)=ax(i)
     tempy(i)=ay(i)
     tempz(i)=az(i)
     ax(i)=0.d0
     ay(i)=0.d0
     az(i)=0.d0
     do j=1,n
        if (j==i) cycle
        r=dist(posx(i),posy(i),posz(i),posx(j),posy(j),posz(j))
        ax(i)=ax(i)+ Fx(r,posx(i)-posx(j),posy(i)-posy(j),posz(i)-posz(j),eps)/m
        ay(i)=ay(i)+ Fy(r,posx(i)-posx(j),posy(i)-posy(j),posz(i)-posz(j),eps)/m
        az(i)=az(i)+ Fz(r,posx(i)-posx(j),posy(i)-posy(j),posz(i)-posz(j),eps)/m
     enddo
  enddo
end subroutine acc

!Updating Position
subroutine pos_upd(n,posx,posy,posz,velx,vely,velz,ax,ay,az,dt,xlen,ylen,zlen,rad)
  implicit none

  integer n,i
  double precision,dimension(n) :: posx,posy,posz,velx,vely,velz,ax,ay,az
  double precision :: dt,xlen,ylen,zlen,rad
     do i=1,n
        posx(i)=posx(i) + velx(i)*dt + (ax(i)*(dt**2))/2
        posy(i)=posy(i) + vely(i)*dt + (ay(i)*(dt**2))/2
        posz(i)=posz(i) + velz(i)*dt + (az(i)*(dt**2))/2

        if(posx(i)<rad .OR. posx(i)>xlen-rad) velx(i)=-velx(i)
        if(posy(i)<rad .OR. posy(i)>ylen-rad) vely(i)=-vely(i)
        if(posz(i)<rad .OR. posz(i)>zlen-rad) velz(i)=-velz(i)
     enddo
end subroutine pos_upd

!Updating Velocity
subroutine vel_upd(n,velx,vely,velz,ax,ay,az,dt,tempx,tempy,tempz)

  implicit none
  integer n,i
  double precision,dimension(n) :: velx,vely,velz,ax,ay,az,tempx,tempy,tempz
  double precision :: dt
  do i=1,n
    velx(i)=velx(i) + 0.5d0*(ax(i)+tempx(i))*dt
    vely(i)=vely(i) + 0.5d0*(ay(i)+tempy(i))*dt
    velz(i)=velz(i) + 0.5d0*(az(i)+tempz(i))*dt
  enddo
end subroutine vel_upd

function pot(posx,posy,posz,eps,n)
  implicit none
  double precision pot,r,dist,eps
  integer i,j,n
  double precision, dimension(n) :: posx,posy,posz
  pot=0.d0
  do i=1,n
  do j=1,n
  if(i==j) cycle
  r=dist(posx(i),posy(i),posz(i),posx(j),posy(j),posz(j))
  pot= pot + (4.d0*eps*((1.d0/r)**12 - (1.d0/r)**6))                 !r is relative distance. x,y,z are components of r.
  enddo
  enddo
end function pot

 function kin(velx,vely,velz,m,n)
 implicit none
  double precision :: kin,m
  integer :: i,n
  double precision, dimension(n) :: velx,vely,velz
  kin=0.d0
  do i=1,n
  kin = kin + ((velx(i)**2 + vely(i)**2 + velz(i)**2)/(2*m))
  enddo
end function kin

function dist(x1,y1,z1,x2,y2,z2)
  implicit none
  double precision :: dist,x1,y1,z1,x2,y2,z2
  dist = sqrt((x1-x2)**2 + (y1-y2)**2 + (z1-z2)**2)
end function dist

function Fx(r,x,y,z,eps)
  implicit none
  double precision :: Fx,r,x,y,z,eps
  Fx = 4.d0*eps*((12.d0/r**14) - (6.d0/r**8))*x
end function Fx

function Fy(r,x,y,z,eps)
  implicit none
  double precision :: Fy,r,x,y,z,eps
  Fy = 4.d0*eps*((12.d0/r**14) - (6.d0/r**8))*y
end function Fy

function Fz(r,x,y,z,eps)
  implicit none
  double precision :: Fz,r,x,y,z,eps
  Fz = 4.d0*eps*((12.d0/r**14) - (6.d0/r**8))*z
end function Fz
1个回答

您可以尝试使用自适应时间步长，dt如果速度太大，则在每一步调整位置。然而，Verlet 方法是辛的，所以如果总能量不是大致守恒的，那么就有问题了。如果两个粒子靠得太近，没有合理的时间步长会使方法稳定，用双精度浮点数精确计算粒子之间的力是不可能的。

虽然看起来您正在随机初始化粒子位置以使它们彼此不会太近，但您似乎没有对粒子速度进行任何防护以防止粒子在第一次更新时过于靠近。您使用的 Lenard-Jones 势在近距离时非常僵硬，因此如果您在初始化系统时不小心，您可能会立即让一个粒子飞到无穷远。

相反，尝试将系统初始化为“冷”状态：将所有速度设置为零并将粒子置于六边形密排排列中，距离等于您选择的势能的平衡距离。你会看到一些关于真正平衡配置的振荡。这将是一个更容易解决的系统，以便您可以验证算法的正确性和稳定性。从那里，您可以慢慢增加初始排列的速度，直到系统达到液态或气态。

您没有询问编码风格，但同样，我建议将您的代码组织成模块并声明intent每个过程参数。Fortran Best Practices提供了有关此主题的更多信息。

其它你可能感兴趣的问题
上一篇凸二次问题求解器给出不同的答案？ 下一篇自适应网格细化基本概念问题