Why Do Argon Atoms in a Toroidal Box Show Unusually Large Positional Values?

[llello]

So I'm working on a problem where you have 100 argon atoms in a box that obey the Lennard Jones potential. We're using the Verlet method to calculate the position of each particle at the next time step. The issue I'm having has to do with our boundary conditions. We assume that atoms are on a torus, that way we don't have to worry about binging into the walls.

So my problem arises when I try to calculate their positions at a later time. We're coding in Fortran. So the position of the ith particle is x(i), y(i) and its contained in an LxL box. I go through and calculate what the final positions are and then update those positions so that they're still inside the box. Just an fyi, I also have to keep around the old position for the Verlet method to work. The way I do this updating is as follows:

! adjust for torus
if (x(i) .gt. L) then
a = floor(x(i)/L)
x(i) = x(i) - a*L
xold(i) = xold(i) - a*L
endif
if (x(i) .lt. 0.0) then
a = ceiling(x(i)/L)
x(i) = x(i) + a*L
xold(i) = xold(i) + a*L
endif

if (y(i) .gt. L) then
a = floor(y(i)/L)
y(i) = y(i) - a*L
yold(i) = yold(i) - a*L
endif
if (y(i) .lt. 0.0) then
a = ceiling(y(i)/L)
y(i) = y(i) + a*L
yold(i) = yold(i) - a*L
endif


Now the length of my box is just 10.0. But for whatever reason, when I output x, y data for the final time, it gives me extremely large values (order of 10^8 and more). Am I oversimplifying the situation when the particle travels from one side of the box to the other or is this the correct way to update it? If it is correct, how is it I can be getting positions on the order of 10^8?

Many thanks to anyone who takes the time to read this.

-Louis-

 

[llello]

Shameless bump of thread.

 

[llello]

Final shameless bump before abandoning hope for a reply.

 

[Vanadium 50]

It can take more than 94 minutes for an answer.

 

[llello]

So I've convinced myself that my updating of positions is valid. However, I am still receiving final positions outside of the box. This is my entire code. Perhaps someone can find the error. Note, I intend to randomize the starting positions and velocities once I get this to actually work.!
! Program to calculate properties of an argon gas
! using the Verlet method in molecular dynamics.
!

program assn6no1
implicit none

double precision x(100), y(100), velx(100), vely(100), xold(100), yold(100), energy(10000), velcm(10000)
double precision storex(10000,100), storey(10000,100), time(10000), accelx(100), accely(100)
double precision xnew(100), ynew(100)
double precision L, dt, tempx, tempy, loopx, loopy, loopdistance, cmxtemp, cmytemp, cosine, sine
double precision lineardistance, distance, mass, answer, energytemp, xaccel, yaccel
double precision xdistance, ydistance, p, q, r, s
integer N, i, j, k, timesteps, a

! n = number of argon atoms, L = size of box side (10sigma, sigma=1)

N = 100
L = 10.0
dt = 0.005
mass = 1.0
timesteps = 1000

!
! initial velocities
!

do i=1, N
velx(i) = 0.0
vely(i) = 0.0
enddo

!
! Making the initial positions of 100 argon atoms and x/y at time -dt
!

do j=1, 10
do i=1, 10
x(i+(j-1)*10) = (i+0.0) - 1.0/2.0
y(j+(i-1)*10) = (i+0.0) - 1.0/2.0
enddo
enddo

!
! Slightly offsetting each position by 1/10*sigma and calculating x/y(0)
!

do i=1, N
if (mod(i,2) .eq. 1) then
x(i) = x(i) + 1.0/100.0
y(i) = y(i) + 1.0/100.0
else
x(i) = x(i) - 1.0/100.0
y(i) = y(i) - 1.0/100.0
endif
enddo

do i=1, N
xold(i) = x(i) - velx(i)*dt
yold(i) = y(i) - vely(i)*dt
enddo

!
! begin calculating motion (time loop)
!
do k=1, timesteps

!
! Calculate the acceleration, first initialize
!
do i=1, N
accelx(i) = 0.0
accely(i) = 0.0
enddo

! first loop is for the ith particle's acceleration
do i=1, N

! calculating a from every other particle
xaccel = 0.0
yaccel = 0.0
tempx = 0.0
tempy = 0.0
do j=1, N

! if statement prevents caclulating acceleration from the same particle
if (i .ne. j) then
cosine = minDist(x(i),x(j))/sqrt( minDist( x(i),x(j) )**2 + minDist( y(i),y(j) )**2 )
sine = minDist(y(i),y(j))/sqrt( minDist( x(i),x(j) )**2 + minDist( y(i),y(j) )**2 )
tempx = xaccel + 1.0/mass*force(x(i),x(j),y(i),y(j))*cosine
tempy = yaccel + 1.0/mass*force(x(i),x(j),y(i),y(j))*sine
xaccel = tempx
yaccel = tempy
endif
enddo
accelx(i) = xaccel
accely(i) = yaccel
enddo

!
! Calculate the new position of ith particle
!

do i=1, N

! computing new x/y
xnew(i) = 2.0*x(i) - xold(i) + accelx(i)*dt**2
ynew(i) = 2.0*y(i) - yold(i) + accely(i)*dt**2

! compute velcocity before adjusting for torus
velx(i) = (xnew(i) - xold(i))/(2*dt)
vely(i) = (ynew(i) - yold(i))/(2*dt)

enddo

! adjusting each particle for torus boundary

do i=1, N

if (xnew(i) .gt. L) then
a = floor(xnew(i)/L)
xnew(i) = xnew(i) - a*L
x(i) = x(i) - a*L
xold(i) = xold(i) - a*L
endif

if (xnew(i) .lt. 0.0) then
a = ceiling(xnew(i)/L)
xnew(i) = xnew(i) - a*L
x(i) = x(i) - a*L
xold(i) = xold(i) - a*L
endif

if (ynew(i) .gt. L) then
a = floor(ynew(i)/L)
ynew(i) = ynew(i) - a*L
y(i) = y(i) - a*L
yold(i) = yold(i) - a*L
endif

if (ynew(i) .lt. 0.0) then
a = ceiling(ynew(i)/L)
ynew(i) = ynew(i) - a*L
y(i) = y(i) - a*L
yold(i) = yold(i) - a*L
endif

enddo


! updates values

do i=1, N
xold(i) = x(i)
x(i) = xnew(i)
enddo

! store time and velocity data
do i=1, N
time(k) = k*dt
storex(k,i) = velx(i)
storey(k,i) = vely(i)
enddo

enddo

!
! calculating energy as function of time
!

do k=1, timesteps
energytemp = 0.0
cmxtemp = 0.0
cmytemp = 0.0
do i=1, N
energytemp = energytemp + 1.0/2.0*(velx(i)**2+vely(i)**2)
cmxtemp = cmxtemp + velx(i)
cmytemp = cmytemp + vely(i)
enddo
energy(k) = energytemp
velcm(k) = sqrt(cmxtemp**2 + cmytemp**2)
enddo!
! Outputting data
!

open (unit=9, file='positionData.txt', access='append', status='new')

do i=1, N
write(9,*) x(i), y(i)
9 continue
enddo

close(unit=9)
contains

!
! function to calculate the minimum distance on toroid
!

double precision function minDist(p,q) result(answer)

double precision, intent(in) :: p
double precision, intent(in) :: q

answer = min( abs(p-q),L-abs(p-q) )

end function minDist!
! function to calculate force between ith/jth particles
!

double precision function force(p,q,r,t) result(answer)

double precision, intent(in) :: p
double precision, intent(in) :: q
double precision, intent(in) :: r
double precision, intent(in) :: t

xdistance = minDist(p,q)
ydistance = minDist(r,t)
distance = sqrt(xdistance**2 + ydistance**2)

! if particles are too far, make force zero

if (distance .le. 3.0) then
answer = 24.0*(2.0/(distance**13)-1.0/(distance**7))
else
answer = 0.0
endif

end function force

end program