4
$\begingroup$

The answer to this question implies that reducing the time step would make it more stable. However I have tried reducing the time step but the system is still unstable(the total energy increases to a very large value). What I understand is that if two balls are very close and the end of some step, in the next step they will have very large velocity, thus making the system unstable. How can I choose the time step so that the system remains stable? I would also appreciate if someone could look into the code I provide to figure out if there is some error in implementation of the algorithm?

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
Improve this question
$\endgroup$

1 Answer 1

Reset to default
6
$\begingroup$

You can try using an adaptive timestep, where dt is adjusted at each step if the velocities get too large. However, the Verlet method is symplectic, so if the total energy is not roughly conserved then something is wrong. If two particles get too close, no reasonable timestep will make the method stable, and accurately computing the forces between particles with double-precision floats will be impossible.

While it looks like you're randomly initializing the particle positions so that they are not too close to each other, you don't appear to have any guards on the particle velocities to keep the particles from approaching too close at the first update. The Lenard-Jones potential that you're using is very stiff at close distances, so if you're not careful in how you initialize the system you can have one particle fly off to infinity right away.

Instead, try initializing the system in a "cold" state: set all the velocities to zero and put the particles in a hexagonally close-packed arrangement, with the distance equal to the equilibrium distance for the potential you've chosen. You'll see some amount of oscillation about the true equilibrium configuration. This will be a much easier system to solve so that you can verify the correctness and stability of your algorithm. From there, you can slowly increase the velocities of the initial arrangement until the system reaches the liquid or gaseous state.

You didn't ask about coding style, but all the same, I recommend organizing your code into modules and declaring the intent of every procedure argument. Fortran Best Practices has more good information on this subject.

Improve this answer
$\endgroup$
4
  • $\begingroup$ you don't appear to have any guards on the particle velocities to keep them from approaching too close at the next timestep Did you mean appraching too high? I thought potential should take care of velocity updation. Wouldn't changing velocities by hand change the simulation altogether? $\endgroup$
    – Yogesh Yadav
    May 23, 2015 at 18:16
  • $\begingroup$ In the main loop, you update the positions before updating the velocities or accelerations. While you have taken precautions to keep the particles far enough apart at the start, your first update for the particle positions could throw that out the window if the initial velocities are too big. You'd have to work out what "too big" is given the parameters of your simulation. The simplest thing is to set the velocities to zero and make sure that works first. In principle, the potential should take care of the issue, provided the system is initialized in a state that isn't too unstable. $\endgroup$
    – Daniel Shapero
    May 23, 2015 at 18:55
  • $\begingroup$ I tried setting the initial velocities to zero. It doesn't improve the situation. $\endgroup$
    – Yogesh Yadav
    May 23, 2015 at 19:25
  • 3
    $\begingroup$ Alright, you're still starting with 200 randomly-placed particles, which is fairly complex. Try the simplest thing possible -- 3 or 4 particles, very far apart. Then move them closer together so they actually interact. After that, add one particle at a time, with regular spacing between them. Then add randomization. At every stage, see what the maximum timestep is for your system to be stable. $\endgroup$
    – Daniel Shapero
    May 23, 2015 at 20:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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