半步蛙跳算法程序

C *******************************************************************
C ** THIS FORTRAN CODE IS INTENDED TO ILLUSTRATE POINTS MADE IN **
C ** THE TEXT. TO OUR KNOWLEDGE IT WORKS CORRECTLY. HOWEVER IT IS **
C ** THE RESPONSIBILITY OF THE USER TO TEST IT, IF IT IS USED IN A **
C ** RESEARCH APPLICATION. **
C *******************************************************************

C *******************************************************************
C ** FICHE F.3 **
C ** TWO SEPARATE PARTS: FORTRAN AND BASIC VERSIONS. **
C *******************************************************************



C *******************************************************************
C ** FICHE F.3 - PART A **
C ** FORTRAN PROGRAM USING THE LEAPFROG ALGORITHM. **
C *******************************************************************



PROGRAM FROGGY

COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ

C *******************************************************************
C ** FORTRAN PROGRAM TO CONDUCT MOLECULAR DYNAMICS OF ATOMS. **
C ** **
C ** A SPECIAL LOW-STORAGE FORM OF THE LEAPFROG VERLET ALGORITHM **
C ** IS USED AND WE TAKE THE LENNARD-JONES POTENTIAL AS AN EXAMPLE.**
C ** **
C ** REFERENCE: **
C ** **
C ** FINCHAM AND HEYES, CCP5 QUARTERLY, 6, 4, 1982. **
C ** **
C ** ROUTINES REFERENCED: **
C ** **
C ** SUBROUTINE READCN ( CNFILE ) **
C ** READS IN CONFIGURATION **
C ** SUBROUTINE FORCE ( DT, SIGMA, RCUT, NEWV, NEWVC, NEWW ) **
C ** CALCULATES THE ACCELERATIONS AND ADVANCES THE VELOCITIES **
C ** FROM T - DT/2 TO T + DT/2. ALSO CALCULATES POTENTIAL **
C ** ENERGY AND VIRIAL AT TIMESTEP T. **
C ** SUBROUTINE MOVE ( DT ) **
C ** ADVANCES POSITIONS FROM T TO T + DT. **
C ** SUBROUTINE KINET ( NEWK ) **
C ** CALCULATES KINETIC ENERGY. **
C ** SUBROUTINE WRITCN ( CNFILE ) **
C ** WRITES OUT CONFIGURATION **
C ** **
C ** PRINCIPAL VARI

ABLES: **
C ** **
C ** INTEGER N NUMBER OF ATOMS **
C ** REAL DT TIMESTEP **
C ** REAL RX(N),RY(N),RZ(N) ATOMIC POSITIONS **
C ** REAL VX(N),VY(N),VZ(N) ATOMIC VELOCITIES **
C ** REAL ACV,ACK ETC. AVERAGE VALUE ACCUMULATORS **
C ** REAL AVV,AVK ETC. AVERAGE VALUES **
C ** REAL ACVSQ, ACKSQ ETC. AVERAGE SQUARED VALUE ACCUMULATORS **
C ** REAL FLV,FLK ETC. FLUCTUATION AVERAGES **
C ** **
C ** USAGE: **
C ** **
C ** THE LEAPFROG ALGORITHM IS OF THE FORM **
C ** VX(T + DT/2) = VX(T - DT/2) + DT * AX(T) (SIMILARLY FOR Y,Z) **
C ** RX(T + DT) = RX(T) + DT * VX(T + DT/2) (SIMILARLY FOR Y,Z) **
C ** TO SAVE STORAGE IN THIS PROGRAM WE ACCUMULATE VALUES AX(T) **
C ** DIRECTLY ONTO THE VELOCITIES VX(T - DT/2) IN THE FORCE LOOP. **
C ** THIS MEANS THAT AN APPROXIMATION MUST BE USED FOR THE KINETIC **
C ** ENERGY AT EACH TIME STEP: **
C ** K = ( OLDK + NEWK ) / 2 + ( OLDV - 2 * V + NEWV ) / 8 **
C ** WHERE K = K( T ), V = V( T ) **
C ** OLDK = K( T - DT/2 ), OLDV = V( T - DT ) **
C ** NEWK = K( T + DT/2 ), NEWV = V( T + DT ) **
C ** AT THE START OF A STEP THE FOLLOWING VARIABLES ARE STORED: **
C ** R : R(STEP) V : V(STEP+1/2) **
C ** OLDV : V(STEP-2) V : V(STEP-1) NEWV : V(STEP) **
C ** OLDK : K(STEP-1/2) NEWK : K(STEP+1/2) **
C ** THIS PROGRAM USES UNIT 10 FOR CONFIGURATION INPUT AND OUTPUT **
C ** **
C ** UNITS: **
C ** **
C ** THE PROGRAM USES LENNARD-JONES REDUCED UNITS FOR USER INPUT **
C ** AND OUTPUT BUT CONDUCTS SIMULATION IN A BOX OF UNIT LENGTH. **
C ** SUMMARY FOR BOX LENGTH L, ATOMIC MASS M, AND LENNARD-JONES **
C ** POTENTIAL PARAMETERS SIGMA AND EPSILON: **
C ** OUR PROGRAM LENNARD-JONES SYSTEM **
C ** LENGTH L SIGMA **
C ** MASS M M **
C ** ENERGY EPSILON EPSILON **
C ** TIME SQRT(M*L**2/EPSILON) SQRT(M*SIGMA**2/E

PSILON) **
C ** VELOCITY SQRT(EPSILON/M) SQRT(EPSILON/M) **
C ** PRESSURE EPSILON/L**3 EPSILON/SIGMA**3 **
C *******************************************************************

INTEGER N
PARAMETER ( N = 108 )

REAL RX(N), RY(N), RZ(N)
REAL VX(N), VY(N), VZ(N)
INTEGER STEP, NSTEP, IPRINT
REAL DENS, NORM, RCUT, DT, SIGMA
REAL V, K, E, W
REAL OLDK, NEWK, OLDV, NEWV, NEWW
REAL OLDVC, NEWVC, VC, KC, EC
REAL VN, KN, EN, ECN, PRES, TEMP
REAL ACV, ACK, ACE, ACEC, ACP, ACT
REAL AVV, AVK, AVE, AVEC, AVP, AVT
REAL ACVSQ, ACKSQ, ACESQ, ACECSQ, ACPSQ, ACTSQ
REAL FLV, FLK, FLE, FLEC, FLP, FLT
REAL SR3, SR9, VLRC, WLRC, PI, SIGCUB
CHARACTER CNFILE*30, TITLE*80
REAL FREE

PARAMETER ( FREE = 3.0 )
PARAMETER ( PI = 3.1415927 )

C *******************************************************************

C ** READ IN INITIAL PARAMETERS **

WRITE(*,'(1H1,'' **** PROGRAM FROGGY **** '')')
WRITE(*,'(//, '' MOLECULAR DYNAMICS OF LENNARD-JONES ATOMS '')')
WRITE(*,'( '' LEAPFROG ALGORITHM WITH MINIMAL STORAGE '')')

WRITE(*,'('' ENTER RUN TITLE '')')
READ (*,'(A)') TITLE
WRITE(*,'('' ENTER NUMBER OF STEPS '')')
READ (*,*) NSTEP
WRITE(*,'('' ENTER NUMBER OF STEPS BETWEEN OUTPUT LINES '')')
READ (*,*) IPRINT
WRITE(*,'('' ENTER CONFIGURATION FILENAME '')')
READ (*,'(A)') CNFILE
WRITE(*,'('' ENTER THE FOLLOWING IN LENNARD-JONES UNITS '')')
WRITE(*,'('' ENTER DENSITY '')')
READ (*,*) DENS
WRITE(*,'('' ENTER POTENTIAL CUTOFF DISTANCE '')')
READ (*,*) RCUT
WRITE(*,'('' ENTER TIME STEP '')')
READ (*,*) DT

WRITE(*,'(//1X,A)') TITLE
WRITE(*,'('' NUMBER OF ATOMS = '',I10 )') N
WRITE(*,'('' NUMBER OF STEPS = '',I10 )') NSTEP
WRITE(*,'('' OUTPUT FREQUENCY = '',I10 )') IPRINT
WRITE(*,'('' POTENTIAL CUTOFF = '',F10.4)') RCUT
WRITE(*,'('' DENSITY = '',F10.4)') DENS
WRITE(*,'('' TIME STEP = '',F10.6)') DT

C ** READ CONFIGURATION INTO COMMON / BLOCK1 / VARIABLES **

CALL READCN ( CNFILE )

C ** CONVERT INPUT DATA TO PROGRAM UNITS **

SIGMA = ( DENS / REAL ( N ) ) ** ( 1.0 / 3.0 )
RCUT = RCUT * SIGMA
DT = DT * SIGMA
DENS = DENS / ( SIGMA ** 3 )

C ** CALCULATE LONG-RANGE CORRECTIONS **
C ** NOTE: SPECIFIC TO LENNARD-JONES **

SR3 = ( SIGMA / RCUT ) ** 3

SR9 = SR3 ** 3
SIGCUB = SIGMA ** 3
VLRC = ( 8.0 /9.0 ) * PI * DENS * SIGCUB * REAL ( N )
: * ( SR9 - 3.0 * SR3 )
WLRC = ( 16.0 / 9.0 ) * PI * DENS * SIGCUB * REAL ( N )
: * ( 2.0 * SR9 - 3.0 * SR3 )

C ** ZERO ACCUMULATORS **

ACV = 0.0
ACK = 0.0
ACE = 0.0
ACEC = 0.0
ACP = 0.0
ACT = 0.0

ACVSQ = 0.0
ACKSQ = 0.0
ACESQ = 0.0
ACECSQ = 0.0
ACPSQ = 0.0
ACTSQ = 0.0

FLV = 0.0
FLK = 0.0
FLE = 0.0
FLEC = 0.0
FLP = 0.0
FLT = 0.0

C ** CALCULATE INITIAL VALUES **

CALL FORCE ( -DT, SIGMA, RCUT, NEWV, NEWVC, NEWW )
CALL MOVE ( -DT )
CALL FORCE ( -DT, SIGMA, RCUT, V, VC, W )
CALL FORCE ( DT, SIGMA, RCUT, V, VC, W )
CALL KINET ( OLDK )
CALL MOVE ( DT )
CALL FORCE ( DT, SIGMA, RCUT, NEWV, NEWVC, NEWW )
CALL KINET ( NEWK )

C ** INCLUDE LONG-RANGE CORRECTIONS **

V = V + VLRC
W = W + WLRC
NEWV = NEWV + VLRC
NEWW = NEWW + WLRC

IF ( IPRINT .LE. 0 ) IPRINT = NSTEP + 1

WRITE(*,'(//1X,''**** START OF DYNAMICS ****'')')
WRITE(*,10001)

C *******************************************************************
C ** MAIN LOOP BEGINS **
C *******************************************************************

DO 1000 STEP = 1, NSTEP

C ** IMPLEMENT ALGORITHM **

CALL MOVE ( DT )

OLDV = V
V = NEWV
OLDVC = VC
VC = NEWVC
W = NEWW

CALL FORCE ( DT, SIGMA, RCUT, NEWV, NEWVC, NEWW )

C ** INCLUDE LONG-RANGE CORRECTIONS **

NEWV = NEWV + VLRC
NEWW = NEWW + WLRC

C ** CALCULATE KINETIC ENERGY AT CURRENT STEP **

K = ( NEWK + OLDK ) / 2.0
: + ( NEWV - 2.0 * V + OLDV ) / 8.0
KC = ( NEWK + OLDK ) / 2.0
: + ( NEWVC - 2.0 * VC + OLDVC ) / 8.0
OLDK = NEWK

CALL KINET ( NEWK )

C ** CALCULATE INSTANTANEOUS VALUES **

E = K + V
EC = KC + VC
VN = V / REAL ( N )
KN = K / REAL ( N )
EN = E / REAL ( N )
ECN = EC / REAL ( N )
TEMP = 2.0 * KN / FREE
PRES = DENS * TEMP + W

C ** CONVERT TO LENNARD-JONES UNITS **

PRES = PRES * SIGMA ** 3

C ** INCREMENT ACCUMULATORS **

ACE = ACE + EN
ACEC = ACEC + ECN
ACK = ACK + KN
ACV = ACV + VN
ACP = ACP + PRES

ACESQ = ACESQ + EN ** 2
ACECSQ = ACECSQ + ECN ** 2
ACKSQ = ACKSQ + KN ** 2
ACVSQ = ACVS

Q + VN ** 2
ACPSQ = ACPSQ + PRES ** 2

C ** OPTIONALLY PRINT INFORMATION **

IF ( MOD( STEP, IPRINT ) .EQ. 0 ) THEN

WRITE(*,'(1X,I8,6(2X,F10.4))')
: STEP, EN, ECN, KN, VN, PRES, TEMP
ENDIF

1000 CONTINUE

C *******************************************************************
C ** MAIN LOOP ENDS **
C *******************************************************************

WRITE(*,'(/1X,''**** END OF DYNAMICS **** ''//)')

C ** OUTPUT RUN AVERAGES **

NORM = REAL ( NSTEP )

AVE = ACE / NORM
AVEC = ACEC / NORM
AVK = ACK / NORM
AVV = ACV / NORM
AVP = ACP / NORM

ACESQ = ( ACESQ / NORM ) - AVE ** 2
ACECSQ = ( ACECSQ / NORM ) - AVEC ** 2
ACKSQ = ( ACKSQ / NORM ) - AVK ** 2
ACVSQ = ( ACVSQ / NORM ) - AVV ** 2
ACPSQ = ( ACPSQ / NORM ) - AVP ** 2

IF ( ACESQ .GT. 0.0 ) FLE = SQRT ( ACESQ )
IF ( ACECSQ .GT. 0.0 ) FLEC = SQRT ( ACECSQ )
IF ( ACKSQ .GT. 0.0 ) FLK = SQRT ( ACKSQ )
IF ( ACVSQ .GT. 0.0 ) FLV = SQRT ( ACVSQ )
IF ( ACPSQ .GT. 0.0 ) FLP = SQRT ( ACPSQ )

AVT = AVK * 2.0 / FREE
FLT = FLK * 2.0 / FREE

WRITE(*,'('' AVERAGES'',6(2X,F10.5))')
: AVE, AVEC, AVK, AVV, AVP, AVT
WRITE(*,'('' FLUCTS '',6(2X,F10.5))')
: FLE, FLEC, FLK, FLV, FLP, FLT

C ** WRITE OUT FINAL CONFIGURATION **

CALL WRITCN ( CNFILE )

STOP

10001 FORMAT(//1X,'TIMESTEP ..ENERGY.. CUTENERGY.'
: ' ..KINETIC. ..POTENT..',
: ' .PRESSURE. ..TEMPER..'/)
END



SUBROUTINE FORCE ( DT, SIGMA, RCUT, V, VC, W )

COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ

C *******************************************************************
C ** ROUTINE TO COMPUTE FORCES WITH CUTOFF & MINIMUM IMAGING. **
C ** **
C ** THE ROUTINE ACTUALLY RETURNS TWO DIFFERENT POTENTIAL ENERGIES.**
C ** V IS CALCULATED USING THE UNSHIFTED LENNARD-JONES POTENTIAL. **
C ** WHEN LONG-RANGE TAIL CORRECTIONS ARE ADDED, THIS MAY BE USED **
C ** TO CALCULATE THERMODYNAMIC INTERNAL ENERGY ETC. **
C ** VC IS CALCULATED USING THE SHIFTED LENNARD-JONES POTENTIAL **
C ** WITH NO DISCONTINUITY AT THE CUTOFF. THIS MAY BE USED TO **
C ** CHECK ENERGY CONSERVATION. **
C ** PROGRAM UNITS MAKE EPSILON = MASS = 1, BUT SIGMA IS NOT UNITY **
C ** SINCE WE TAKE A BOX OF UNIT LENGTH. **
C ** THE ROUTINE IS COMBINED WITH THE LEAPFROG ALGORITHM FOR **
C ** ADVANCING VELOCITIES, I.E. FORCES ARE ACCUMULATED DIRECTLY **
C ** O

NTO VELOCITIES. **
C *******************************************************************

INTEGER N
PARAMETER ( N = 108 )
REAL RX(N), RY(N), RZ(N)
REAL VX(N), VY(N), VZ(N)
REAL V, RCUT, DT, SIGMA, W, VC

INTEGER I, J, NCUT
REAL RCUTSQ, VCUT, SIGSQ
REAL RXI, RYI, RZI, VXI, VYI, VZI
REAL RXIJ, RYIJ, RZIJ, RIJSQ
REAL SR2, SR6, SR12
REAL VIJ, WIJ, VELIJ, DVX, DVY, DVZ

C *******************************************************************

C ** CALCULATE SQUARED DISTANCES FOR INNER LOOP **

SIGSQ = SIGMA ** 2
RCUTSQ = RCUT ** 2

C ** CONDUCT OUTER LOOP OVER ATOMS **

V = 0.0
W = 0.0
NCUT = 0

DO 1000 I = 1, N - 1

RXI = RX(I)
RYI = RY(I)
RZI = RZ(I)
VXI = VX(I)
VYI = VY(I)
VZI = VZ(I)

C ** CONDUCT INNER LOOP OVER ATOMS **

DO 900 J = I + 1, N

RXIJ = RXI - RX(J)
RXIJ = RXIJ - ANINT ( RXIJ )

IF ( ABS ( RXIJ ) .LT. RCUT ) THEN

RYIJ = RYI - RY(J)
RYIJ = RYIJ - ANINT ( RYIJ )
RIJSQ = RXIJ ** 2 + RYIJ ** 2

IF ( RIJSQ .LT. RCUTSQ ) THEN

RZIJ = RZI - RZ(J)
RZIJ = RZIJ - ANINT ( RZIJ )
RIJSQ = RIJSQ + RZIJ ** 2

IF ( RIJSQ .LT. RCUTSQ ) THEN

SR2 = SIGSQ / RIJSQ
SR6 = SR2 ** 3
SR12 = SR6 ** 2
VIJ = 4.0 * ( SR12 - SR6 )
WIJ = 24.0 * ( 2.0 * SR12 - SR6 )
VELIJ = WIJ * DT / RIJSQ
DVX = VELIJ * RXIJ
DVY = VELIJ * RYIJ
DVZ = VELIJ * RZIJ
VXI = VXI + DVX
VYI = VYI + DVY
VZI = VZI + DVZ
VX(J) = VX(J) - DVX
VY(J) = VY(J) - DVY
VZ(J) = VZ(J) - DVZ
V = V + VIJ
W = W + WIJ
NCUT = NCUT + 1

ENDIF

ENDIF

ENDIF

900 CONTINUE

C ** END OF INNER LOOP **

VX(I) = VXI
VY(I) = VYI
VZ(I) = VZI

1000 CONTINUE

C ** END OF OUTER LOOP **

C ** CALCULATE POTENTIAL SHIFT **

SR2 = SIGSQ / RCUTSQ
SR6 = SR2 * SR2 * SR2
SR12 = SR6 * SR6
VIJ = 4.0 * ( SR12 - SR6 )
VC = V - REAL ( NCUT ) * VIJ

W = W / 3.0

RETURN
END



SUBROUTI

NE READCN ( CNFILE )

COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ

C *******************************************************************
C ** SUBROUTINE TO READ IN INITIAL CONFIGURATION FROM UNIT 10 **
C *******************************************************************

INTEGER N
PARAMETER ( N = 108 )
CHARACTER CNFILE*(*)
REAL RX(N), RY(N), RZ(N)
REAL VX(N), VY(N), VZ(N)

INTEGER CNUNIT, NN
PARAMETER ( CNUNIT = 10 )

C ******************************************************************

OPEN ( UNIT = CNUNIT, FILE = CNFILE, STATUS = 'OLD',
: FORM = 'UNFORMATTED' )

READ ( CNUNIT ) NN
IF ( NN .NE. N ) STOP ' INCORRECT NUMBER OF ATOMS '
READ ( CNUNIT ) RX, RY, RZ
READ ( CNUNIT ) VX, VY, VZ

CLOSE ( UNIT = CNUNIT )

RETURN
END



SUBROUTINE WRITCN ( CNFILE )

COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ

C *******************************************************************
C ** ROUTINE TO WRITE OUT FINAL CONFIGURATION TO UNIT 10 **
C *******************************************************************

INTEGER N
PARAMETER ( N = 108 )
CHARACTER CNFILE*(*)
REAL RX(N), RY(N), RZ(N)
REAL VX(N), VY(N), VZ(N)

INTEGER CNUNIT
PARAMETER ( CNUNIT = 10 )

C ****************************************************************

OPEN ( UNIT = CNUNIT, FILE = CNFILE, STATUS = 'UNKNOWN',
: FORM = 'UNFORMATTED' )

WRITE ( CNUNIT ) N
WRITE ( CNUNIT ) RX, RY, RZ
WRITE ( CNUNIT ) VX, VY, VZ

CLOSE ( UNIT = CNUNIT )

RETURN
END



SUBROUTINE MOVE ( DT )

COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ

C *******************************************************************
C ** LEAPFROG VERLET ALGORITHM FOR ADVANCING POSITIONS **
C *******************************************************************

INTEGER N
PARAMETER ( N = 108 )
REAL RX(N), RY(N), RZ(N)
REAL VX(N), VY(N), VZ(N)
REAL DT

INTEGER I

C *******************************************************************

DO 1000 I = 1, N

RX(I) = RX(I) + VX(I) * DT
RY(I) = RY(I) + VY(I) * DT
RZ(I) = RZ(I) + VZ(I) * DT

1000 CONTINUE

RETURN
END



SUBROUTINE KINET ( K )

COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ

C *******************************************************************
C ** COMPUTES KINETIC ENERGY **
C *******************************************************************

INTEGER N
PARAMETER ( N = 108 )

REAL RX(N), RY(N), RZ(N)
REAL VX(N), VY(N), VZ(N)
REAL K

INTEGER I

C *******************************************************************

K = 0.0

DO 100 I = 1, N

K = K + VX(I) ** 2 + VY(I) ** 2 + VZ(I) ** 2

100 CONTINUE

K = K * 0.5

RETURN
END