## Argon gas in a Box (torus)

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:

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-

 PhysOrg.com science news on PhysOrg.com >> Front-row seats to climate change>> Attacking MRSA with metals from antibacterial clays>> New formula invented for microscope viewing, substitutes for federally controlled drug
 Final shameless bump before abandoning hope for a reply.

Mentor

## Argon gas in a Box (torus)

It can take more than 94 minutes for an answer.

 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