      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
