C    *******************************************************************        
C    ** FICHE F.31                                                    **        
C    ** CONSTANT-NPT MOLECULAR DYNAMICS USING CONSTRAINT ALGORITHM.   **        
C    *******************************************************************        
                                                                                
        PROGRAM EVAMOR                                                          
                                                                                
        COMMON / BLOCK1 / RX, RY, RZ, RX1, RY1, RZ1,                            
     :                    RX2, RY2, RZ2, RX3, RY3, RZ3,                         
     :                    PX, PY, PZ, PX1, PY1, PZ1,                            
     :                    PX2, PY2, PZ2, PX3, PY3, PZ3,                         
     :                    FX, FY, FZ                                            
        COMMON / BLOCK2 / BOX, BOX1, BOX2, BOX3                                 
                                                                                
C    *******************************************************************        
C    ** CONSTANT-NPT MOLECULAR DYNAMICS USING CONSTRAINT ALGORITHM.   **        
C    **                                                               **        
C    ** THE MODIFIED EQUATIONS OF MOTION ARE AS FOLLOWS:              **        
C    **       R1 = P/M + CHI*R                                        **        
C    **       P1 = F - CHI*P - XI*P                                   **        
C    **       V1 = 3*V*CHI                                            **        
C    ** WHERE R1 IS THE TIME DERIVATIVE OF POSITION R                 **        
C    ** P1 IS THE TIME DERIVATIVE OF MOMENTUM P                       **        
C    ** V1 IS THE TIME DERIVATIVE OF VOLUME V                         **        
C    ** AND CHI AND XI ARE LAGRANGE MULTIPLIERS                       **        
C    **       CHI = ( - SUM( (R.P) X / (R.R) ) ) / ( SUM(X) + 9PV )   **        
C    **       XI  = SUM(P.F) / SUM(P.P) - CHI                         **        
C    ** HERE TERMS LIKE (P.F) ARE SCALAR PRODUCTS OVER PAIR TERMS AND **        
C    ** PV STANDS FOR PRESSURE * VOLUME.                              **        
C    ** WE SOLVE THESE EQUATIONS BY A GEAR 4-VALUE METHOD FOR FIRST   **        
C    ** ORDER DIFFERENTIAL EQUATIONS                                  **        
C    **                                                               **        
C    ** REFERENCES:                                                   **        
C    **                                                               **        
C    ** EVANS AND MORRISS, CHEM PHYS 77, 63, 1983.                    **        
C    ** EVANS AND MORRISS, COMPUT PHYS REP 1, 297, 1984.              **        
C    **                                                               **        
C    ** PRINCIPAL VARIABLES:                                          **        
C    **                                                               **        
C    ** INTEGER N                             NUMBER OF MOLECULES     **        
C    ** REAL    DT                            TIMESTEP                **        
C    ** REAL    RX(N),RY(N),RZ(N)             POSITIONS               **        
C    ** REAL    RX1(N),RY1(N),RZ1(N)          FIRST DERIVATIVES       **        
C    ** REAL    RX2(N),RY2(N),RZ2(N)          SECOND DERIVATIVES      **        
C    ** REAL    RX3(N),RY3(N),RZ3(N)          THIRD DERIVATIVES       **        
C    ** REAL    PX(N),PY(N),PZ(N)             MOMENTA                 **        
C    ** REAL    PX1(N),PY1(N),PZ1(N)          FIRST DERIVATIVES       **        
C    ** REAL    PX2(N),PY2(N),PZ2(N)          SECOND DERIVATIVES      **        
C    ** REAL    PX3(N),PY3(N),PZ3(N)          THIRD DERIVATIVES       **        
C    ** REAL    BOX,BOX1,BOX2,BOX3            BOX LENGTH AND DERIVS   **        
C    ** REAL    FX(N),FY(N),FZ(N)             TOTAL FORCES            **        
C    **                                                               **        
C    ** ROUTINES REFERENCED                                           **        
C    **                                                               **        
C    ** SUBROUTINE READCN ( CNFILE )                                  **        
C    **    READS IN CONFIGURATION AND BOX VARIABLES                   **        
C    ** SUBROUTINE FORCE ( RCUT, V, W, X, RPX, PF )                   **        
C    **    CALCULATES FORCES, POTENTIAL, VIRIAL, HYPERVIRIAL ETC.     **        
C    ** SUBROUTINE KINET ( K )                                        **        
C    **    CALCULATES KINETIC ENERGY                                  **        
C    ** SUBROUTINE PREDIC ( DT )                                      **        
C    **    PREDICTOR ROUTINE FOR CONFIGURATION AND BOX VARIABLES      **        
C    ** SUBROUTINE CORREC ( DT, CHI, CHIPXI )                         **        
C    **    CORRECTOR ROUTINE FOR CONFIGURATION AND BOX VARIABLES      **        
C    ** SUBROUTINE WRITCN ( CNFILE )                                  **        
C    **    WRITES OUT CONFIGURATION AND BOX VARIABLES                 **        
C    ** SUBROUTINE SCALE ( TEMPER, PRESUR, WLRC0, XLRC0, RCUT )       **        
C    **    SCALES MOMENTA AND BOX SIZE TO GIVE DESIRED T AND P        **        
C    *******************************************************************        
                                                                                
        INTEGER     N                                                           
        PARAMETER ( N = 108 )                                                   
                                                                                
        REAL        RX(N), RY(N), RZ(N)                                         
        REAL        RX1(N), RY1(N), RZ1(N)                                      
        REAL        RX2(N), RY2(N), RZ2(N)                                      
        REAL        RX3(N), RY3(N), RZ3(N)                                      
        REAL        PX(N), PY(N), PZ(N)                                         
        REAL        PX1(N), PY1(N), PZ1(N)                                      
        REAL        PX2(N), PY2(N), PZ2(N)                                      
        REAL        PX3(N), PY3(N), PZ3(N)                                      
        REAL        FX(N), FY(N), FZ(N)                                         
        REAL        BOX, BOX1, BOX2, BOX3                                       
                                                                                
        INTEGER     STEP, NSTEP, IPRINT, ISCALE                                 
        REAL        ACV, ACK, ACE, ACP, ACT, ACH, ACD                           
        REAL        ACVSQ, ACKSQ, ACESQ, ACPSQ, ACTSQ, ACHSQ, ACDSQ             
        REAL        AVV, AVK, AVE, AVP, AVT, AVH, AVD                           
        REAL        FLV, FLK, FLE, FLP, FLT, FLH, FLD                           
        REAL        DT, DENS, TEMP, RCUT, VOL, PRES, H, NORM                    
        REAL        TEMPER, PRESUR                                              
        REAL        K, V, W, E                                                  
        REAL        KN, VN, EN, HN                                              
        REAL        X, RPX, PF, PP, PV, CHI, CHIPXI                             
        REAL        SR3, SR9, VLRC, WLRC, XLRC, VLRC0, WLRC0, XLRC0, PI         
        CHARACTER   TITLE*80, CNFILE*30                                         
        REAL        FREE                                                        
                                                                                
        PARAMETER ( FREE = 3.0 )                                                
        PARAMETER ( PI = 3.1415927 )                                            
                                                                                
C    *******************************************************************        
                                                                                
        WRITE(*,'(1H1,'' **** PROGRAM EVAMOR ****                  '')')        
        WRITE(*,'(//1X,''MOLECULAR DYNAMICS OF LENNARD-JONES ATOMS '')')        
        WRITE(*,'(1X,''CONSTANT-NPT ALGORITHM OF EVANS AND MORRISS '')')        
                                                                                
C    ** BASIC SIMULATION PARAMETERS **                                          
                                                                                
        WRITE(*,'('' ENTER RUN TITLE                               '')')        
        READ (*,'(A)') TITLE                                                    
        WRITE(*,'('' ENTER CONFIGURATION FILENAME                  '')')        
        READ (*,'(A)') CNFILE                                                   
        WRITE(*,'('' ENTER NUMBER OF STEPS                         '')')        
        READ (*,*) NSTEP                                                        
        WRITE(*,'('' ENTER INTERVAL BETWEEN PRINTS                 '')')        
        READ (*,*) IPRINT                                                       
        WRITE(*,'('' ENTER INTERVAL BETWEEN BOX AND MOMENTUM SCALES'')')        
        READ (*,*) ISCALE                                                       
        WRITE(*,'('' ENTER THE FOLLOWING IN L-J REDUCED UNITS      '')')        
        WRITE(*,'('' ENTER TIMESTEP                                '')')        
        READ (*,*) DT                                                           
        WRITE(*,'('' ENTER POTENTIAL CUTOFF                        '')')        
        READ (*,*) RCUT                                                         
        WRITE(*,'('' ENTER DESIRED TEMPERATURE                     '')')        
        READ (*,*) TEMPER                                                       
        WRITE(*,'('' ENTER DESIRED PRESSURE                        '')')        
        READ (*,*) PRESUR                                                       
                                                                                
        WRITE(*,'(//1X,A)') TITLE                                               
        WRITE(*,'('' CONFIGURATION FILENAME '',A)') CNFILE                      
        WRITE(*,'('' NUMBER OF STEPS  = '',I6   )') NSTEP                       
        WRITE(*,'('' PRINT INTERVAL   = '',I6   )') IPRINT                      
        WRITE(*,'('' SCALE INTERVAL   = '',I6   )') ISCALE                      
        WRITE(*,'('' TIMESTEP         = '',F10.5)') DT                          
        WRITE(*,'('' POTENTIAL CUTOFF = '',F10.5)') RCUT                        
        WRITE(*,'('' DESIRED TEMP.    = '',F10.5)') TEMPER                      
        WRITE(*,'('' DESIRED PRES.    = '',F10.5)') PRESUR                      
                                                                                
C    ** READCN MUST READ IN INITIAL CONFIGURATION    **                         
C    ** AND ASSIGN VALUES TO BOX AND ITS DERIVATIVES **                         
                                                                                
        CALL READCN ( CNFILE )                                                  
        VOL = BOX ** 3                                                          
        DENS = REAL ( N ) / VOL                                                 
                                                                                
        WRITE(*,'('' INITIAL DENS.    = '',F10.5)') DENS                        
                                                                                
        IF ( IPRINT .LE. 0 ) IPRINT = NSTEP + 1                                 
        IF ( ISCALE .LE. 0 ) ISCALE = NSTEP + 1                                 
                                                                                
C    ** PREPARE FACTORS FOR LONG-RANGE CORRECTIONS **                           
C    ** NB: SPECIFIC TO LENNARD-JONES POTENTIAL    **                           
                                                                                
        SR3 = ( 1.0 / RCUT ) ** 3                                               
        SR9 = SR3 ** 3                                                          
        VLRC0 =   8.0 * PI * REAL ( N ) * SR9 / 9.0                             
     :          - 8.0 * PI * REAL ( N ) * SR3 / 3.0                             
        WLRC0 =  32.0 * PI * REAL ( N ) * SR9 / 9.0                             
     :         - 16.0 * PI * REAL ( N ) * SR3 / 3.0                             
        XLRC0 = 128.0 * PI * REAL ( N ) * SR9 / 9.0                             
     :         - 32.0 * PI * REAL ( N ) * SR3 / 3.0                             
                                                                                
C    ** ZERO ACCUMULATORS **                                                    
                                                                                
        ACV = 0.0                                                               
        ACK = 0.0                                                               
        ACE = 0.0                                                               
        ACP = 0.0                                                               
        ACT = 0.0                                                               
        ACH = 0.0                                                               
        ACD = 0.0                                                               
                                                                                
        ACVSQ = 0.0                                                             
        ACKSQ = 0.0                                                             
        ACESQ = 0.0                                                             
        ACPSQ = 0.0                                                             
        ACTSQ = 0.0                                                             
        ACHSQ = 0.0                                                             
        ACDSQ = 0.0                                                             
                                                                                
        FLV = 0.0                                                               
        FLK = 0.0                                                               
        FLE = 0.0                                                               
        FLP = 0.0                                                               
        FLT = 0.0                                                               
        FLH = 0.0                                                               
        FLD = 0.0                                                               
                                                                                
        WRITE(*,'(//1X,''**** START OF DYNAMICS ****'')')                       
        WRITE(*,10001)                                                          
                                                                                
C    *******************************************************************        
C    ** MAIN LOOP STARTS                                              **        
C    *******************************************************************        
                                                                                
        DO 1000 STEP = 1, NSTEP                                                 
                                                                                
C       ** IMPLEMENT ALGORITHM **                                               
                                                                                
           CALL PREDIC ( DT )                                                   
                                                                                
           CALL FORCE ( RCUT, V, W, X, RPX, PF )                                
           CALL KINET ( K )                                                     
                                                                                
C       ** INCLUDE LONG-RANGE CORRECTIONS IN ALGORITHM **                       
                                                                                
           KN   = K / REAL ( N )                                                
           VOL  = BOX ** 3                                                      
           DENS = REAL ( N ) / VOL                                              
           WLRC = WLRC0 * DENS                                                  
           XLRC = XLRC0 * DENS                                                  
           TEMP = 2.0 * KN / FREE                                               
           PRES = DENS * TEMP + ( W + WLRC ) / VOL                              
           PV   = PRES * VOL                                                    
           PP   = 2.0 * K                                                       
                                                                                
           CHI    = - RPX / 9.0 / ( X + XLRC + PV )                             
           CHIPXI = PF / PP                                                     
                                                                                
           CALL CORREC ( DT, CHI, CHIPXI )                                      
                                                                                
C       ** CALCULATE INSTANTANEOUS VALUES   **                                  
C       ** INCLUDING LONG-RANGE CORRECTIONS **                                  
                                                                                
           VOL  = BOX ** 3                                                      
           DENS = REAL ( N ) / VOL                                              
           VLRC = VLRC0 * DENS                                                  
           WLRC = WLRC0 * DENS                                                  
           XLRC = XLRC0 * DENS                                                  
           V    = V + VLRC                                                      
           W    = W + WLRC                                                      
           X    = X + XLRC                                                      
           E    = K + V                                                         
           VN   = V / REAL ( N )                                                
           EN   = E / REAL ( N )                                                
           TEMP = 2.0 * KN / FREE                                               
           PRES = DENS * TEMP + W / VOL                                         
           H    = E + PRES * VOL                                                
           HN   = H / REAL ( N )                                                
                                                                                
C       ** INCREMENT ACCUMULATORS **                                            
                                                                                
           ACE = ACE + EN                                                       
           ACK = ACK + KN                                                       
           ACV = ACV + VN                                                       
           ACP = ACP + PRES                                                     
           ACH = ACH + HN                                                       
           ACD = ACD + DENS                                                     
                                                                                
           ACESQ = ACESQ + EN ** 2                                              
           ACKSQ = ACKSQ + KN ** 2                                              
           ACVSQ = ACVSQ + VN ** 2                                              
           ACPSQ = ACPSQ + PRES ** 2                                            
           ACHSQ = ACHSQ + HN ** 2                                              
           ACDSQ = ACDSQ + DENS ** 2                                            
                                                                                
C       ** OPTIONALLY PRINT INFORMATION **                                      
                                                                                
           IF ( MOD ( STEP, IPRINT ) .EQ. 0 ) THEN                              
                                                                                
              WRITE(*,'(1X,I8,7(2X,F10.5))')                                    
     :                 STEP, EN, HN, KN, VN, PRES, TEMP, DENS                   
                                                                                
           ENDIF                                                                
                                                                                
C       ** OPTIONALLY SCALE MOMENTA AND BOX SIZE **                             
                                                                                
           IF ( MOD ( STEP, ISCALE ) .EQ. 0 ) THEN                              
                                                                                
              CALL SCALE ( TEMPER, PRESUR, WLRC0, XLRC0, RCUT )                 
                                                                                
           ENDIF                                                                
                                                                                
1000    CONTINUE                                                                
                                                                                
C    *******************************************************************        
C    ** MAIN LOOP ENDS                                                **        
C    *******************************************************************        
                                                                                
        WRITE(*,'(/1X,''**** END OF DYNAMICS **** ''//)')                       
                                                                                
C    ** WRITE OUT FINAL CONFIGURATION **                                        
C    ** INCLUDING BOX,BOX1,BOX2,BOX3  **                                        
                                                                                
        CALL WRITCN ( CNFILE )                                                  
                                                                                
C    ** WRITE OUT FINAL AVERAGES **                                             
                                                                                
        NORM = REAL ( NSTEP )                                                   
        AVE  = ACE / NORM                                                       
        AVK  = ACK / NORM                                                       
        AVV  = ACV / NORM                                                       
        AVP  = ACP / NORM                                                       
        AVH  = ACH / NORM                                                       
        AVD  = ACD / NORM                                                       
                                                                                
        ACESQ = ( ACESQ / NORM ) - AVE ** 2                                     
        ACKSQ = ( ACKSQ / NORM ) - AVK ** 2                                     
        ACVSQ = ( ACVSQ / NORM ) - AVV ** 2                                     
        ACPSQ = ( ACPSQ / NORM ) - AVP ** 2                                     
        ACHSQ = ( ACHSQ / NORM ) - AVH ** 2                                     
        ACDSQ = ( ACDSQ / NORM ) - AVD ** 2                                     
                                                                                
        IF ( ACESQ .GT. 0.0 ) FLE = SQRT ( ACESQ )                              
        IF ( ACKSQ .GT. 0.0 ) FLK = SQRT ( ACKSQ )                              
        IF ( ACVSQ .GT. 0.0 ) FLV = SQRT ( ACVSQ )                              
        IF ( ACPSQ .GT. 0.0 ) FLP = SQRT ( ACPSQ )                              
        IF ( ACHSQ .GT. 0.0 ) FLH = SQRT ( ACHSQ )                              
        IF ( ACDSQ .GT. 0.0 ) FLD = SQRT ( ACDSQ )                              
                                                                                
        AVT = AVK / 1.5                                                         
        FLT = FLK / 1.5                                                         
                                                                                
        WRITE(*,'('' AVERAGES'',7(2X,F10.5))')                                  
     :             AVE, AVH, AVK, AVV, AVP, AVT, AVD                            
        WRITE(*,'('' FLUCTS  '',7(2X,F10.5))')                                  
     :             FLE, FLH, FLK, FLV, FLP, FLT, FLD                            
                                                                                
        STOP                                                                    
                                                                                
10001   FORMAT(//1X,'TIMESTEP  ..ENERGY..  .ENTHALPY.  ..KINETIC.',             
     :          '  ..POTENT..  .PRESSURE.  ..TEMPER..  ..DENSITY.')             
        END                                                                     
                                                                                
                                                                                
                                                                                
        SUBROUTINE FORCE ( RCUT, V, W, X, RPX, PF )                             
                                                                                
        COMMON / BLOCK1 / RX, RY, RZ, RX1, RY1, RZ1,                            
     :                    RX2, RY2, RZ2, RX3, RY3, RZ3,                         
     :                    PX, PY, PZ, PX1, PY1, PZ1,                            
     :                    PX2, PY2, PZ2, PX3, PY3, PZ3,                         
     :                    FX, FY, FZ                                            
        COMMON / BLOCK2 / BOX, BOX1, BOX2, BOX3                                 
                                                                                
C    *******************************************************************        
C    ** LENNARD-JONES FORCE ROUTINE IN REDUCED UNITS                  **        
C    **                                                               **        
C    ** THE POTENTIAL IS V(R) = 4*((1/R)**12-(1/R)**6)                **        
C    ** WE INCLUDE SPHERICAL CUTOFF AND MINIMUM IMAGING IN CUBIC BOX. **        
C    **                                                               **        
C    ** PRINCIPAL VARIABLES:                                          **        
C    **                                                               **        
C    ** INTEGER N                 NUMBER OF MOLECULES                 **        
C    ** REAL    RX(N),RY(N),RZ(N) MOLECULAR POSITIONS                 **        
C    ** REAL    PX(N),PY(N),PZ(N) MOLECULAR MOMENTA                   **        
C    ** REAL    FX(N),FY(N),FZ(N) MOLECULAR FORCES                    **        
C    ** REAL    BOX               SIMULATION BOX LENGTH               **        
C    ** REAL    RCUT              PAIR POTENTIAL CUTOFF               **        
C    ** REAL    V                 POTENTIAL ENERGY                    **        
C    ** REAL    W                 VIRIAL FUNCTION                     **        
C    ** REAL    X                 HYPERVIRIAL FUNCTION                **        
C    ** REAL    RPX               ADDITIONAL PAIR SUMMATION           **        
C    ** REAL    PF                ANOTHER PAIR SUMMATION              **        
C    *******************************************************************        
                                                                                
        INTEGER     N                                                           
        PARAMETER ( N = 108 )                                                   
                                                                                
        REAL        RX(N), RY(N), RZ(N)                                         
        REAL        RX1(N), RY1(N), RZ1(N)                                      
        REAL        RX2(N), RY2(N), RZ2(N)                                      
        REAL        RX3(N), RY3(N), RZ3(N)                                      
        REAL        PX(N), PY(N), PZ(N)                                         
        REAL        PX1(N), PY1(N), PZ1(N)                                      
        REAL        PX2(N), PY2(N), PZ2(N)                                      
        REAL        PX3(N), PY3(N), PZ3(N)                                      
        REAL        FX(N), FY(N), FZ(N)                                         
        REAL        BOX, BOX1, BOX2, BOX3                                       
                                                                                
        INTEGER     I, J                                                        
        REAL        RCUT, V, W, X, RPX, PF                                      
        REAL        BOXINV, RCUTSQ                                              
        REAL        RXI, RYI, RZI, RXIJ, RYIJ, RZIJ, RIJSQ                      
        REAL        PXI, PYI, PZI, PXIJ, PYIJ, PZIJ                             
        REAL        FXI, FYI, FZI, FXIJ, FYIJ, FZIJ                             
        REAL        SR2, SR6, SR12, VIJ, WIJ, FIJ, XIJ, RP                      
                                                                                
C    *******************************************************************        
                                                                                
C    ** USEFUL QUANTITIES **                                                    
                                                                                
        BOXINV = 1.0 / BOX                                                      
        RCUTSQ = RCUT ** 2                                                      
                                                                                
C    ** ZERO FORCES, POTENTIAL, VIRIAL **                                       
                                                                                
        DO 100 I = 1, N                                                         
                                                                                
           FX(I) = 0.0                                                          
           FY(I) = 0.0                                                          
           FZ(I) = 0.0                                                          
                                                                                
100     CONTINUE                                                                
                                                                                
        V   = 0.0                                                               
        W   = 0.0                                                               
        X   = 0.0                                                               
        RPX = 0.0                                                               
        PF  = 0.0                                                               
                                                                                
C    ** OUTER LOOP BEGINS **                                                    
                                                                                
        DO 200 I = 1, N - 1                                                     
                                                                                
           RXI = RX(I)                                                          
           RYI = RY(I)                                                          
           RZI = RZ(I)                                                          
           FXI = FX(I)                                                          
           FYI = FY(I)                                                          
           FZI = FZ(I)                                                          
           PXI = PX(I)                                                          
           PYI = PY(I)                                                          
           PZI = PZ(I)                                                          
                                                                                
C       ** INNER LOOP BEGINS **                                                 
                                                                                
           DO 199 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 ** 2 + RYIJ ** 2 + RZIJ ** 2                        
                                                                                
              IF ( RIJSQ .LT. RCUTSQ ) THEN                                     
                                                                                
                 SR2   = 1.0 / RIJSQ                                            
                 SR6   = SR2 * SR2 * SR2                                        
                 SR12  = SR6 ** 2                                               
                 VIJ   = SR12 - SR6                                             
                 V     = V + VIJ                                                
                 WIJ   = VIJ + SR12                                             
                 W     = W + WIJ                                                
                 FIJ   = WIJ * SR2                                              
                 FXIJ  = FIJ * RXIJ                                             
                 FYIJ  = FIJ * RYIJ                                             
                 FZIJ  = FIJ * RZIJ                                             
                 FXI   = FXI + FXIJ                                             
                 FYI   = FYI + FYIJ                                             
                 FZI   = FZI + FZIJ                                             
                 FX(J) = FX(J) - FXIJ                                           
                 FY(J) = FY(J) - FYIJ                                           
                 FZ(J) = FZ(J) - FZIJ                                           
                 XIJ   = WIJ + 2.0 * SR12                                       
                 X     = X + XIJ                                                
                                                                                
                 PXIJ  = PXI - PX(J)                                            
                 PYIJ  = PYI - PY(J)                                            
                 PZIJ  = PZI - PZ(J)                                            
                 RP    = RXIJ * PXIJ + RYIJ * PYIJ + RZIJ * PZIJ                
                 RPX   = RPX + RP * XIJ * SR2                                   
                 PF    = PF + PXIJ * FXIJ + PYIJ * FYIJ + PZIJ * FZIJ           
                                                                                
              ENDIF                                                             
                                                                                
199        CONTINUE                                                             
                                                                                
C       ** INNER LOOP ENDS **                                                   
                                                                                
           FX(I) = FXI                                                          
           FY(I) = FYI                                                          
           FZ(I) = FZI                                                          
                                                                                
200     CONTINUE                                                                
                                                                                
C    ** OUTER LOOP ENDS **                                                      
                                                                                
C    ** MULTIPLY RESULTS BY NUMERICAL FACTORS **                                
                                                                                
        DO 300 I = 1, N                                                         
                                                                                
           FX(I) = FX(I) * 24.0                                                 
           FY(I) = FY(I) * 24.0                                                 
           FZ(I) = FZ(I) * 24.0                                                 
                                                                                
300     CONTINUE                                                                
                                                                                
        V   = V * 4.0                                                           
        W   = W * 24.0 / 3.0                                                    
        X   = X * 144.0 / 9.0                                                   
        RPX = RPX * 144.0                                                       
        PF  = PF * 24.0                                                         
                                                                                
        RETURN                                                                  
        END                                                                     
                                                                                
                                                                                
                                                                                
        SUBROUTINE KINET ( K )                                                  
                                                                                
        COMMON / BLOCK1 / RX, RY, RZ, RX1, RY1, RZ1,                            
     :                    RX2, RY2, RZ2, RX3, RY3, RZ3,                         
     :                    PX, PY, PZ, PX1, PY1, PZ1,                            
     :                    PX2, PY2, PZ2, PX3, PY3, PZ3,                         
     :                    FX, FY, FZ                                            
        COMMON / BLOCK2 / BOX, BOX1, BOX2, BOX3                                 
                                                                                
C    *******************************************************************        
C    ** ROUTINE TO COMPUTE KINETIC ENERGY                             **        
C    *******************************************************************        
                                                                                
        INTEGER     N                                                           
        PARAMETER ( N = 108 )                                                   
                                                                                
        REAL        RX(N), RY(N), RZ(N)                                         
        REAL        RX1(N), RY1(N), RZ1(N)                                      
        REAL        RX2(N), RY2(N), RZ2(N)                                      
        REAL        RX3(N), RY3(N), RZ3(N)                                      
        REAL        PX(N), PY(N), PZ(N)                                         
        REAL        PX1(N), PY1(N), PZ1(N)                                      
        REAL        PX2(N), PY2(N), PZ2(N)                                      
        REAL        PX3(N), PY3(N), PZ3(N)                                      
        REAL        FX(N), FY(N), FZ(N)                                         
        REAL        BOX, BOX1, BOX2, BOX3                                       
                                                                                
        REAL        K                                                           
        INTEGER     I                                                           
                                                                                
C    *******************************************************************        
                                                                                
        K = 0.0                                                                 
                                                                                
        DO 1000 I = 1, N                                                        
                                                                                
           K = K + PX(I) ** 2 + PY(I) ** 2 + PZ(I) ** 2                         
                                                                                
1000    CONTINUE                                                                
                                                                                
        K = 0.5 * K                                                             
                                                                                
        RETURN                                                                  
        END                                                                     
                                                                                
                                                                                
                                                                                
        SUBROUTINE READCN ( CNFILE )                                            
                                                                                
        COMMON / BLOCK1 / RX, RY, RZ, RX1, RY1, RZ1,                            
     :                    RX2, RY2, RZ2, RX3, RY3, RZ3,                         
     :                    PX, PY, PZ, PX1, PY1, PZ1,                            
     :                    PX2, PY2, PZ2, PX3, PY3, PZ3,                         
     :                    FX, FY, FZ                                            
        COMMON / BLOCK2 / BOX, BOX1, BOX2, BOX3                                 
                                                                                
C    *******************************************************************        
C    ** SUBROUTINE TO READ IN INITIAL CONFIGURATION FROM UNIT 10      **        
C    *******************************************************************        
                                                                                
        INTEGER     N                                                           
        PARAMETER ( N = 108 )                                                   
                                                                                
        REAL        RX(N), RY(N), RZ(N)                                         
        REAL        RX1(N), RY1(N), RZ1(N)                                      
        REAL        RX2(N), RY2(N), RZ2(N)                                      
        REAL        RX3(N), RY3(N), RZ3(N)                                      
        REAL        PX(N), PY(N), PZ(N)                                         
        REAL        PX1(N), PY1(N), PZ1(N)                                      
        REAL        PX2(N), PY2(N), PZ2(N)                                      
        REAL        PX3(N), PY3(N), PZ3(N)                                      
        REAL        FX(N), FY(N), FZ(N)                                         
        REAL        BOX, BOX1, BOX2, BOX3                                       
                                                                                
        CHARACTER   CNFILE*(*)                                                  
                                                                                
        INTEGER     CNUNIT, NN                                                  
        PARAMETER ( CNUNIT = 10 )                                               
                                                                                
C     ******************************************************************        
                                                                                
        OPEN ( UNIT = CNUNIT, FILE = CNFILE,                                    
     :         STATUS = 'OLD', FORM = 'UNFORMATTED' )                           
                                                                                
        READ ( CNUNIT ) NN, BOX, BOX1, BOX2, BOX3                               
        IF ( NN .NE. N ) STOP 'INCORRECT VALUE OF N'                            
        READ ( CNUNIT ) RX, RY, RZ                                              
        READ ( CNUNIT ) RX1, RY1, RZ1                                           
        READ ( CNUNIT ) RX2, RY2, RZ2                                           
        READ ( CNUNIT ) RX3, RY3, RZ3                                           
        READ ( CNUNIT ) PX, PY, PZ                                              
        READ ( CNUNIT ) PX1, PY1, PZ1                                           
        READ ( CNUNIT ) PX2, PY2, PZ2                                           
        READ ( CNUNIT ) PX3, PY3, PZ3                                           
                                                                                
        CLOSE ( UNIT = CNUNIT )                                                 
                                                                                
        RETURN                                                                  
        END                                                                     
                                                                                
                                                                                
                                                                                
        SUBROUTINE WRITCN ( CNFILE )                                            
                                                                                
        COMMON / BLOCK1 / RX, RY, RZ, RX1, RY1, RZ1,                            
     :                    RX2, RY2, RZ2, RX3, RY3, RZ3,                         
     :                    PX, PY, PZ, PX1, PY1, PZ1,                            
     :                    PX2, PY2, PZ2, PX3, PY3, PZ3,                         
     :                    FX, FY, FZ                                            
        COMMON / BLOCK2 / BOX, BOX1, BOX2, BOX3                                 
                                                                                
C    *******************************************************************        
C    ** ROUTINE TO WRITE OUT FINAL CONFIGURATION TO UNIT 10           **        
C    *******************************************************************        
                                                                                
        INTEGER     N                                                           
        PARAMETER ( N = 108 )                                                   
                                                                                
        REAL        RX(N), RY(N), RZ(N)                                         
        REAL        RX1(N), RY1(N), RZ1(N)                                      
        REAL        RX2(N), RY2(N), RZ2(N)                                      
        REAL        RX3(N), RY3(N), RZ3(N)                                      
        REAL        PX(N), PY(N), PZ(N)                                         
        REAL        PX1(N), PY1(N), PZ1(N)                                      
        REAL        PX2(N), PY2(N), PZ2(N)                                      
        REAL        PX3(N), PY3(N), PZ3(N)                                      
        REAL        FX(N), FY(N), FZ(N)                                         
        REAL        BOX, BOX1, BOX2, BOX3                                       
                                                                                
        CHARACTER   CNFILE*(*)                                                  
                                                                                
        INTEGER     CNUNIT                                                      
        PARAMETER ( CNUNIT = 10 )                                               
                                                                                
C    *******************************************************************        
                                                                                
        OPEN ( UNIT = CNUNIT, FILE = CNFILE,                                    
     :         STATUS = 'OLD', FORM = 'UNFORMATTED' )                           
                                                                                
        WRITE ( CNUNIT ) N, BOX, BOX1, BOX2, BOX3                               
        WRITE ( CNUNIT ) RX, RY, RZ                                             
        WRITE ( CNUNIT ) RX1, RY1, RZ1                                          
        WRITE ( CNUNIT ) RX2, RY2, RZ2                                          
        WRITE ( CNUNIT ) RX3, RY3, RZ3                                          
        WRITE ( CNUNIT ) PX, PY, PZ                                             
        WRITE ( CNUNIT ) PX1, PY1, PZ1                                          
        WRITE ( CNUNIT ) PX2, PY2, PZ2                                          
        WRITE ( CNUNIT ) PX3, PY3, PZ3                                          
                                                                                
        CLOSE ( UNIT = CNUNIT )                                                 
                                                                                
        RETURN                                                                  
        END                                                                     
                                                                                
                                                                                
                                                                                
        SUBROUTINE PREDIC ( DT )                                                
                                                                                
        COMMON / BLOCK1 / RX, RY, RZ, RX1, RY1, RZ1,                            
     :                    RX2, RY2, RZ2, RX3, RY3, RZ3,                         
     :                    PX, PY, PZ, PX1, PY1, PZ1,                            
     :                    PX2, PY2, PZ2, PX3, PY3, PZ3,                         
     :                    FX, FY, FZ                                            
        COMMON / BLOCK2 / BOX, BOX1, BOX2, BOX3                                 
                                                                                
C    *******************************************************************        
C    ** STANDARD TAYLOR SERIES PREDICTORS                             **        
C    *******************************************************************        
                                                                                
        INTEGER     N                                                           
        PARAMETER ( N = 108 )                                                   
                                                                                
        REAL        RX(N), RY(N), RZ(N)                                         
        REAL        RX1(N), RY1(N), RZ1(N)                                      
        REAL        RX2(N), RY2(N), RZ2(N)                                      
        REAL        RX3(N), RY3(N), RZ3(N)                                      
        REAL        PX(N), PY(N), PZ(N)                                         
        REAL        PX1(N), PY1(N), PZ1(N)                                      
        REAL        PX2(N), PY2(N), PZ2(N)                                      
        REAL        PX3(N), PY3(N), PZ3(N)                                      
        REAL        FX(N), FY(N), FZ(N)                                         
        REAL        BOX, BOX1, BOX2, BOX3                                       
        REAL        DT                                                          
                                                                                
        INTEGER     I                                                           
        REAL        C1, C2, C3                                                  
                                                                                
C    *******************************************************************        
                                                                                
        C1 = DT                                                                 
        C2 = C1 * DT / 2.0                                                      
        C3 = C2 * DT / 3.0                                                      
                                                                                
        DO 100 I = 1, N                                                         
                                                                                
           RX(I)  = RX(I)  + C1*RX1(I) + C2*RX2(I) + C3*RX3(I)                  
           RY(I)  = RY(I)  + C1*RY1(I) + C2*RY2(I) + C3*RY3(I)                  
           RZ(I)  = RZ(I)  + C1*RZ1(I) + C2*RZ2(I) + C3*RZ3(I)                  
           RX1(I) = RX1(I) + C1*RX2(I) + C2*RX3(I)                              
           RY1(I) = RY1(I) + C1*RY2(I) + C2*RY3(I)                              
           RZ1(I) = RZ1(I) + C1*RZ2(I) + C2*RZ3(I)                              
           RX2(I) = RX2(I) + C1*RX3(I)                                          
           RY2(I) = RY2(I) + C1*RY3(I)                                          
           RZ2(I) = RZ2(I) + C1*RZ3(I)                                          
                                                                                
           PX(I)  = PX(I)  + C1*PX1(I) + C2*PX2(I) + C3*PX3(I)                  
           PY(I)  = PY(I)  + C1*PY1(I) + C2*PY2(I) + C3*PY3(I)                  
           PZ(I)  = PZ(I)  + C1*PZ1(I) + C2*PZ2(I) + C3*PZ3(I)                  
           PX1(I) = PX1(I) + C1*PX2(I) + C2*PX3(I)                              
           PY1(I) = PY1(I) + C1*PY2(I) + C2*PY3(I)                              
           PZ1(I) = PZ1(I) + C1*PZ2(I) + C2*PZ3(I)                              
           PX2(I) = PX2(I) + C1*PX3(I)                                          
           PY2(I) = PY2(I) + C1*PY3(I)                                          
           PZ2(I) = PZ2(I) + C1*PZ3(I)                                          
                                                                                
100     CONTINUE                                                                
                                                                                
        BOX  = BOX  + C1*BOX1 + C2*BOX2 + C3*BOX3                               
        BOX1 = BOX1 + C1*BOX2 + C2*BOX3                                         
        BOX2 = BOX2 + C1*BOX3                                                   
                                                                                
        RETURN                                                                  
        END                                                                     
                                                                                
                                                                                
                                                                                
        SUBROUTINE CORREC ( DT, CHI, CHIPXI )                                   
                                                                                
        COMMON / BLOCK1 / RX, RY, RZ, RX1, RY1, RZ1,                            
     :                    RX2, RY2, RZ2, RX3, RY3, RZ3,                         
     :                    PX, PY, PZ, PX1, PY1, PZ1,                            
     :                    PX2, PY2, PZ2, PX3, PY3, PZ3,                         
     :                    FX, FY, FZ                                            
        COMMON / BLOCK2 / BOX, BOX1, BOX2, BOX3                                 
                                                                                
C    *******************************************************************        
C    ** GEAR CORRECTOR ALGORITHM.                                     **        
C    **                                                               **        
C    ** FOR TIMESTEP-SCALED VARIABLES GEAR COEFFICIENTS WOULD BE AS   **        
C    ** FOLLOWS (4-VALUE METHOD, 1ST-ORDER D.E.): 3/8, 1, 3/4, 1/6.   **        
C    ** CHIPXI IS SHORT FOR ( CHI + XI ) (SEE CHAPTER 7).             **        
C    *******************************************************************        
                                                                                
        INTEGER     N                                                           
        PARAMETER ( N = 108 )                                                   
                                                                                
        REAL        RX(N), RY(N), RZ(N)                                         
        REAL        RX1(N), RY1(N), RZ1(N)                                      
        REAL        RX2(N), RY2(N), RZ2(N)                                      
        REAL        RX3(N), RY3(N), RZ3(N)                                      
        REAL        PX(N), PY(N), PZ(N)                                         
        REAL        PX1(N), PY1(N), PZ1(N)                                      
        REAL        PX2(N), PY2(N), PZ2(N)                                      
        REAL        PX3(N), PY3(N), PZ3(N)                                      
        REAL        FX(N), FY(N), FZ(N)                                         
        REAL        BOX, BOX1, BOX2, BOX3                                       
        REAL        DT, CHI, CHIPXI                                             
                                                                                
        INTEGER     I                                                           
        REAL        C1, C2, C3, COEFF0, COEFF2, COEFF3                          
        REAL        CORR, CORRX, CORRY, CORRZ, CORPX, CORPY, CORPZ              
        REAL        RX1I, RY1I, RZ1I, PX1I, PY1I, PZ1I                          
                                                                                
        REAL        GEAR0, GEAR2, GEAR3                                         
        PARAMETER ( GEAR0 = 3.0 / 8.0,                                          
     :              GEAR2 = 3.0 / 4.0,                                          
     :              GEAR3 = 1.0 / 6.0 )                                         
                                                                                
C    *******************************************************************        
                                                                                
        C1 = DT                                                                 
        C2 = C1 * DT / 2.0                                                      
        C3 = C2 * DT / 3.0                                                      
                                                                                
        COEFF0 = GEAR0 * C1                                                     
        COEFF2 = GEAR2 * C1 / C2                                                
        COEFF3 = GEAR3 * C1 / C3                                                
                                                                                
        DO 400 I = 1, N                                                         
                                                                                
           RX1I = PX(I) + CHI * RX(I)                                           
           RY1I = PY(I) + CHI * RY(I)                                           
           RZ1I = PZ(I) + CHI * RZ(I)                                           
           CORRX = RX1I - RX1(I)                                                
           CORRY = RY1I - RY1(I)                                                
           CORRZ = RZ1I - RZ1(I)                                                
                                                                                
           RX(I)  = RX(I)  + COEFF0 * CORRX                                     
           RY(I)  = RY(I)  + COEFF0 * CORRY                                     
           RZ(I)  = RZ(I)  + COEFF0 * CORRZ                                     
           RX1(I) = RX1I                                                        
           RY1(I) = RY1I                                                        
           RZ1(I) = RZ1I                                                        
           RX2(I) = RX2(I) + COEFF2 * CORRX                                     
           RY2(I) = RY2(I) + COEFF2 * CORRY                                     
           RZ2(I) = RZ2(I) + COEFF2 * CORRZ                                     
           RX3(I) = RX3(I) + COEFF3 * CORRX                                     
           RY3(I) = RY3(I) + COEFF3 * CORRY                                     
           RZ3(I) = RZ3(I) + COEFF3 * CORRZ                                     
                                                                                
           PX1I = FX(I) - CHIPXI * PX(I)                                        
           PY1I = FY(I) - CHIPXI * PY(I)                                        
           PZ1I = FZ(I) - CHIPXI * PZ(I)                                        
           CORPX = PX1I - PX1(I)                                                
           CORPY = PY1I - PY1(I)                                                
           CORPZ = PZ1I - PZ1(I)                                                
                                                                                
           PX(I)  = PX(I)  + COEFF0 * CORPX                                     
           PY(I)  = PY(I)  + COEFF0 * CORPY                                     
           PZ(I)  = PZ(I)  + COEFF0 * CORPZ                                     
           PX1(I) = PX1I                                                        
           PY1(I) = PY1I                                                        
           PZ1(I) = PZ1I                                                        
           PX2(I) = PX2(I) + COEFF2 * CORPX                                     
           PY2(I) = PY2(I) + COEFF2 * CORPY                                     
           PZ2(I) = PZ2(I) + COEFF2 * CORPZ                                     
           PX3(I) = PX3(I) + COEFF3 * CORPX                                     
           PY3(I) = PY3(I) + COEFF3 * CORPY                                     
           PZ3(I) = PZ3(I) + COEFF3 * CORPZ                                     
                                                                                
400     CONTINUE                                                                
                                                                                
        CORR = CHI * BOX - BOX1                                                 
        BOX  = BOX  + COEFF0 * CORR                                             
        BOX1 = CHI * BOX                                                        
        BOX2 = BOX2 + COEFF2 * CORR                                             
        BOX3 = BOX3 + COEFF3 * CORR                                             
                                                                                
        RETURN                                                                  
        END                                                                     
                                                                                
                                                                                
                                                                                
        SUBROUTINE SCALE ( TEMPER, PRESUR, WLRC0, XLRC0, RCUT )                 
                                                                                
        COMMON / BLOCK1 / RX, RY, RZ, RX1, RY1, RZ1,                            
     :                    RX2, RY2, RZ2, RX3, RY3, RZ3,                         
     :                    PX, PY, PZ, PX1, PY1, PZ1,                            
     :                    PX2, PY2, PZ2, PX3, PY3, PZ3,                         
     :                    FX, FY, FZ                                            
        COMMON / BLOCK2 / BOX, BOX1, BOX2, BOX3                                 
                                                                                
C    *******************************************************************        
C    ** SCALES MOMENTA AND BOX SIZE TO GIVE DESIRED VALUES OF T AND P.**        
C    **                                                               **        
C    ** WE USE DIRECT MOMENTUM SCALING TO GIVE THE TEMPERATURE AND A  **        
C    ** NEWTON-RAPHSON METHOD FOR THE BOX SIZE.                       **        
C    ** WE USE THE PRESSURE DERIVATIVE FUNCTION                       **        
C    ** BOX * DP/D(BOX) = 3 * VOL * DP/D(VOL) = 3 * ( - P - X / VOL ) **        
C    ** WHERE P IS THE PRESSURE AND X THE HYPERVIRIAL FUNCTION.       **        
C    ** LONG-RANGE CORRECTIONS ARE TAKEN INTO ACCOUNT.                **        
C    *******************************************************************        
                                                                                
        INTEGER     N                                                           
        PARAMETER ( N = 108 )                                                   
                                                                                
        REAL        RX(N), RY(N), RZ(N)                                         
        REAL        RX1(N), RY1(N), RZ1(N)                                      
        REAL        RX2(N), RY2(N), RZ2(N)                                      
        REAL        RX3(N), RY3(N), RZ3(N)                                      
        REAL        PX(N), PY(N), PZ(N)                                         
        REAL        PX1(N), PY1(N), PZ1(N)                                      
        REAL        PX2(N), PY2(N), PZ2(N)                                      
        REAL        PX3(N), PY3(N), PZ3(N)                                      
        REAL        FX(N), FY(N), FZ(N)                                         
        REAL        BOX, BOX1, BOX2, BOX3                                       
                                                                                
        INTEGER     I                                                           
        REAL        TEMPER, PRESUR, WLRC0, XLRC0, RCUT                          
        REAL        K, KN, V, W, X, RPX, PF, FACTOR                             
        REAL        VOL, DENS, WLRC, XLRC, TEMP, PRES                           
        REAL        FREE, TOL                                                   
        PARAMETER ( FREE = 3.0, TOL = 0.01 )                                    
                                                                                
C    *******************************************************************        
                                                                                
C    ** SCALE MOMENTA **                                                        
                                                                                
        CALL KINET ( K )                                                        
        KN = K / REAL ( N )                                                     
        TEMP = 2.0 * KN / FREE                                                  
        FACTOR = SQRT ( TEMPER / TEMP )                                         
                                                                                
        DO 100 I = 1, N                                                         
                                                                                
           PX(I) = PX(I) * FACTOR                                               
           PY(I) = PY(I) * FACTOR                                               
           PZ(I) = PZ(I) * FACTOR                                               
                                                                                
100     CONTINUE                                                                
                                                                                
        TEMP = TEMPER                                                           
                                                                                
C    ** SCALE PRESSURE USING NEWTON-RAPHSON PROCEDURE **                        
                                                                                
        CALL FORCE ( RCUT, V, W, X, RPX, PF )                                   
        VOL  = BOX ** 3                                                         
        DENS = REAL ( N ) / VOL                                                 
        WLRC = WLRC0 * DENS                                                     
        XLRC = XLRC0 * DENS                                                     
        W    = W + WLRC                                                         
        X    = X + XLRC                                                         
        PRES = DENS * TEMP + W / VOL                                            
                                                                                
1000    IF ( ABS ( PRES - PRESUR ) .GT. TOL ) THEN                              
                                                                                
           FACTOR = 1.0 + ( PRES - PRESUR ) / ( PRES + X / VOL ) / 3.0          
                                                                                
           DO 200 I = 1, N                                                      
                                                                                
              RX(I) = RX(I) * FACTOR                                            
              RY(I) = RY(I) * FACTOR                                            
              RZ(I) = RZ(I) * FACTOR                                            
                                                                                
200        CONTINUE                                                             
                                                                                
           BOX = BOX * FACTOR                                                   
                                                                                
           CALL FORCE ( RCUT, V, W, X, RPX, PF )                                
           VOL  = BOX ** 3                                                      
           DENS = REAL ( N ) / VOL                                              
           WLRC = WLRC0 * DENS                                                  
           XLRC = XLRC0 * DENS                                                  
           W    = W + WLRC                                                      
           X    = X + XLRC                                                      
           PRES = DENS * TEMP + W / VOL                                         
                                                                                
           GOTO 1000                                                            
                                                                                
        ENDIF                                                                   
                                                                                
        RETURN                                                                  
        END                                                                     
