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                                                                     
