Try   HackMD

Soft Disk Lennard-Jones 2D fluid

In this proble we try to simulate soft-disk fluid that consist of small particle interact with each other through Lennard-Jones potential in which are repulsive at short distance and attractive at longer distance.

(1)V(r)=4ϵ[(σr)12(σr)6]
And the force is the gradient of the potential plus a minus sign
(2)F(r)=24ϵrσr6[2(σr)61]

Verlet and Predictor-Corrector method

There are two easy way that we apply to solve the equation of motion by this force, since this system required relatively large amount of particle so we need a method that is both efficiency and accuracy. Both verlet and predictor-corrector are the two method that is on the top of the shelf. We first use velocity verlet method and did convergence test on both energy and momentum. The result are pretty decent, as the time step reduce the value converges and both energy and momentum have
relatively small error at

dt=104.
Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

We then take a step further by apply predictor-corrector method in which the central idea is that at each step it try to converge the difference between predictor and corrector until it reach a desire tolerance. Below figure show that the difference between predictor and corrector after at each time of correction. The difference decrease down to machine accuracy after roughly 8 iteration. So in our program, there are two criteria that we can consider the corrector are good enough: first is it need to at least iterate over at least 4 times, secondly it the difference needs to be less than the tolerance. If the the value didn't converges after 8 iteration, we will stop the iteration.

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Similarly, we did the convergence test with predictor-corrector method using difference tolerance which show a better accuracy with the same time step compare to verlet method. Follow table summarize the pros and cons of two method.

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

As you can see from above two figure about convergence of two method, I also show the temperature of the system. After a while the temperature approach equilibrium. I then measure the probability distribution of the speed

P(vx),
P(vy)
and
P(v)
at different temperature. Below three column are different initial velocity which result to different equilibrium temperature.

Periodic Boundary Condition

Since we are considering a infinite 2d system, we can use peridic boundary condition to achieve that. When particle hit one side of the wall, it will reappear on the other side; and the force we need to calculate the nearest force considering the image of periodic boundary condition.

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Maxwell-Boltzmann distribution

As you can see, according to Maxwell-Boltzmann distribution when the temperature are lower the distribution are more center at low speed with a higher peak, but if the temperature increase the distribution spreads out and the height also decrease. Since our system are smaller (only 49 particles), the distribution looks very discrete. If we increase the number of particles, the equilibrium temperature increase and the graph become smoother. Since the temperature also higher the distribution didn't seem much difference, but we can still see that the peak of each distribution shift to higher velocity.

NP=49
Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

NP=100
Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Temperature and Pressure

The formula for measuring temperature and pressure are given below.

(3)kBT=mdv2
and
(4)P=1A[NkBT12k=1Nj<kdV(r)drrjk]

Using above two equation we can monitor system temperature and
pressure at all time. Below show that after some time both
temperature and pressure goes to equilibrium.
NP=49

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

NP=100
Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Heat Capacity

If we then measure the energy after equilibrium and then by varying the initial velocity using following formula we can get roughly the desire equilibrium temperature.

(5)v=vTdesireTold
And the heat capacity
cv
is given by
dEdT
. Below left figure show the
E(t)
relation. The energy is negative is because it is govern by potential energy. In the case of
NP=49
the heat capacity is roughly 58.
Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

There is an alternatively way of calculating heat capacity which is given by

(6)T2T2=N(kBT)2[1NkBcv].
Using this formula we then measure the same system for heat capacity, the heat capacity are different by some degree when
NP=100
. But if we look at same method for two different number
of particles, you can see that
cv
roughly double as the particle number double which is what we expected.
NP=49

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

NP=100
Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Pair correlation of LJ fluid

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Code

!======================================================== ! purpose: MD simulation ! ! methods: velocity verlet + predictor-corrector ! ! input: NP, L, maxt, vmax, h ! ! output: position, velocity, T, K, U, P, Cv ! ! date programer description of change ! ==== ========= ===================== ! 11/19/20 morrisH original code !======================================================== program MD implicit none integer :: NP, i, n, j double precision :: L, maxt, h, time, tol double precision, dimension(:,:), allocatable :: r, v, rbef real :: start, finish double precision :: T, vmax, P double precision :: TF ! parameter time = 0.d0 NP = 100 L = 8.d0 maxt = 2.d-1 h = 1.d-4 tol = 1.d-17 vmax = 0.d0 TF = 1000.d0 n = nint(maxt / h) allocate(r(NP, 2)) allocate(rbef(NP, 2)) allocate(v(NP, 2)) call init(NP, L, r, v, vmax) ! call output(NP, L, r, v, time, T, P) call cpu_time(start) ! move first step rbef(i, :) = r(i, :) call verlet(NP, L, r, v, h) do do i = 1, n ! call verlet(NP, L, r, v, h) call predcorr(NP, L, r, v, h, rbef, tol) ! call output(NP, L, r, v, time, T, P) call pstatus(n, i, time, T, P) time = time + h end do call output(NP, L, r, v, time, T, P) ! rescale velocity do j = 1, NP v(j, :) = dsqrt((T + 1.d1) / T) * v(j, :) end do if (T.gt.TF) then; exit; end if end do call cpu_time(finish) write(*, *) 'cost :', finish-start end program MD subroutine init(NP, L, r, v, vmax) implicit none integer :: NP, i, j, num double precision :: L, r(NP, 2), v(NP, 2), vmax, dx, dy call random_seed() num = int(sqrt(dble(NP))) dx = L / (dble(num) + 1.d0) dy = L / (dble(num) + 1.d0) ! random position do i = 1, num do j = 1, num r(num*(i-1)+j, 1) = dble(i) * dx r(num*(i-1)+j, 2) = dble(j) * dy end do end do ! random velocity call random_number(v) v = (v - .5d0) * 2.d0 * vmax end subroutine init subroutine output(NP, L, r, v, time, T, P) implicit none integer :: NP, i, j double precision :: r(NP, 2), v(NP, 2), L, K, U, E, Px, Py, T, P, Cv double precision :: r1(2), r2(2), d, dx, dy, time double precision :: TT ! conpute some physical quantity K = 0.d0 U = 0.d0 Px = 0.d0 Py = 0.d0 T = 0.d0 TT = 0.d0 P = 0.d0 do i = 1, NP r1 = r(i, :) K = K + 0.5d0 * (v(i, 1)**2 + v(i, 2)**2) Px = Px + v(i, 1) Py = Py + v(i, 2) T = T + (v(i, 1)**2 + v(i, 2)**2) TT = TT + (v(i, 1)**2 + v(i, 2)**2)**2 ! calculate potential energy do j = 1, i if (i.ne.j) then r2 = r(j, :) d = dsqrt((r1(1) - r2(1))**2 + (r1(2) - r2(2))**2) if (d.gt.(L/2.d0)) then ! find nearest image dx = r2(1) - r1(1) dy = r2(2) - r1(2) if (dx.gt.(L/2.d0)) then; dx = dx - L; end if if (dy.gt.(L/2.d0)) then; dy = dy - L; end if if (dx.lt.(-L/2.d0)) then; dx = dx + L; end if if (dy.lt.(-L/2.d0)) then; dy = dy + L; end if d = dsqrt(dx ** 2 + dy ** 2) U = U + 4.d0 * ( (1.d0 / d)**12 - (1.d0 / d)*6 ) P = P + 0.5d0 * (24.d0 / d) * (2.d0 * (1 / d)**12 - (1 / d)**6) * d else U = U + 4.d0 * ( (1.d0 / d)**12 - (1.d0 / d)*6 ) P = P + 0.5d0 * (24.d0 / d) * (2.d0 * (1 / d)**12 - (1 / d)**6) * d end if end if end do end do T = T / NP TT = TT / NP E = K + U P = (P + NP * T) / L**2 Cv = NP / (1.d0 - (TT - T**2) / (NP * T**2)) ! open(66, file='data/np49/Cv/r.dat', access='append', status='unknown') ! write(66, '(100e40.30)') r(:, 1) ! write(66, '(100e40.30)') r(:, 2) ! close(66) ! open(77, file='data/np49/Cv/v.dat', access='append', status='unknown') ! write(77, '(100e40.30)') v(:, 1) ! write(77, '(100e40.30)') v(:, 2) ! close(77) open(66, file='data/np100/Cv/info_alt.dat', access='append', status='unknown') write(66, '(9e40.30)') time, K, U, E, Px, Py, T, P, Cv close(66) end subroutine output subroutine verlet(NP, L, r, v, h) implicit none integer :: NP, i, j double precision :: L, r(NP, 2), v(NP, 2), rnew(NP, 2), vnew(NP, 2) double precision :: a1(2), a2(2), d, dx, dy, r1(2), r2(2) double precision :: h do i = 1, NP r1 = r(i, :) ! compute a(t) a1 = 0.d0 do j = 1, NP if (i.ne.j) then r2 = r(j, :) d = dsqrt((r1(1) - r2(1))**2 + (r1(2) - r2(2))**2) if (d.gt.(L/2.d0)) then ! find nearest image dx = r2(1) - r1(1) dy = r2(2) - r1(2) if (dx.gt.(L/2.d0)) then; dx = dx - L; end if if (dy.gt.(L/2.d0)) then; dy = dy - L; end if if (dx.lt.(-L/2.d0)) then; dx = dx + L; end if if (dy.lt.(-L/2.d0)) then; dy = dy + L; end if d = dsqrt(dx ** 2 + dy ** 2) a1 = a1 + (24.d0 / d) * (2.d0 * (1 / d)**12 - (1 / d)**6) * & (/ -dx, -dy /) / d else a1 = a1 + (24.d0 / d) * (2.d0 * (1 / d)**12 - (1 / d)**6) * & (r1 - r2) / d end if end if end do rnew(i, :) = r(i, :) + v(i, :) * h + .5d0 * a1 * h ** 2 ! compute a(t+h) a2 = 0.d0 r1 = rnew(i, :) do j = 1, NP if (i.ne.j) then r2 = r(j, :) d = dsqrt((r1(1) - r2(1))**2 + (r1(2) - r2(2))**2) if (d.gt.(L/2.d0)) then ! find nearest image dx = r2(1) - r1(1) dy = r2(2) - r1(2) if (dx.gt.(L/2.d0)) then; dx = dx - L; end if if (dy.gt.(L/2.d0)) then; dy = dy - L; end if if (dx.lt.(-L/2.d0)) then; dx = dx + L; end if if (dy.lt.(-L/2.d0)) then; dy = dy + L; end if d = dsqrt(dx ** 2 + dy ** 2) a2 = a2 + (24.d0 / d) * (2.d0 * (1 / d)**12 - (1 / d)**6) * & (/ -dx, -dy /) / d else a2 = a2 + (24.d0 / d) * (2.d0 * (1 / d)**12 - (1 / d)**6) * & (r1 - r2) / d end if end if end do vnew(i, :) = v(i, :) + h / 2.d0 * (a1 + a2) ! periodic bd if (rnew(i, 1) > L) then rnew(i, 1) = rnew(i, 1) - L elseif (rnew(i, 1) < 0.d0) then rnew(i, 1) = rnew(i, 1) + L end if if (rnew(i, 2) > L) then rnew(i, 2) = rnew(i, 2) - L elseif (rnew(i, 2) < 0.d0) then rnew(i, 2) = rnew(i, 2) + L end if end do r = rnew v = vnew end subroutine verlet subroutine predcorr(NP, L, r, v, h, rbef, tol) implicit none integer :: NP, i, j, k double precision :: L, r(NP, 2), v(NP, 2), rnew(NP, 2), vnew(NP, 2) double precision :: a1(2), a2(2), d, dx, dy, r1(2), r2(2) double precision :: h double precision :: rbef(NP, 2), rp(2), rc(2), vc(2), diff, tol do i = 1, NP ! predictor rp = rbef(i, :) + 2.d0 * v(i, :) * h r1 = r(i, :) ! compute a(t) a1 = 0.d0 do j = 1, NP if (i.ne.j) then r2 = r(j, :) d = dsqrt((r1(1) - r2(1))**2 + (r1(2) - r2(2))**2) if (d.gt.(L/2.d0)) then ! find nearest image dx = r2(1) - r1(1) dy = r2(2) - r1(2) if (dx.gt.(L/2.d0)) then; dx = dx - L; end if if (dy.gt.(L/2.d0)) then; dy = dy - L; end if if (dx.lt.(-L/2.d0)) then; dx = dx + L; end if if (dy.lt.(-L/2.d0)) then; dy = dy + L; end if d = dsqrt(dx ** 2 + dy ** 2) a1 = a1 + (24.d0 / d) * (2.d0 * (1 / d)**12 - (1 / d)**6) * & (/ -dx, -dy /) / d else a1 = a1 + (24.d0 / d) * (2.d0 * (1 / d)**12 - (1 / d)**6) * & (r1 - r2) / d end if end if end do ! calculate a^p k = 0 do k = k + 1 ! compute a(t+h) a2 = 0.d0 r1 = rp do j = 1, NP if (i.ne.j) then r2 = r(j, :) d = dsqrt((r1(1) - r2(1))**2 + (r1(2) - r2(2))**2) if (d.gt.(L/2.d0)) then ! find nearest image dx = r2(1) - r1(1) dy = r2(2) - r1(2) if (dx.gt.(L/2.d0)) then; dx = dx - L; end if if (dy.gt.(L/2.d0)) then; dy = dy - L; end if if (dx.lt.(-L/2.d0)) then; dx = dx + L; end if if (dy.lt.(-L/2.d0)) then; dy = dy + L; end if d = dsqrt(dx ** 2 + dy ** 2) a2 = a2 + (24.d0 / d) * (2.d0 * (1 / d)**12 - (1 / d)**6) * & (/ -dx, -dy /) / d else a2 = a2 + (24.d0 / d) * (2.d0 * (1 / d)**12 - (1 / d)**6) * & (r1 - r2) / d end if end if end do ! corrector vc = v(i, :) + .5d0 * h * (a1 + a2) rc = r(i, :) + .5d0 * h * (v(i, :) + vc) ! compute error diff = dsqrt( (rc(1) - rp(1))**2 + (rc(2) - rp(2))**2 ) ! output diff ! open(66, file='data/diff.dat', access='append', status='unknown') ! write(66, '(2e40.30)') dble(k), diff ! close(66) ! print *, k, diff if ((diff.le.tol).or.(k.ge.8)) then exit else rp = rc end if end do rnew(i, :) = rc vnew(i, :) = vc ! periodic bd if (rnew(i, 1) > L) then rnew(i, 1) = rnew(i, 1) - L elseif (rnew(i, 1) < 0.d0) then rnew(i, 1) = rnew(i, 1) + L end if if (rnew(i, 2) > L) then rnew(i, 2) = rnew(i, 2) - L elseif (rnew(i, 2) < 0.d0) then rnew(i, 2) = rnew(i, 2) + L end if end do rbef = r r = rnew v = vnew end subroutine predcorr subroutine pstatus(n, i, time, T, P) implicit none integer :: n, i double precision :: time, T, P write(*, '(1a, 1f6.2, 1a5, 1a5, 1f4.1, 1a, 1a, 1a5, 1f6.2, 1a5, 1f8.2)') & 't:', time, '', 'left:', dble(n-i)/dble(n)*100.d0, '%', & '', 'T=', T, 'P=', P end subroutine pstatus