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