Фортран 90

program mol_dyn

implicit real(8) (a-h, o-z) 
integer, parameter :: ndim = 1000 
!
! Массивы координат,  проекций скоростей и ускорений частиц 
!
real(8), dimension(1:ndim) :: x, y, vx, vy, ax, ay
!
!  Задание начальной конфигурации частиц
!
call start(x, y, vx, vy, N, Sx, Sy, dt, dt2, nsnap, ntime) 
!
!  Вычисление ускорений частиц
!
call accel(x, y, ax, ay, N, Sx, Sy, zpe)

zpe =0.0d0
!
!  Энергия выводится через nsnap шагов
!
do isnap = 1, nsnap
!
! Внутренний цикл - по шагам по времени
!
do itime = 1, ntime
!
!  Выполняется сдвиг частиц
!
call move(x, y, vx, vy, ax, ay,N, Sx, Sy, dt, dt2, zke, zpe) 
end do
!
! Вывод значений энергии:  кинетической,  потенциальной и полной
!
call output(zke,   zpe,   Sx,   Sy,  dt,  N,  ntime)

end do

end program !
!====================================================
!
subroutine start(x,  y,  vx,  VY,  N,   Sx,   Sy,  dt,  dt2,  nsnap,  ntime)

implicit real(8) (a-h, o-z) integer,  parameter :: ndim = 1000
real(8), dimension(1:ndim) ::  x,  y, vx, vy, ax, ay
!
!  Число частиц 
!
n = 200
!
!  Размер области моделирования
!
sx = 100.0d0
sy = 100.0d0 !
! Шаг по времени
!
dt = 1.0d-11 dt2 = dt**2 !
! Максимальное значение скорости
!
vmax = 10.0d0 !
!   Число строк в таблице вывода и количество шагов по времени между ними
!
nsnap = 200
ntime = 500
!
!  Формируется начальная хаотическая конфигурация
!
do i = 1, N
call rndm(ranx) 
x(i) = sx * ranx 
call rndm(rany) 
y(i) = sy  * rany 
call rndm(ranx) 
call rndm(rany) 
vx(i) = vmax * (2  * ranx - 1) vy(i)  = vmax * (2 * rany - 1) end do

do i = 1, N
vxcum = vxcum + vx(i) vycum = vycum + vy(i) 
end do

vxcum = vxcum / N vycum = vycum / N

do i = 1, N
vx(i) = vx(i) - vxcum vy(i) = vy(i) - vycum 
end do

end subroutine

!
!====================================================
!
subroutine move(x, y, vx, vy, ax, ay, N, Sx, Sy, dt, dt2, zke, zpe)

implicit real(8) (a-h, o-z) 
integer, parameter :: ndim = 1000
real(8), dimension(1:ndim) ::  x,   y,   vx,   vy,   ax,   ay

do i = 1, N
xnew = x(i) + vx(i) * dt + 0.5d0 * ax(i) * dt2 
ynew = y(i) + vy(i) * dt + 0.5d0 * ay(i) * dt2
!
! Учет периодических граничных условий
!
call cellp(xnew,  ynew, vx(i), vy(i), sx, sy) 
x(i) = xnew 
y(i) = ynew
vx(i) = vx(i) + 0.5d0 * ax(i) * dt vy(i) = vy(i) + 0.5d0 * ay(i) * dt 
end do

call accel(x, y, ax,  ay,  N,   Sx,   Sy,   zpe)

do i = 1,  n
vx(i) = vx(i)  + 0.5d0  * dt * ax(i)
vy(i) = vy(i)  + 0.5d0  * dt * ay(i)
zke = zke + vx(i)**2 + vy(i)**2 end do

end subroutine
!
!====================================================
!
subroutine accel(x, Y, ax, ay, N, sx, sy, zpe)

implicit real(8) (a-h, o-z)
integer,  parameter  ::  ndim = 1000
real(8),  dimension(1:ndim) :: x, y, ax, ay

do i = 1,  n
 ax(i) = 0.0d0
ay(i) = 0.0d0 end do

do i = 1,  n 
do j = 1,  n 
 if(i /= j) 
  then 
  dx = x(i) - x(j)
dY = y(i) - y(j)

if(dabs(dx)>0.5 * sx) dx = dx - sign(sx, dx) if(dabs(dy)> 0.5 * sy)  dy = dy - sign(sy,  dy)

!  Вычисление силы,  действующей на j-ю частицу
!
r = dsqrt(dx**2 + dy**2) 
ri = 1.0d0 / r 
ri6 = ri**6
g = 24.0d0 * ri6 * ri * (2.0d0 * ri6 - 1)
force = 0.5d0 * g * ri
pot = 4.0d0  * ri6  * (ri6**2 - 1.0d0)

ax(i) = ax(i) +	force	* dx
ay(i) = ay(i) +	force	* dy
ax(j) = ax(j) -	force	* dx
ay(j) = ay(j) -	force	* dy
zpe = zpe + pot
end if

end do 
end do

end subroutine !
!====================================================
!
subroutine cellp(xnew,  ynew,  vx,  vy,   Sx,   Sy)

implicit real(8)   (a-h,   o-z)
if(xnew < 0) xnew = xnew + sx
if(xnew > sx) xnew = xnew - sx
if(ynew < 0) ynew = ynew + sy
if(ynew > sy)ynew = ynew - sy

end subroutine
!
!====================================================
!
subroutine output(zke, zpe, sx, sy,dt, n, ntime) 
implicit real(8)   (a-h,   o-z)
data iff/0/

if(iff == 0) then iff = 1
write(6,*) '	ke pe tot'
end if

zke = 0.5d0 * zke / ntime zpe = zpe  / ntime 
tot = zke + zpe

write(6, "(6(1x,  e13.6))") zke, zpe,tot

zke = 0.0d0
zpe = 0.0d0

end subroutine 
!
!  Генератор псевдослучайных чисел 
!
subroutine rndm(ran) implicit real*8(a-h,o-z)
data a31/2147483647.d0/,  mul/16807.d0/,  x/268435451.d0/ 
t = dmod(x * mul, a31) 
ra = t
x = t
ran = ra / a31 end subroutine