|
older-version
|
|
a.tgz,
f.00,
f.1,
f.10,
f.11,
f.12,
f.13,
f.14,
f.15,
f.16,
f.17,
f.18,
f.19,
f.2,
f.20,
f.21,
f.22,
f.23,
f.24,
f.25,
f.26,
f.27,
f.28,
f.29,
f.3,
f.30,
f.31,
f.32,
f.33,
f.34,
f.35,
f.36,
f.37,
f.4,
f.5,
f.6,
f.7,
f.8,
f.9
|
|
|
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.12 **
C ** MONTE CARLO SIMULATION IN THE CONSTANT-NPT ENSEMBLE. **
C *******************************************************************
PROGRAM MCNPT
COMMON / BLOCK1 / RX, RY, RZ
C *******************************************************************
C ** MONTE CARLO SIMULATION IN THE CONSTANT-NPT ENSEMBLE. **
C ** **
C ** THIS PROGRAM TAKES A CONFIGURATION OF LENNARD JONES ATOMS **
C ** AND PERFORMS A MONTE CARLO SIMULATION AT CONSTANT NPT. THE **
C ** BOX IS IN UNITS OF SIGMA THE LENNARD JONES DIAMETER. **
C ** THERE ARE NO LOOKUP TABLES INCLUDED. **
C ** **
C ** REFERENCE: **
C ** **
C ** MCDONALD, CHEM. PHYS. LETT. 3, 241, 1969. **
C ** **
C ** PRINCIPAL VARIABLES: **
C ** **
C ** INTEGER N NUMBER OF MOLECULES **
C ** REAL RX(N),RY(N),RZ(N) POSITIONS **
C ** REAL VOL VOLUME **
C ** REAL BOX BOX LENGTH **
C ** REAL DENS REDUCED DENSITY **
C ** REAL TEMP REDUCED TEMPERATURE **
C ** REAL SIGMA REDUCED LJ DIAMETER **
C ** REAL DRMAX MAXIMUM DISPLACEMENT **
C ** REAL V THE POTENTIAL ENERGY **
C ** REAL W THE VIRIAL **
C ** REAL PRESUR REQUIRED PRESSURE **
C ** REAL DBOXMX MAX CHANGE IN BOX **
C ** **
C ** USAGE: **
C ** **
C ** CONDUCTS MONTE CARLO SIMULATION AT CONSTANT PRESSURE FOR A **
C ** SPECIFIED NUMBER OF CYCLES FROM A GIVEN INITIAL CONFIGURATION.**
C ** **
C ** UNITS: **
C ** **
C ** THIS PROGRAM USES THE USUAL REDUCED LJ UNITS. IN PARTICULAR **
C ** THE BOX LENGTH IS IN UNITS OF SIGMA. **
C ** **
C ** ROUTINES REFERENCED: **
C ** **
C ** SUBROUTINE ENERGY ( RXI, RYI, RZI, I, RCUT, BOX, V12, V6, **
C ** : W12, W6 ) **
C ** CALCULATES THE ENERGY AND VIRIAL FOR ATOM I IN THE FLUID **
C ** REAL FUNCTION RANF ( DUMMY ) **
C ** RETURNS A UNIFORM RANDOM NUMBER BETWEEN ZERO AND ONE **
C ** SUBROUTINE READCN ( CNFILE, BOX ) **
C ** READS IN CONFIGURATION AND BOX VARIABLES **
C ** SUBROUTINE SUMUP ( RCUT, RMIN, OVRLAP, BOX, V12, V6, W12, W6 )**
C ** CALCULATES POTENTIAL AND VIRIAL FOR A CONFIGURATION **
C ** SUBROUTINE WRITCN ( CNFILE, BOX ) **
C ** WRITES OUT CONFIGURATION AND BOX VARIABLES **
C *******************************************************************
INTEGER N
PARAMETER ( N = 108 )
REAL RX(N), RY(N), RZ(N)
INTEGER STEP, NSTEP, IPRINT, IRATIO, IRATB, I
REAL ACV, ACP, ACD, ACM, ACATMA, ACBOXA, NORM
REAL ACVSQ, ACPSQ, ACDSQ
REAL AVV, AVP, AVD
REAL FLV, FLP, FLD
REAL DENS, TEMP, RCUT, RMIN, PRESUR, BOX, VOL, PRES, VN
REAL BOXINV, BOXNEW, RATBOX, RAT12, RAT6, DVOL, DPV
REAL DELTHB, DRMAX, DBOXMX, BETA, RANF, DUMMY, RATIO
REAL RRBOX, BRATIO, DELTVB, RCUTN
REAL RXIOLD, RYIOLD, RZIOLD, RXINEW, RYINEW, RZINEW
REAL V12OLD, V6OLD, V12NEW, V6NEW, VS
REAL W12OLD, W6OLD, W12NEW, W6NEW, WS, PS
REAL DELV12, DELV6, DELW12, DELW6, DELTV
REAL V6, W6, V12, W12, VLRC, VLRCN, WLRC, WLRCN
REAL SR3, SR9, VLRC6, WLRC6, VLRC12, WLRC12, PI
CHARACTER TITLE*80, CNFILE*30
LOGICAL OVRLAP
PARAMETER ( PI = 3.1415927 )
C *******************************************************************
WRITE(*,'(1H1,'' **** PROGRAM MCNPT **** '')')
WRITE(*,'(//'' MONTE CARLO IN A CONSTANT-NPT ENSEMBLE '')')
WRITE(*,'( '' LENNARD-JONES ATOMS '')')
C ** BASIC SIMULATION PARAMETERS **
WRITE(*,'('' ENTER RUN TITLE '')')
READ (*,'(A)') TITLE
WRITE(*,'('' ENTER CONFIGURATION FILENAME '')')
READ (*,'(A)') CNFILE
WRITE(*,'('' ENTER NUMBER OF CYCLES '')')
READ (*,*) NSTEP
WRITE(*,'('' ENTER INTERVAL BETWEEN PRINTS IN CYCLES '')')
READ (*,*) IPRINT
WRITE(*,'('' ENTER INTERVAL FOR UPDATE OF MAXIMUM '')')
WRITE(*,'('' DISPLACEMENT OF ATOMS IN CYCLES '')')
READ (*,*) IRATIO
WRITE(*,'('' ENTER INTERVAL FOR UPDATE OF MAXIMUM '')')
WRITE(*,'('' DISPLACEMENT OF THE BOX IN CYCLES '')')
READ (*,*) IRATB
WRITE(*,'(/'' ENTER THE FOLLOWING IN L-J REDUCED UNITS ''/)')
WRITE(*,'('' ENTER POTENTIAL CUTOFF '')')
READ (*,*) RCUT
WRITE(*,'('' ENTER DESIRED PRESSURE '')')
READ (*,*) PRESUR
WRITE(*,'('' ENTER DESIRED TEMPERATURE '')')
READ (*,*) TEMP
WRITE(*,'(//1X,A)') TITLE
WRITE(*,'('' CONFIGURATION FILENAME '',A)') CNFILE
WRITE(*,'('' NUMBER OF CYCLES = '',I6)') NSTEP
WRITE(*,'('' PRINT INTERVAL = '',I6)') IPRINT
WRITE(*,'('' RATIO UPDATE INTERVAL FOR ATOMS = '',I6)') IRATIO
WRITE(*,'('' RATIO UPDATE INTERVAL FOR BOX = '',I6)') IRATB
WRITE(*,'('' POTENTIAL CUTOFF = '',F10.5)')RCUT
WRITE(*,'('' DESIRED PRES. = '',F10.5)')PRESUR
WRITE(*,'('' DESIRED TEMP = '',F10.5)')TEMP
C ** READCN READS IN INITIAL CONFIGURATION **
C ** ORIGIN SHOULD BE AT CENTRE OF BOX **
CALL READCN ( CNFILE, BOX )
C ** SET DEPENDENT PARAMETERS **
VOL = BOX ** 3
BOXINV = 1.0 / BOX
DENS = REAL ( N ) / VOL
IF ( RCUT .GT. ( 0.5 * BOX ) ) STOP 'CUT-OFF TOO LARGE'
DBOXMX = BOX / 40.0
DRMAX = 0.15
RMIN = 0.70
BETA = 1.0 / TEMP
WRITE(*,'('' INITIAL DENSITY = '',F10.5)') DENS
C ** CALCULATE LONG-RANGE CORRECTIONS FOR LJ POTENTIAL. **
C ** 6 IS FOR ATTRACTIVE CONTRIBUTIONS 12 IS FOR REPULSIVE **
SR3 = ( 1.0 / RCUT ) ** 3
SR9 = SR3 ** 3
VLRC12 = 8.0 * PI * DENS * REAL ( N ) * SR9 / 9.0
VLRC6 = - 8.0 * PI * DENS * REAL ( N ) * SR3 / 3.0
WLRC12 = 4.0 * VLRC12
WLRC6 = 2.0 * VLRC6
VLRC = VLRC12 + VLRC6
WLRC = WLRC12 + WLRC6
C ** ZERO ACCUMULATORS **
ACM = 0.0
ACATMA = 0.0
ACBOXA = 0.0
ACV = 0.0
ACP = 0.0
ACD = 0.0
ACVSQ = 0.0
ACPSQ = 0.0
ACDSQ = 0.0
FLV = 0.0
FLP = 0.0
FLD = 0.0
C ** CALCULATE INITIAL ENERGY AND VIRIAL **
CALL SUMUP ( RCUT, RMIN, OVRLAP, BOX, V12, V6, W12, W6 )
IF ( OVRLAP ) STOP 'OVERLAP IN INITIAL CONFIGURATION'
C ** CALCULATE THE INITIAL ENERGY AND VIRIAL **
VS = ( V12 + V6 + VLRC ) / REAL ( N )
WS = ( W12 + W6 + WLRC ) / REAL ( N )
PS = DENS * TEMP + ( W12 + W6 + WLRC ) / VOL
C ** ADD LONG RANGE CORRECTIONS **
C ** INTO THE ENERGY AND VIRIAL **
V12 = V12 + VLRC12
V6 = V6 + VLRC6
W12 = W12 + WLRC12
W6 = W6 + WLRC6
WRITE(*,'('' INITIAL V/N = '',F10.6)') VS
WRITE(*,'('' INITIAL W/N = '',F10.6)') WS
WRITE(*,'('' INITIAL P = '',F10.6)') PS
WRITE(*,'(//'' **** START OF MARKOV CHAIN ****'')')
WRITE(*,10001)
C *******************************************************************
C ** MAIN LOOP STARTS **
C *******************************************************************
DO 100 STEP = 1, NSTEP
C ** LOOP OVER ATOMS STARTS **
DO 97 I = 1, N
RXIOLD = RX(I)
RYIOLD = RY(I)
RZIOLD = RZ(I)
C ** CALCULATE V FOR AN ATOM IN OLD STATE **
CALL ENERGY ( RXIOLD, RYIOLD, RZIOLD, I, RCUT, BOX,
: V12OLD, V6OLD, W12OLD, W6OLD )
C ** MOVE ATOM I AND PICKUP CENTRAL IMAGE **
RXINEW = RXIOLD + ( 2.0 * RANF ( DUMMY ) - 1.0 ) * DRMAX
RYINEW = RYIOLD + ( 2.0 * RANF ( DUMMY ) - 1.0 ) * DRMAX
RZINEW = RZIOLD + ( 2.0 * RANF ( DUMMY ) - 1.0 ) * DRMAX
RXINEW = RXINEW - ANINT ( RXINEW * BOXINV ) * BOX
RYINEW = RYINEW - ANINT ( RYINEW * BOXINV ) * BOX
RZINEW = RZINEW - ANINT ( RZINEW * BOXINV ) * BOX
C ** CALCULATE V FOR ATOM IN NEW STATE **
CALL ENERGY ( RXINEW, RYINEW, RZINEW, I, RCUT, BOX,
: V12NEW, V6NEW, W12NEW, W6NEW )
C ** CHECK FOR ACCEPTANCE **
DELV12 = V12NEW - V12OLD
DELV6 = V6NEW - V6OLD
DELW12 = W12NEW - W12OLD
DELW6 = W6NEW - W6OLD
DELTV = DELV12 + DELV6
DELTVB = BETA * DELTV
IF ( DELTVB .LT. 75.0 ) THEN
IF ( DELTV .LE. 0.0 ) THEN
V12 = V12 + DELV12
V6 = V6 + DELV6
W12 = W12 + DELW12
W6 = W6 + DELW6
RX(I) = RXINEW
RY(I) = RYINEW
RZ(I) = RZINEW
ACATMA = ACATMA + 1.0
ELSEIF ( EXP ( - DELTVB ) .GT. RANF ( DUMMY ) ) THEN
V12 = V12 + DELV12
V6 = V6 + DELV6
W12 = W12 + DELW12
W6 = W6 + DELW6
RX(I) = RXINEW
RY(I) = RYINEW
RZ(I) = RZINEW
ACATMA = ACATMA + 1.0
ENDIF
ENDIF
VN = ( V12 + V6 ) / REAL ( N )
PRES = DENS * TEMP + ( W12 + W6 ) / VOL
C ** INCREMENT ACCUMULATORS **
ACM = ACM + 1.0
ACV = ACV + VN
ACP = ACP + PRES
ACD = ACD + DENS
ACVSQ = ACVSQ + VN ** 2
ACPSQ = ACPSQ + PRES ** 2
ACDSQ = ACDSQ + DENS ** 2
97 CONTINUE
C ** ENDS LOOP OVER ATOMS IN ONE CYCLE **
C ** ATTEMPT A BOX MOVE **
BOXNEW = BOX + ( 2.0 * RANF ( DUMMY ) - 1.0 ) * DBOXMX
RATBOX = BOX / BOXNEW
RRBOX = 1.0 / RATBOX
RCUTN = RCUT * RRBOX
C ** CALCULATE SCALING PARAMETERS **
RAT6 = RATBOX ** 6
RAT12 = RAT6 * RAT6
C ** SCALE ENERGY, AND VIRIAL INCLUDING LRC **
V12NEW = V12 * RAT12
V6NEW = V6 * RAT6
W12NEW = W12 * RAT12
W6NEW = W6 * RAT6
C ** CALCULATE CHANGE IN ENERGY AND VOLUME **
DELTV = V12NEW + V6NEW - V12 - V6
DPV = PRESUR * ( BOXNEW ** 3 - VOL )
DVOL = 3.0 * TEMP * REAL ( N ) * ALOG ( RATBOX )
DELTHB = BETA * ( DELTV + DPV + DVOL )
C ** CHECK FOR ACCEPTANCE **
IF ( DELTHB .LT. 75.0 ) THEN
IF ( DELTHB .LE. 0.0 ) THEN
V12 = V12NEW
V6 = V6NEW
W12 = W12NEW
W6 = W6NEW
DO 98 I = 1, N
RX(I) = RX(I) * RRBOX
RY(I) = RY(I) * RRBOX
RZ(I) = RZ(I) * RRBOX
98 CONTINUE
BOX = BOXNEW
RCUT = RCUTN
ACBOXA = ACBOXA + 1.0
ELSEIF ( EXP ( - DELTHB ) .GT. RANF ( DUMMY ) ) THEN
V12 = V12NEW
V6 = V6NEW
W12 = W12NEW
W6 = W6NEW
DO 99 I = 1, N
RX(I) = RX(I) * RRBOX
RY(I) = RY(I) * RRBOX
RZ(I) = RZ(I) * RRBOX
99 CONTINUE
BOX = BOXNEW
RCUTN = RCUT
ACBOXA = ACBOXA + 1.0
ENDIF
ENDIF
BOXINV = 1.0 / BOX
VOL = BOX ** 3
DENS = REAL ( N ) / VOL
C ** CALCULATE ENERGY AND PRESSURE **
VN = ( V12 + V6 ) / REAL ( N )
PRES = DENS * TEMP + ( W12 + W6 ) / VOL
C ** INCREMENT ACCUMULATORS **
ACM = ACM + 1.0
ACV = ACV + VN
ACP = ACP + PRES
ACD = ACD + DENS
ACVSQ = ACVSQ + VN ** 2
ACPSQ = ACPSQ + PRES ** 2
ACDSQ = ACDSQ + DENS ** 2
C ** ENDS ATTEMPTED BOX MOVE **
C ** PERFORM PERIODIC OPERATIONS **
IF ( MOD ( STEP, IRATIO ) .EQ. 0 ) THEN
C ** ADJUST MAXIMUM DISPLACEMENT FOR ATOMS **
RATIO = ACATMA / REAL ( N * IRATIO )
IF ( RATIO .GT. 0.5 ) THEN
DRMAX = DRMAX * 1.05
ELSE
DRMAX = DRMAX * 0.95
ENDIF
ACATMA = 0.0
ENDIF
IF ( MOD ( STEP, IRATB ) .EQ. 0 ) THEN
C ** ADJUST MAXIMUM DISPLACEMENT FOR THE BOX **
BRATIO = ACBOXA / REAL ( IRATB )
IF ( BRATIO .GT. 0.5 ) THEN
DBOXMX = DBOXMX * 1.05
ELSE
DBOXMX = DBOXMX * 0.95
ENDIF
ACBOXA = 0.0
ENDIF
IF ( MOD ( STEP, IPRINT ) .EQ. 0 ) THEN
C ** OPTIONALLY PRINT INFORMATION **
WRITE(*,'(1X,I8,5(2X,F10.5))')
: STEP, VN, PRES, DENS, RATIO, BRATIO
ENDIF
100 CONTINUE
C *******************************************************************
C ** MAIN LOOP ENDS **
C *******************************************************************
WRITE(*,'(/1X,''**** END OF MARKOV CHAIN **** ''//)')
C ** WRITE OUT FINAL CONFIGURATION AND BOXLENGTH **
CALL WRITCN ( CNFILE, BOX )
C ** WRITE OUT FINAL AVERAGES **
NORM = REAL ( ACM )
AVV = ACV / NORM
AVP = ACP / NORM
AVD = ACD / NORM
ACVSQ = ( ACVSQ / NORM ) - AVV ** 2
ACPSQ = ( ACPSQ / NORM ) - AVP ** 2
ACDSQ = ( ACDSQ / NORM ) - AVD ** 2
IF ( ACVSQ .GT. 0.0 ) FLV = SQRT ( ACVSQ )
IF ( ACPSQ .GT. 0.0 ) FLP = SQRT ( ACPSQ )
IF ( ACDSQ .GT. 0.0 ) FLD = SQRT ( ACDSQ )
WRITE(*, 10002)
WRITE(*,'('' AVERAGES'',3(2X,F10.5))') AVV, AVP, AVD
WRITE(*,'('' FLUCTS '',3(2X,F10.5))') FLV, FLP, FLD
STOP
10001 FORMAT(//1X,' CYCLE ..POTENT.. ..PRESSURE.. ..DENSITY..',
: ' ..RATIO.. ..BRATIO..')
10002 FORMAT(//1X,' CYCLE ..POTENT.. ..PRESSURE.. ..DENSITY..')
END
SUBROUTINE SUMUP ( RCUT, RMIN, OVRLAP, BOX, V12, V6, W12, W6 )
COMMON / BLOCK1 / RX, RY, RZ
C *******************************************************************
C ** CALCULATES THE TOTAL POTENTIAL ENERGY FOR A CONFIGURATION. **
C ** **
C ** PRINCIPAL VARIABLES: **
C ** **
C ** INTEGER N THE NUMBER OF ATOMS **
C ** REAL RX(N(,RY(N),RZ(N) THE POSITIONS OF THE ATOMS **
C ** REAL V THE POTENTIAL ENERGY **
C ** REAL W THE VIRIAL **
C ** LOGICAL OVRLAP TRUE FOR SUBSTANTIAL ATOM OVERLAP **
C ** **
C ** USAGE: **
C ** **
C ** THE SUBROUTINE RETURNS THE TOTAL POTENTIAL ENERGY AT THE **
C ** BEGINNING AND END OF THE RUN. **
C *******************************************************************
INTEGER N
PARAMETER ( N = 108 )
REAL RX(N), RY(N), RZ(N)
REAL RMIN, RCUT, V12, V6, W12, W6, BOX
LOGICAL OVRLAP
REAL RCUTSQ, RMINSQ, VIJ12, VIJ6
REAL RXI, RYI, RZI, RXIJ, RYIJ, RZIJ
REAL SR2, SR6, RIJSQ, BOXINV
INTEGER I, J
C *******************************************************************
OVRLAP = .FALSE.
RCUTSQ = RCUT * RCUT
RMINSQ = RMIN * RMIN
BOXINV = 1.0 / BOX
V12 = 0.0
V6 = 0.0
W12 = 0.0
W6 = 0.0
C ** LOOP OVER ALL THE PAIRS IN THE LIQUID **
DO 100 I = 1, N - 1
RXI = RX(I)
RYI = RY(I)
RZI = RZ(I)
DO 99 J = I + 1, N
RXIJ = RXI - RX(J)
RYIJ = RYI - RY(J)
RZIJ = RZI - RZ(J)
RXIJ = RXIJ - ANINT ( RXIJ * BOXINV ) * BOX
RYIJ = RYIJ - ANINT ( RYIJ * BOXINV ) * BOX
RZIJ = RZIJ - ANINT ( RZIJ * BOXINV ) * BOX
RIJSQ = RXIJ * RXIJ + RYIJ * RYIJ + RZIJ * RZIJ
IF ( RIJSQ .LT. RMINSQ ) THEN
OVRLAP = .TRUE.
RETURN
ELSEIF ( RIJSQ .LT. RCUTSQ ) THEN
SR2 = 1.0 / RIJSQ
SR6 = SR2 * SR2 * SR2
VIJ12 = SR6 * SR6
VIJ6 = - SR6
V12 = V12 + VIJ12
V6 = V6 + VIJ6
W12 = W12 + VIJ12
W6 = W6 + VIJ6 * 0.5
ENDIF
99 CONTINUE
100 CONTINUE
V12 = 4.0 * V12
V6 = 4.0 * V6
W12 = 48.0 * W12 / 3.0
W6 = 48.0 * W6 / 3.0
RETURN
END
SUBROUTINE ENERGY ( RXI, RYI, RZI, I, RCUT, BOX,
: V12, V6, W12, W6 )
COMMON / BLOCK1 / RX, RY, RZ
C *******************************************************************
C ** CALCULATES THE POTENTIAL ENERGY OF I WITH ALL OTHER ATOMS **
C ** **
C ** PRINCIPAL VARIABLES: **
C ** **
C ** INTEGER I THE ATOM OF INTEREST **
C ** INTEGER N THE NUMBER OF ATOMS **
C ** REAL RX(N),RY(N),RZ(N) THE ATOM POSITIONS **
C ** REAL RXI,RYI,RZI THE COORDINATES OF ATOM I **
C ** REAL V THE POTENTIAL ENERGY OF ATOM I **
C ** REAL W THE VIRIAL OF ATOM I **
C ** **
C ** USAGE: **
C ** **
C ** THIS SUBROUTINE IS USED TO CALCULATE THE CHANGE OF ENERGY **
C ** DURING A TRIAL MOVE OF ATOM I. IT IS CALLED BEFORE AND **
C ** AFTER THE RANDOM DISPLACEMENT OF I. **
C *******************************************************************
INTEGER N
PARAMETER ( N = 108 )
REAL RX(N), RY(N), RZ(N)
REAL RCUT, BOX, RXI, RYI, RZI, V12, V6, W12, W6
INTEGER I
REAL RCUTSQ, VIJ12, VIJ6
REAL RXIJ, RYIJ, RZIJ, RIJSQ, SR2, SR6, BOXINV
INTEGER J
C ******************************************************************
RCUTSQ = RCUT * RCUT
BOXINV = 1.0 / BOX
V12 = 0.0
V6 = 0.0
W12 = 0.0
W6 = 0.0
C ** LOOP OVER ALL MOLECULES EXCEPT I **
DO 100 J = 1, N
IF ( I .NE. J ) THEN
RXIJ = RXI - RX(J)
RYIJ = RYI - RY(J)
RZIJ = RZI - RZ(J)
RXIJ = RXIJ - ANINT ( RXIJ * BOXINV ) * BOX
RYIJ = RYIJ - ANINT ( RYIJ * BOXINV ) * BOX
RZIJ = RZIJ - ANINT ( RZIJ * BOXINV ) * BOX
RIJSQ = RXIJ * RXIJ + RYIJ * RYIJ + RZIJ * RZIJ
IF ( RIJSQ .LT. RCUTSQ ) THEN
SR2 = 1.0 / RIJSQ
SR6 = SR2 * SR2 * SR2
VIJ12 = SR6 * SR6
VIJ6 = - SR6
V12 = V12 + VIJ12
V6 = V6 + VIJ6
W12 = W12 + VIJ12
W6 = W6 + VIJ6 * 0.5
ENDIF
ENDIF
100 CONTINUE
V12 = 4.0 * V12
V6 = 4.0 * V6
W12 = 48.0 * W12 / 3.0
W6 = 48.0 * W6 / 3.0
RETURN
END
SUBROUTINE READCN ( CNFILE, BOX )
COMMON / BLOCK1 / RX, RY, RZ
C *******************************************************************
C ** SUBROUTINE TO READ IN THE CONFIGURATION FROM UNIT 10 **
C *******************************************************************
INTEGER N
PARAMETER ( N = 108 )
CHARACTER CNFILE*(*)
REAL RX(N), RY(N), RZ(N)
REAL BOX
INTEGER CNUNIT
PARAMETER ( CNUNIT = 10 )
INTEGER NN
C ********************************************************************
OPEN ( UNIT = CNUNIT, FILE = CNFILE, STATUS = 'OLD',
: FORM = 'UNFORMATTED' )
READ ( CNUNIT ) NN, BOX
IF ( NN .NE. N ) STOP 'N ERROR IN READCN'
READ ( CNUNIT ) RX, RY, RZ
CLOSE ( UNIT = CNUNIT )
RETURN
END
SUBROUTINE WRITCN ( CNFILE, BOX )
COMMON / BLOCK1 / RX, RY, RZ
C *******************************************************************
C ** SUBROUTINE TO WRITE OUT THE CONFIGURATION TO UNIT 10 **
C *******************************************************************
INTEGER N
PARAMETER ( N = 108 )
CHARACTER CNFILE*(*)
REAL RX(N), RY(N), RZ(N)
REAL BOX
INTEGER CNUNIT
PARAMETER ( CNUNIT = 10 )
C ********************************************************************
OPEN ( UNIT = CNUNIT, FILE = CNFILE, STATUS = 'UNKNOWN',
: FORM = 'UNFORMATTED' )
WRITE ( CNUNIT ) N, BOX
WRITE ( CNUNIT ) RX, RY, RZ
CLOSE ( UNIT = CNUNIT )
RETURN
END
REAL FUNCTION RANF ( DUMMY )
C *******************************************************************
C ** RETURNS A UNIFORM RANDOM VARIATE IN THE RANGE 0 TO 1. **
C ** **
C ** *************** **
C ** ** WARNING ** **
C ** *************** **
C ** **
C ** GOOD RANDOM NUMBER GENERATORS ARE MACHINE SPECIFIC. **
C ** PLEASE USE THE ONE RECOMMENDED FOR YOUR MACHINE. **
C *******************************************************************
INTEGER L, C, M
PARAMETER ( L = 1029, C = 221591, M = 1048576 )
INTEGER SEED
REAL DUMMY
SAVE SEED
DATA SEED / 0 /
C *******************************************************************
SEED = MOD ( SEED * L + C, M )
RANF = REAL ( SEED ) / M
RETURN
END
|
