PROGRAM METROPOL C---------------------------------------------------------------------- C C USE DOUBLE-LENGTH REAL NUMBERS AND LONG INTEGERS IN THIS PROGRAM C C---------------------------------------------------------------------- C THIS IS A CELL-STRUCTURE MONTE CARLO PROGRAM IN WHICH THE CELL C CONTAINS 125 MOLECULES. THE MOLECULES ARE FIRST PLACED ON A C 5 X 5 X 5 REGULAR GRID AND THEN DISPLACED RANDOMLY BY UP TO D/20 C IN EACH PRINCIPAL DIRECTION WHERE D IS THE INITIAL SEPARATION. C AFTER 20000 CONFIGURATIONS, TO ALLOW THE SYSTEM TO SETTLE DOWN, C EVERY 50TH MOLECULAR CONFIGURATION IS SAMPLED TO FIND THE QUANTITY C SUM(F(I,J)r(I,J))/3 AND ALSO THE RADIAL DENSITY FUNCTION FOR THE C NEXT 80000 CONFIGURATIONS. C THE PROGRAM TRANSFORMS EVERY QUANTITY INTO S.I. UNITS. SOME FORM OF C REDUCED UNITS WOULD KEEP NUMBERS IN A MORE CONVENIENT RANGE BUT MAKE C THE PROGRAM MORE DIFFICULT TO FOLLOW. C C THE ALGORITHM IS BASED ON THAT IN N. Metropolis, A.W.Rosenbluth, C M.N. Rosenbluth and A.H. Teller, J. Chem. Phys. 21, 1087. C DIMENSION X(125,3),POT(125,125),TEMPOT(125),RR(4), +SEPAR(125,4),TAB(0:100),DD(0:100) C DATA FOR THE RANDOM-NUMBER GENERATOR. IR IS THE SEED DATA IR,IX,IY,IM/199,171,11213,53125/ C DATA FOR THE LENNARD-JONES FORCE. SIG IS IN nm AND C EPS AS EPS/k WHERE k IS BOLTZMANN'S CONSTANT. THE FIGURES C ARE FOR ARGON. DATA SIG,EPS/0.345,120/ C CONVERT TO S.I. SIG=SIG*1.0E-9 EPS=EPS*1.38E-23 C T IS THE TEMPERATURE OF THE LIQUID. DATA T/329/ C THE RANGE OF THE FORCE, DLIM, IS SET AT 2.4 X SIG CORRESPONDING TO C WHERE THE FORCE FALLS TO 1% OF ITS MAXIMUM VALUE. THE SIDE OF THE C CUBICAL CELL IS 5 X SIG FOR V* = 1. THE VALUE OF V* IS READ IN C WHICH GIVES THE CELL EDGE, EL. SIG2=SIG*SIG TFEPS=24.0*EPS C BOLTZMANN CONSTANT CAY=1.38E-23 C NUMBER OF GENERATED CONFIGURATIONS NTOT=100000 PI=4.0*ATAN(1.0) DLIM2=(2.4*SIG)**2 WRITE(6,'('' READ IN THE VALUE OF V* '')') READ(5,*)VSTAR EL=5.0*SIG*VSTAR**(1.0/3.0) GRID=0.2*EL VIR=0 C THE MOLECULES ARE INITIALLY SET ON A UNIFORM GRID. C CLEAR TABLE FOR RADIAL DISTRIBUTION DO 17 I=1,100 TAB(I)=0 17 CONTINUE NUM=0 DO 1 I=1,5 DO 1 J=1,5 DO 1 K=1,5 NUM=NUM+1 X(NUM,1)=(I-3)*GRID X(NUM,2)=(J-3)*GRID X(NUM,3)=(K-3)*GRID 1 CONTINUE C CLEAR POTENTIAL ENERGY TABLE DO 15 II=1,125 DO 15 JJ=1,125 POT(II,JJ)=0 15 CONTINUE C CALCULATE INITIAL POTENTIAL ENERGY TABLE DO 16 II=1,124 DO 16 JJ=II+1,125 CALL DIS(X,SEPAR,EL,II,JJ) IF(SEPAR(JJ,4).GT.DLIM2)GOTO 16 A=(SIG2/SEPAR(JJ,4))**3 POT(II,JJ)=4*EPS*A*(A-1.0) POT(JJ,II)=POT(II,JJ) 16 CONTINUE NOSTEP=0 100 NOSTEP=NOSTEP+1 IF((NOSTEP/50)*50.EQ.NOSTEP)WRITE(6,125)NOSTEP 125 FORMAT(15H STARTING STEP ,I5) C CHOOSE A MOLECULE IR=MOD(IR*IX+IY,IM) MOL=INT(NUM*FLOAT(IR)/FLOAT(IM))+1 C CHOOSE A DIRECTION 98 DO 19 K=1,3 IR=MOD(IR*IX+IY,IM) RR(K)=2*FLOAT(IR)/FLOAT(IM)-1.0 19 CONTINUE RR(4)=RR(1)**2+RR(2)**2+RR(3)**2 IF(RR(4).GT.1.0)GOTO 98 RS=SQRT(RR(4)) C CHOOSE A DISTANCE FOR THE MOLECULE TO BE MOVED IR=MOD(IR*IX+IY,IM) DEL=0.25*GRID*FLOAT(IR)/FLOAT(IM) C CLEAR TEMPORARY POTENTIAL TABLE DO 27 J=1,125 TEMPOT(J)=0 27 CONTINUE C CALCULATE NEW POTENTIAL CONTRIBUTIONS DO 21 J=1,125 IF(J.EQ.MOL)GOTO 21 DO 22 K=1,3 SEPAR(J,K)=X(MOL,K)+DEL*RR(K)/RS-X(J,K) IF(ABS(SEPAR(J,K)+EL).LT.ABS(SEPAR(J,K)))SEPAR(J,K)=SEPAR(J,K)+EL IF(ABS(SEPAR(J,K)-EL).LT.ABS(SEPAR(J,K)))SEPAR(J,K)=SEPAR(J,K)-EL 22 CONTINUE SEPAR(J,4)=SEPAR(J,1)**2+SEPAR(J,2)**2+SEPAR(J,3)**2 IF(SEPAR(J,4).GT.DLIM2)GOTO 21 A=(SIG2/SEPAR(J,4))**3 TEMPOT(J)=4*EPS*A*(A-1.0) 21 CONTINUE C FIND DELTAPHI SUM1=0 SUM2=0 DO 28 J=1,125 SUM1=SUM1+POT(MOL,J) SUM2=SUM2+TEMPOT(J) 28 CONTINUE DELPHI=SUM2-SUM1 C TEST DELPHI IF(DELPHI.LT.0)THEN CALL CHANGE(POT,TEMPOT,X,RR,DEL,EL,MOL) GOTO 55 ENDIF C DELPHI IS POSITIVE. DECIDE ON NEXT CONFIGURATION Z=EXP(-DELPHI/CAY/T) IR=MOD(IR*IX+IY,IM) R=FLOAT(IR)/FLOAT(IM) IF(Z.GT.R)THEN CALL CHANGE(POT,TEMPOT,X,RR,DEL,EL,MOL) ENDIF 55 IF(NOSTEP.LE.50000)GOTO 100 IF((NOSTEP/10)*10.NE.NOSTEP)GOTO 100 C CALCULATE CONTRIBUTIONS TO THE VIRIAL TERM AND TO THE RADIAL C DISTRIBUTION FUNCTION. DO 40 I=1,124 DO 41 J=I+1,125 CALL DIS(X,SEPAR,EL,I,J) IF(SEPAR(J,4).GT.DLIM2)GOTO 47 A=(SIG2/SEPAR(J,4))**3 C=TFEPS*A*(1.0-2.0*A) C ADD CONTRIBUTION TO VIRIAL TERM VIR=VIR-C/3.0 C ADD RADIAL DISTANCE TO TABLE 47 L=INT(20.0*SQRT(ABS(SEPAR(J,4)))/SIG) IF(L.GT.48)GOTO 41 TAB(L)=TAB(L)+1 41 CONTINUE 40 CONTINUE IF(NOSTEP.LT.NTOT)GOTO 100 C TAKE THE AVERAGE OF THE VIRIAL TERM OVER 1600 TIMESTEPS VIR=VIR/5000.0 C NOW CALCULATE THE FACTOR 1 + VIR/NkT FACTOR=1.0+VIR/CAY/T/FLOAT(NUM) OPEN(UNIT=9,FILE='LPT1') C OUTPUT TEMPERATURE AND VSTAR WRITE(9,350)T,VSTAR 350 FORMAT(15H TEMPERATURE = ,F6.1,9H VSTAR = ,F7.2) WRITE(9,'('' '')') C OUTPUT FACTOR WRITE(9,200)FACTOR 200 FORMAT(21H THE FACTOR PV/NkT = ,F7.3) C OUTPUT THE RADIAL DENSITY DISTRIBUTION 4*PI*R**2*RO(R) ON AN C ABSOLUTE SCALE RELATIVE TO THE AVERAGE DENSITY. DO 90 I=0,47 TAB(I)=6.161E-3*VSTAR*TAB(I)/((I+1.0)**3-I**3) DD(I)=0.05*(I+0.5) 90 CONTINUE WRITE(9,'('' '')') WRITE(9,'('' RADIAL DENSITY FUNCTION ON AN ABOLUTE SCALE '')') WRITE(9,300)(DD(I),TAB(I),I=0,47) 300 FORMAT(2(F10.3,F12.6)) STOP END SUBROUTINE DIS(X,S,EL,I,J) DIMENSION X(125,3),S(125,4) DO 20 K=1,3 S(J,K)=X(I,K)-X(J,K) IF(ABS(S(J,K)+EL).LT.ABS(S(J,K)))S(J,K)=S(J,K)+EL IF(ABS(S(J,K)-EL).LT.ABS(S(J,K)))S(J,K)=S(J,K)-EL 20 CONTINUE S(J,4)=S(J,1)**2+S(J,2)**2+S(J,3)**2 RETURN END SUBROUTINE CHANGE(P,T,X,R,D,E,M) DIMENSION P(125,125),T(125),X(125,3),R(4) RS=SQRT(R(4)) DO 2 J=1,125 P(M,J)=T(J) P(J,M)=T(J) 2 CONTINUE C MOVE MOLECULE M AND PLACE IN CELL DO 3 K=1,3 X(M,K)=X(M,K)+D*R(K)/RS X(M,K)=AMOD(X(M,K)+9.5*E,E)-0.5*E 3 CONTINUE RETURN END