- #1
llello
- 32
- 0
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-
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-