Программа решения на языке программирования Фортран-4

 

PROGRAM md

DIMENSION x(2000),y(2000),z(2000),vx(2000),vy(2000),vz(2000)

DIMENSION ax(2000),ay(2000),az(2000)

character cnfile*30

data dt3 /0/

write(*,'(" enter end data file name ")')

read (*,'(A)') cnfile

open(2,file = cnfile,status = 'unknown')

CALL start(x,y,z,vx,vy,vz,N,Sx,Sy,Sz,dt,dt2,nsnap,ntime)

CALL accel(x,y,z,ax,ay,az,N,Sx,Sy,Sz,virial,zpe)

virial =0.0

zpe =0.0

DO 100 isnap = 1,nsnap

DO 10 itime = 1,ntime

CALL move(x,y,z,vx,vy,vz,ax,ay,az,N,Sx,Sy,Sz,dt,dt2,flx,fly,flz,virial,zke,zpe)

10 CONTINUE

dt3 = dt3 + dt

CALL output(flx,fly,flz,virial,zke,zpe,Sx,Sy,Sz,dt,N,ntime, dt3)

100 CONTINUE

close(2)

STOP

END

! подпрограмма

SUBROUTINE start(x,y,z,vx,vy,vz,N,Sx,Sy,Sz,dt,dt2,nsnap,ntime)

DIMENSION x(2000),y(2000),z(2000),vx(2000),vy(2000),vz(2000)

WRITE(6,*) 'Input number of particles'!'введите число частиц'

READ(5,*) N

WRITE(6,*) 'Input box lenght: Sx, Sy, Sz'!'введите размеры ящика'

READ(5,*) Sx,Sy,Sz

WRITE(6,*) 'Input time step: dt'!'введите шаг по времени'

READ(5,*) dt

dt2 = dt*dt

WRITE(6,*) 'Input max. velocity: vmax'!'введите максимальное значение скорости'

READ(5,*) vmax

WRITE(6,*) 'Input number of snapshots and number of steps between their: nsnap and ntime'!'введите число снимков и число шагов между ними'

READ(5,*) nsnap,ntime

WRITE(6,*) 'Input random negative number: iseed'!'введите отрицательное случайное начальное число'

READ(5,*) iseed

WRITE(6,*) 'Input number of configation: 1=old, 2=hot, 3=could'!'введите конфигурацию: 1=старая, 2=горячая, 3=хоподная'

READ(5,*) icf

IF (icf.eq.1) then

! чтение старой конфигурации

DO 10 i=1,N

READ(8,12) x(i),y(i),z(i),vx(i),vy(i),vz(i)

12 FORMAT(4(2x,f10.5))

10 CONTINUE

ELSEIF (icf.eq.3) then

! упорядоченная (холодная) начальная конфигурация

area1 = Sx*Sy*Sz/N

ys = 0.5*sqrt(3.0)

a = sqrt(area1/ys)

Ly = 2*int(0.5*(1.0 + Sy/(a*ys)))

Lz = 2*int(0.5*(1.0 + Sy/(a*ys)))

Lx = N/(Ly*Lz)

DO 30 ix = 1,Lx

DO 20 iy = 1,Ly

do 19 iz = 1,Lz

z(i) = (iz - 0.5)*a*ys

i=(iy - 1)*Ly + ix

y(i) = (iy - 0.5)*a*ys

IF (mod(iy,2).eq.0) then

x(i) = (ix - 0.25)*a

ELSE

x(i) = (ix - 0.75)*a

END IF

vx(i) = vmax*(2*ran(iseed) - 1)

vy(i) = vmax*(2*ran(iseed) - 1)

vz(i) = vmax*(2*ran(iseed) - 1)

19 continue

20 CONTINUE

30 CONTINUE

ELSE

! случайная (горячая) конфигурация

DO 40 i = 1,N

x(i) = Sx*ran(iseed)

y(i) = Sy*ran(iseed)

z(i) = Sz*ran(iseed)

vx(i) = vmax*(2*ran(iseed) - 1)

vy(i) = vmax*(2*ran(iseed) - 1)

vz(i) = vmax*(2*ran(iseed) - 1)

40 CONTINUE

ENDIF

DO 50 i = 1,N

vxcum = vxcum + vx(i)

vycum = vycum + vy(i)

vzcum = vzcum + vz(i)

50 CONTINUE

vxcum = vxcum/N

vycum = vycum/N

vzcum = vzcum/N

DO 60 i = 1,N

vx(i) = vx(i) - vxcum

vy(i) = vy(i) - vycum

vz(i) = vz(i) - vzcum

60 CONTINUE

RETURN

END

 

Подпрограмма MOVE

 

SUBROUTINE move (x,y,z,vx,vy,vz,ax,ay,az,N,Sx,Sy,Sz,dt,dt2,flx,fly,flz,virial,zke,zpe)

DIMENSION x(2000),y(2000),z(2000),vx(2000),vy(2000),vz(2000)

DIMENSION ax(2000),ay(2000),az(2000)

DO 10 i = 1,N

xnew = x(i) + vx(i)*dt + 0.5*ax(i)*dt2

ynew = y(i) + vy(i)*dt + 0.5*ay(i)*dt2

znew = z(i) + vz(i)*dt + 0.5*az(i)*dt2

CALL cellp(xnew,ynew,znew,vx(i),vy(i),vz(i),Sx,Sy,Sz,flx,fly,flz)

x(i) = xnew

y(i) = ynew

z(i) = znew

vx(i) = vx(i) + 0.5*ax(i)*dt

vy(i) = vy(i) + 0.5*ay(i)*dt

vz(i) = vz(i) + 0.5*az(i)*dt

10 CONTINUE

CALL accel(x,y,z,ax,ay,az,N,Sx,Sy,SZ,virial,zpe)

DO 20 i = 1,n

vx(i) = vx(i) + 0.5*dt*ax(i)

vy(i) = vy(i) + 0.5*dt*ay(i)

vz(i) = vz(i) + 0.5*dt*az(i)

zke = zke + vx(i)*vx(i) + vy(i)*vy(i) + vz(i)*vz(i)

virial = virial + ax(i)*x(i) + ay(i)*y(i) + az(i)*z(i)

20 CONTINUE

RETURN

END

Подпрограмма accel

 

SUBROUTINE accel(x,y,z,ax,ay,az,N,Sx,Sy,Sz,virial,zpe)

DIMENSION x(2000),y(2000),z(2000),ax(2000),ay(2000),az(2000)

DO 1 i = 1,n

ax(i) =0.0

ay(i) =0.0

az(i) =0.0

1 CONTINUE

DO 20 i = 1,(n-1)

DO 10 j = (i+1),n

dx = x(i) - x(j)

dy = y(i) - y(j)

dz = z(i) - z(j)

CALL sep(dx,dy,dz,Sx,Sy,Sz)

r = sqrt(dx*dx + dy*dy + dz*dz)

CALL FandU(r,force,pot)

ax(i) = ax(i) + force*dx

ay(i) = ay(i) + force*dy

az(i) = az(i) + force*dz

virial = virial + force*r*r

zpe = zpe + pot

! определение силы, действующей на частицу j по 3-му закону Ньютона

ax(j) = ax(j) - force*dx

ay(j) = ay(j) - force*dy

az(j) = az(j) - force*dz

10 CONTINUE

20 CONTINUE

RETURN

END

 

Подпрограмма FandU

 

SUBROUTINE FandU(r,force,pot)

ri = 1/r

ri3 = ri*ri*ri

ri6 = ri3*ri3

ri12 = ri6*ri6

g = 24*ri*ri6*(2*ri6 - 1)

force = g*ri

pot = 4*(ri12 - ri6)

RETURN

END

 

SUBROUTINE sep(dx,dy,dz,Sx,Sy,Sz)

IF (abs(dx).gt.0.5*Sx) dx = dx - sign(Sx,dx)

IF (abs(dy).gt.0.5*Sy) dy = dy - sign(Sy,dy)

IF (abs(dz).gt.0.5*Sz) dz = dz - sign(Sz,dz)

RETURN

END

 

SUBROUTINE cell(xnew,ynew,znew,Sx,Sy,Sz)

IF (xnew.lt.0) xnew = xnew + Sx

IF (xnew.gt.Sx) xnew = xnew - Sx

IF (ynew.lt.0) ynew = ynew + Sy

IF (ynew.gt.Sy) ynew = ynew - Sy

IF (znew.lt.0) znew = znew + Sz

IF (znew.gt.Sz) znew = znew - Sz

RETURN

END

 

SUBROUTINE cellp(xnew,ynew,znew,vx,vy,vz,Sx,Sy,Sz,flx,fly,flz)

IF (xnew.lt.0) then

xnew = xnew + Sx

flx = flx - vx

END IF

IF (xnew.gt.Sx) then

xnew = xnew - Sx

flx = flx + vx

END IF

IF (ynew.lt.0) then

ynew = ynew + Sy

fly = fly - vy

END IF

IF (ynew.gt.Sy) then

ynew = ynew - Sy

fly = fly + vy

END IF

IF (znew.lt.0) then

znew = znew + Sz

flz = flz - vz

END IF

IF (znew.gt.Sz) then

znew = znew - Sz

flz = flz + vz

END IF

RETURN

END

 

SUBROUTINE output(flx,fly,flz,virial,zke,zpe,Sx,Sy,Sz,dt,N,ntime,dt3)

data iff /0/

IF (iff.eq.0) then

iff = 1

write(6,*)' ke pe tot pflux pvirial pideal'

END IF

pflux = ((flx/Sx) + (fly/Sy) + (flz/Sz))/(3*dt*ntime)

zke = 0.5*zke/ntime

zpe = zpe/ntime

tot = zke + zpe

pideal = zke/(Sx*Sy*Sz)

pvirial = pideal + (0.5/(Sx*Sy*Sz))*virial/ntime

write(6,13) zke,zpe,tot,pflux, pvirial, pideal

write(2,14) dt3,zpe

13 format(6(1x,e12.4))

14 format(6(1x,e11.4))

zke = 0

zpe = 0

fly = 0

flx = 0

flz = 0

virial = 0

RETURN

END