!========================================================
! 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