C    *******************************************************************        
C    ** FICHE F.28                                                    **        
C    ** CONSTANT-NVT MOLECULAR DYNAMICS USING AN EXTENDED SYSTEM.     **        
C    *******************************************************************        
                                                                                
        PROGRAM NOSE                                                            
                                                                                
        COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ, AX, AY, AZ,                   
     :                    BX, BY, BZ, CX, CY, CZ, FX, FY, FZ                    
        COMMON / BLOCK2 / S, SV, SA, SB, SC, SF                                 
                                                                                
C    *******************************************************************        
C    ** CONSTANT-NVT MOLECULAR DYNAMICS USING AN EXTENDED SYSTEM.     **        
C    **                                                               **        
C    ** REFERENCE:                                                    **        
C    **                                                               **        
C    ** NOSE, MOLEC. PHYS. 52, 255, 1984.                             **        
C    **                                                               **        
C    ** ROUTINES REFERENCED:                                          **        
C    **                                                               **        
C    ** SUBROUTINE FORCE ( BOX, RCUT, V, VC, W )                      **        
C    **    CALCULATES ACCELERATIONS, POTENTIAL AND VIRIAL             **        
C    ** SUBROUTINE KINET ( K )                                        **        
C    **    CALCULATES KINETIC ENERGY                                  **        
C    ** SUBROUTINE READCN ( CNFILE, BOX )                             **        
C    **    READS CONFIGURATION                                        **        
C    ** SUBROUTINE WRITCN ( CNFILE, BOX )                             **        
C    **    WRITES CONFIGURATION                                       **        
C    ** SUBROUTINE PREDIC ( DT )                                      **        
C    **    PREDICTS POSITIONS AT TIME T + DT                          **        
C    ** SUBROUTINE CORREC ( DT )                                      **        
C    **    CORRECTS POSITIONS AT TIME T + DT                          **        
C    **                                                               **        
C    ** PRINCIPAL VARIABLES:                                          **        
C    **                                                               **        
C    ** INTEGER N                             NUMBER OF ATOMS         **        
C    ** REAL    DT                            TIMESTEP                **        
C    ** REAL    RX(N),RY(N),RZ(N)             ATOM POSITIONS          **        
C    ** REAL    VX(N),VY(N),VZ(N)             FIRST DERIVATIVES       **        
C    ** REAL    AX(N),AY(N),AZ(N)             SECOND DERIVATIVES      **        
C    ** REAL    BX(N),BY(N),BZ(N)             THIRD DERIVATIVES       **        
C    ** REAL    FX(N),FY(N),FZ(N)             TOTAL FORCES            **        
C    ** REAL    S,SV,SA,SB,SC,SF              S AND DERIVATIVES       **        
C    ** REAL    Q                             THERMAL INERTIA         **        
C    ** REAL    ACV,ACK ETC.                  AV VALUE ACCUMULATORS   **        
C    ** REAL    AVV,AVK ETC.                  AV VALUES               **        
C    ** REAL    ACVSQ, ACKSQ ETC.             AV**2 VAL ACCUMULATORS  **        
C    ** REAL    FLV, FLK ETC.                 FLUCT AVERAGES          **        
C    **                                                               **        
C    ** USAGE:                                                        **        
C    **                                                               **        
C    ** THE MODIFIED EQUATIONS OF MOTION ARE AS FOLLOWS:              **        
C    ** A = F/(M*S**2) - 2*SV*V/S                                     **        
C    ** SA = (SUM (M*S/Q)*V**2 ) - ((3N-3)+1)*KT/(Q*S)                **        
C    ** WHERE R,V,A,.. ARE SUCCESSIVE DERIVATIVES OF POSITION         **        
C    ** AND SV,SA,.... ARE SUCCESSIVE DERIVATIVES OF EXTRA VARIABLE S **        
C    ** F REPRESENTS FORCES, M PARTICLE MASS, Q EXTRA VARIABLE MASS   **        
C    ** KT IS TEMPERATURE*BOLTZMANN'S CONSTANT                        **        
C    ** AND WE HAVE (3N-3) DEGREES OF FREEDOM IF MOMENTUM FIXED.      **        
C    ** IN A SENSE, S IS A TIME-SCALING FACTOR.                       **        
C    ** WE SOLVE THESE EQUATIONS BY A GEAR 5-VALUE METHOD FOR SECOND  **        
C    ** ORDER DIFFERENTIAL EQUATIONS.                                 **        
C    ** THIS PROGRAM USES UNIT 10 FOR CONFIGURATION INPUT AND OUTPUT  **        
C    *******************************************************************        
                                                                                
        INTEGER       N                                                         
        PARAMETER   ( N = 108 )                                                 
                                                                                
        REAL          RX(N), RY(N), RZ(N), VX(N), VY(N), VZ(N)                  
        REAL          AX(N), AY(N), AZ(N), BX(N), BY(N), BZ(N)                  
        REAL          CX(N), CY(N), CZ(N), FX(N), FY(N), FZ(N)                  
        REAL          S, SV, SA, SB, SC, SF                                     
                                                                                
        INTEGER       STEP, NSTEP, IPRINT                                       
        REAL          Q, QK, QV                                                 
        REAL          ACV, ACK, ACE, ACEC, ACP, ACT, ACH, ACS                   
        REAL          AVV, AVK, AVE, AVEC, AVP, AVT, AVH, AVS                   
        REAL          ACVSQ, ACKSQ, ACESQ, ACECSQ, ACPSQ, ACTSQ                 
        REAL          ACHSQ, ACSSQ                                              
        REAL          FLV, FLK, FLE, FLEC, FLP, FLT, FLH, FLS                   
        REAL          DT, DENS, TEMPER, RCUT, BOX, PRES, NORM, TEMP             
        REAL          K, V, VC, W, E, EC, FREE, FREEN, H, VOL                   
        REAL          KN, VN, EN, ECN, HN                                       
        REAL          SR3, SR9, VLRC, WLRC, PI                                  
        CHARACTER     TITLE*80, CNFILE*30                                       
                                                                                
        PARAMETER   ( PI = 3.1415927 )                                          
        PARAMETER   ( FREE = 3.0 )                                              
                                                                                
C    *******************************************************************        
                                                                                
C    ** READ IN INITIAL PARAMETERS **                                           
                                                                                
        WRITE(*,'(1H1,'' **** PROGRAM NOSE ****                    '')')        
        WRITE(*,'(//1X,''MOLECULAR DYNAMICS OF LENNARD-JONES ATOMS '')')        
        WRITE(*,'('' CONSTANT-NVT DYNAMICS BY THE METHOD OF NOSE   '')')        
        WRITE(*,'('' ENTER RUN TITLE                               '')')        
        READ (*,'(A)') TITLE                                                    
        WRITE(*,'('' ENTER NUMBER OF STEPS                         '')')        
        READ (*,*) NSTEP                                                        
        WRITE(*,'('' ENTER INTERVAL BETWEEN PRINTS                 '')')        
        READ (*,*) IPRINT                                                       
        WRITE(*,'('' ENTER CONFIGURATION FILENAME                  '')')        
        READ (*,'(A)') CNFILE                                                   
        WRITE(*,'('' ENTER THE FOLLOWING IN L-J REDUCED UNITS      '')')        
        WRITE(*,'('' ENTER TIMESTEP                                '')')        
        READ (*,*) DT                                                           
        WRITE(*,'('' ENTER TEMPERATURE                             '')')        
        READ (*,*) TEMPER                                                       
        WRITE(*,'('' ENTER POTENTIAL CUTOFF                        '')')        
        READ (*,*) RCUT                                                         
        WRITE(*,'('' ENTER THERMAL INERTIA PARAMETER               '')')        
        READ (*,*) Q                                                            
                                                                                
        WRITE(*,'(//1X,A)') TITLE                                               
        WRITE(*,'('' NUMBER OF ATOMS  = '',I6   )') N                           
        WRITE(*,'('' NUMBER OF STEPS  = '',I6   )') NSTEP                       
        WRITE(*,'('' PRINT INTERVAL   = '',I6   )') IPRINT                      
        WRITE(*,'('' CONFIGURATION FILENAME '',A)') CNFILE                      
        WRITE(*,'('' TIMESTEP         = '',F10.5)') DT                          
        WRITE(*,'('' TEMPERATURE      = '',F10.5)') TEMPER                      
        WRITE(*,'('' POTENTIAL CUTOFF = '',F10.5)') RCUT                        
        WRITE(*,'('' THERMAL INERTIA  = '',F10.2)') Q                           
                                                                                
        FREEN = FREE * REAL ( N )                                               
                                                                                
C    ** READCN MUST READ IN INITIAL CONFIGURATION  **                           
C    ** AND ASSIGN VALUES TO S AND ITS DERIVATIVES **                           
                                                                                
        CALL READCN ( CNFILE, BOX )                                             
                                                                                
        VOL  = BOX ** 3                                                         
        DENS = REAL ( N ) / VOL                                                 
                                                                                
        WRITE(*,'('' DENSITY          = '',F10.5)') DENS                        
                                                                                
C    ** CALCULATE LONG-RANGE CORRECTIONS **                                     
C    ** NB: SPECIFIC TO LENNARD-JONES    **                                     
                                                                                
        SR3  = ( 1.0 / RCUT ) ** 3                                              
        SR9  = SR3 ** 3                                                         
        VLRC = ( 8.0 / 9.0  ) * PI * DENS * REAL ( N ) *                        
     :         ( SR9 - 3.0 * SR3 )                                              
        WLRC = ( 16.0 / 9.0 ) * PI * DENS * REAL ( N ) *                        
     :         ( 2.0 * SR9 - 3.0 * SR3 )                                        
                                                                                
C    ** ZERO ACCUMULATORS **                                                    
                                                                                
        ACV  = 0.0                                                              
        ACK  = 0.0                                                              
        ACE  = 0.0                                                              
        ACEC = 0.0                                                              
        ACP  = 0.0                                                              
        ACT  = 0.0                                                              
        ACH  = 0.0                                                              
        ACS  = 0.0                                                              
                                                                                
        ACVSQ  = 0.0                                                            
        ACKSQ  = 0.0                                                            
        ACESQ  = 0.0                                                            
        ACECSQ = 0.0                                                            
        ACPSQ  = 0.0                                                            
        ACTSQ  = 0.0                                                            
        ACHSQ  = 0.0                                                            
        ACSSQ  = 0.0                                                            
                                                                                
        FLV  = 0.0                                                              
        FLK  = 0.0                                                              
        FLE  = 0.0                                                              
        FLEC = 0.0                                                              
        FLP  = 0.0                                                              
        FLT  = 0.0                                                              
        FLH  = 0.0                                                              
        FLS  = 0.0                                                              
                                                                                
        IF ( IPRINT .LE. 0 ) IPRINT = NSTEP + 1                                 
                                                                                
        WRITE(*,'(//1X,''**** START OF DYNAMICS ****'')')                       
        WRITE(*,10001)                                                          
                                                                                
C    *******************************************************************        
C    ** MAIN LOOP BEGINS                                              **        
C    *******************************************************************        
                                                                                
        DO 1000 STEP = 1, NSTEP                                                 
                                                                                
C       ** IMPLEMENT ALGORITHM **                                               
                                                                                
           CALL PREDIC ( DT )                                                   
           CALL FORCE  ( BOX, RCUT, V, VC, W )                                  
           CALL KINET  ( K )                                                    
                                                                                
           SF = ( 2.0 * K - ( FREEN + 1.0 ) * TEMPER ) / ( S * Q )              
                                                                                
           CALL CORREC ( DT )                                                   
                                                                                
           QK   = 0.5 * Q * ( SV ** 2 )                                         
           QV   = ( FREEN + 1.0 ) * TEMPER * LOG ( S )                          
           V    = V + VLRC                                                      
           W    = W + WLRC                                                      
           E    = K + V                                                         
           EC   = K + VC                                                        
           H    = K + VC + QK + QV                                              
           EN   = E / REAL ( N )                                                
           VN   = V / REAL ( N )                                                
           ECN  = EC / REAL ( N )                                               
           HN   = H / REAL ( N )                                                
           TEMP = KN * 2.0 / FREE                                               
           PRES = DENS * TEMP + W / VOL                                         
                                                                                
C       ** INCREMENT ACCUMULATORS **                                            
                                                                                
           ACE  = ACE  + EN                                                     
           ACEC = ACEC + ECN                                                    
           ACK  = ACK  + KN                                                     
           ACV  = ACV  + VN                                                     
           ACP  = ACP  + PRES                                                   
           ACH  = ACH  + HN                                                     
           ACS  = ACS  + S                                                      
                                                                                
           ACESQ  = ACESQ  + EN ** 2                                            
           ACECSQ = ACECSQ + ECN ** 2                                           
           ACKSQ  = ACKSQ  + KN ** 2                                            
           ACVSQ  = ACVSQ  + VN ** 2                                            
           ACPSQ  = ACPSQ  + PRES ** 2                                          
           ACHSQ  = ACHSQ  + HN ** 2                                            
           ACSSQ  = ACSSQ  + S ** 2                                             
                                                                                
C       ** OPTIONALLY PRINT INFORMATION **                                      
                                                                                
           IF ( MOD ( STEP, IPRINT ) .EQ. 0 ) THEN                              
                                                                                
              WRITE(*,'(1X,I8,8(2X,F10.4))')                                    
     :           STEP, HN, EN, ECN, KN, VN, TEMP, PRES, S                       
                                                                                
           ENDIF                                                                
                                                                                
1000    CONTINUE                                                                
                                                                                
C    *******************************************************************        
C    ** MAIN LOOP ENDS                                                **        
C    *******************************************************************        
                                                                                
        WRITE(*,'(/1X,''**** END OF DYNAMICS **** ''//)')                       
                                                                                
C    ** WRITE OUT FINAL AVERAGES **                                             
                                                                                
        NORM = REAL ( NSTEP )                                                   
                                                                                
        AVE  = ACE  / NORM                                                      
        AVEC = ACEC / NORM                                                      
        AVK  = ACK  / NORM                                                      
        AVV  = ACV  / NORM                                                      
        AVP  = ACP  / NORM                                                      
        AVH  = ACH  / NORM                                                      
        AVS  = ACS  / NORM                                                      
                                                                                
        ACESQ  = ( ACESQ  / NORM ) - AVE  ** 2                                  
        ACECSQ = ( ACECSQ / NORM ) - AVEC ** 2                                  
        ACKSQ  = ( ACKSQ  / NORM ) - AVK  ** 2                                  
        ACVSQ  = ( ACVSQ  / NORM ) - AVV  ** 2                                  
        ACPSQ  = ( ACPSQ  / NORM ) - AVP  ** 2                                  
        ACHSQ  = ( ACHSQ  / NORM ) - AVH  ** 2                                  
        ACSSQ  = ( ACSSQ  / NORM ) - AVS  ** 2                                  
                                                                                
        IF ( ACESQ  .GT. 0.0 ) FLE  = SQRT ( ACESQ  )                           
        IF ( ACECSQ .GT. 0.0 ) FLEC = SQRT ( ACECSQ )                           
        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 ( ACSSQ  .GT. 0.0 ) FLS  = SQRT ( ACSSQ  )                           
                                                                                
        AVT = AVK * 2.0 / FREE                                                  
        FLT = FLK * 2.0 / FREE                                                  
                                                                                
        WRITE(*,'('' AVERAGES'',8(2X,F10.5))')                                  
     :     AVH, AVE, AVEC, AVK, AVV, AVT, AVP, AVS                              
        WRITE(*,'('' FLUCTS  '',8(2X,F10.5))')                                  
     :     FLH, FLE, FLEC, FLK, FLV, FLT, FLP, FLS                              
                                                                                
C    ** WRITE OUT FINAL CONFIGURATION **                                        
                                                                                
        CALL WRITCN ( CNFILE, BOX )                                             
                                                                                
        STOP                                                                    
                                                                                
10001   FORMAT(//1X,'TIMESTEP    ...HAMIL..  ..ENERGY..',                       
     :              '  CUTENERGY.  ..KINETIC.  ..POTENT..',                     
     :              '  ..TEMPER..  .PRESSURE.  ....S.....'/)                    
        END                                                                     
                                                                                
                                                                                
                                                                                
        SUBROUTINE FORCE ( BOX, RCUT, V, VC, W )                                
                                                                                
        COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ, AX, AY, AZ,                   
     :                    BX, BY, BZ, CX, CY, CZ, FX, FY, FZ                    
                                                                                
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    ** TWO POTENTIAL ENERGIES ARE RETURNED.                          **        
C    ** V IS CALCULATED USING THE UNSHIFTED POTENTIAL. WHEN LONG-     **        
C    ** RANGE TAIL CORRECTIONS ARE ADDED, THIS MAY BE USED TO         **        
C    ** CALCULATE THERMODYNAMIC INTERNAL ENERGY ETC.                  **        
C    ** VC IS CALCULATED USING THE SHIFTED POTENTIAL WITH NO          **        
C    ** DISCONTINUITY AT THE CUTOFF.  THIS MAY BE USED TO CHECK THE   **        
C    ** CONSERVATION LAWS.                                            **        
C    **                                                               **        
C    ** PRINCIPAL VARIABLES:                                          **        
C    **                                                               **        
C    ** INTEGER N                 NUMBER OF MOLECULES                 **        
C    ** REAL    RX(N),RY(N),RZ(N) MOLECULAR POSITIONS                 **        
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    VC                SHIFTED POTENTIAL                   **        
C    ** REAL    W                 VIRIAL FUNCTION                     **        
C    ** REAL    VIJ               PAIR POTENTIAL BETWEEN I AND J      **        
C    ** REAL    WIJ               NEGATIVE OF PAIR VIRIAL FUNCTION W  **        
C    *******************************************************************        
                                                                                
        INTEGER     N                                                           
        PARAMETER ( N = 108 )                                                   
                                                                                
        REAL        RCUT, BOX, V, VC, W                                         
        REAL        RX(N), RY(N), RZ(N), VX(N), VY(N), VZ(N)                    
        REAL        AX(N), AY(N), AZ(N), BX(N), BY(N), BZ(N)                    
        REAL        CX(N), CY(N), CZ(N), FX(N), FY(N), FZ(N)                    
                                                                                
        INTEGER     I, J, NCUT                                                  
        REAL        BOXINV, RCUTSQ                                              
        REAL        RXI, RYI, RZI, FXI, FYI, FZI                                
        REAL        RXIJ, RYIJ, RZIJ, FXIJ, FYIJ, FZIJ                          
        REAL        RIJSQ, SR2, SR6, SR12, VIJ, WIJ, FIJ                        
                                                                                
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                                                                
                                                                                
        NCUT = 0                                                                
        V    = 0.0                                                              
        W    = 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)                                                          
                                                                                
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                                           
                 NCUT  = NCUT + 1                                               
                                                                                
              ENDIF                                                             
                                                                                
199        CONTINUE                                                             
                                                                                
C       ** INNER LOOP ENDS **                                                   
                                                                                
           FX(I) = FXI                                                          
           FY(I) = FYI                                                          
           FZ(I) = FZI                                                          
                                                                                
200     CONTINUE                                                                
                                                                                
C    ** OUTER LOOP ENDS **                                                      
                                                                                
C    ** CALCULATE SHIFTED POTENTIAL **                                          
                                                                                
        SR2  = 1.0 / RCUTSQ                                                     
        SR6  = SR2 * SR2 * SR2                                                  
        SR12 = SR6 * SR6                                                        
        VIJ  = SR12 - SR6                                                       
        VC   = V - REAL ( NCUT ) * VIJ                                          
                                                                                
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                                                           
        VC = VC * 4.0                                                           
        W  = W  * 24.0 / 3.0                                                    
                                                                                
        RETURN                                                                  
        END                                                                     
                                                                                
                                                                                
                                                                                
        SUBROUTINE KINET ( K )                                                  
                                                                                
        COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ, AX, AY, AZ,                   
     :                    BX, BY, BZ, CX, CY, CZ, FX, FY, FZ                    
        COMMON / BLOCK2 / S, SV, SA, SB, SC, SF                                 
                                                                                
C    *******************************************************************        
C    ** ROUTINE TO COMPUTE KINETIC ENERGY IN NOSE METHOD              **        
C    *******************************************************************        
                                                                                
        INTEGER     N                                                           
        PARAMETER ( N = 108 )                                                   
                                                                                
        REAL        K                                                           
                                                                                
        REAL        RX(N), RY(N), RZ(N), VX(N), VY(N), VZ(N)                    
        REAL        AX(N), AY(N), AZ(N), BX(N), BY(N), BZ(N)                    
        REAL        CX(N), CY(N), CZ(N), FX(N), FY(N), FZ(N)                    
        REAL        S, SV, SA, SB, SC, SF                                       
                                                                                
        INTEGER     I                                                           
                                                                                
C    *******************************************************************        
                                                                                
        K = 0.0                                                                 
                                                                                
        DO 1000 I = 1, N                                                        
                                                                                
           K = K + VX(I) ** 2 + VY(I) ** 2 + VZ(I) ** 2                         
                                                                                
1000    CONTINUE                                                                
                                                                                
        K = 0.5 * ( S ** 2 ) * K                                                
                                                                                
        RETURN                                                                  
        END                                                                     
                                                                                
                                                                                
                                                                                
                                                                                
        SUBROUTINE READCN ( CNFILE, BOX )                                       
                                                                                
        COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ, AX, AY, AZ,                   
     :                    BX, BY, BZ, CX, CY, CZ, FX, FY, FZ                    
        COMMON / BLOCK2 / S, SV, SA, SB, SC, SF                                 
                                                                                
C    *******************************************************************        
C    ** SUBROUTINE TO READ IN INITIAL CONFIGURATION FROM UNIT 10      **        
C    *******************************************************************        
                                                                                
        INTEGER     N                                                           
        PARAMETER ( N = 108 )                                                   
                                                                                
        CHARACTER   CNFILE*(*)                                                  
        REAL        BOX                                                         
        REAL        RX(N), RY(N), RZ(N), VX(N), VY(N), VZ(N)                    
        REAL        AX(N), AY(N), AZ(N), BX(N), BY(N), BZ(N)                    
        REAL        CX(N), CY(N), CZ(N), FX(N), FY(N), FZ(N)                    
        REAL        S, SV, SA, SB, SC, SF                                       
                                                                                
        INTEGER     CNUNIT, NN                                                  
        PARAMETER ( CNUNIT = 10 )                                               
                                                                                
C     ******************************************************************        
                                                                                
        OPEN ( UNIT = CNUNIT, FILE = CNFILE,                                    
     :         STATUS= 'OLD', FORM = 'UNFORMATTED' )                            
                                                                                
        READ  ( CNUNIT ) NN, BOX, S, SV, SA, SB, SC                             
        IF ( NN .NE. N ) STOP ' INCORRECT NUMBER OF ATOMS '                     
        READ  ( CNUNIT ) RX, RY, RZ                                             
        READ  ( CNUNIT ) VX, VY, VZ                                             
        READ  ( CNUNIT ) AX, AY, AZ                                             
        READ  ( CNUNIT ) BX, BY, BZ                                             
        READ  ( CNUNIT ) CX, CY, CZ                                             
                                                                                
        CLOSE ( UNIT = CNUNIT )                                                 
                                                                                
        RETURN                                                                  
        END                                                                     
                                                                                
                                                                                
                                                                                
        SUBROUTINE WRITCN ( CNFILE, BOX )                                       
                                                                                
        COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ, AX, AY, AZ,                   
     :                    BX, BY, BZ, CX, CY, CZ, FX, FY, FZ                    
        COMMON / BLOCK2 / S, SV, SA, SB, SC, SF                                 
                                                                                
C    *******************************************************************        
C    ** ROUTINE TO WRITE OUT FINAL CONFIGURATION                      **        
C    *******************************************************************        
                                                                                
        INTEGER     N                                                           
        PARAMETER ( N = 108 )                                                   
                                                                                
        CHARACTER   CNFILE*(*)                                                  
        REAL        BOX                                                         
        REAL        RX(N), RY(N), RZ(N), VX(N), VY(N), VZ(N)                    
        REAL        AX(N), AY(N), AZ(N), BX(N), BY(N), BZ(N)                    
        REAL        CX(N), CY(N), CZ(N), FX(N), FY(N), FZ(N)                    
        REAL        S, SV, SA, SB, SC, SF                                       
                                                                                
        INTEGER     CNUNIT                                                      
        PARAMETER ( CNUNIT = 10 )                                               
                                                                                
C    ****************************************************************           
                                                                                
        OPEN ( UNIT = CNUNIT, FILE = CNFILE,                                    
     :         STATUS='OLD', FORM = 'UNFORMATTED' )                             
                                                                                
        WRITE ( CNUNIT ) N, BOX, S, SV, SA, SB, SC                              
        WRITE ( CNUNIT ) RX, RY, RZ                                             
        WRITE ( CNUNIT ) VX, VY, VZ                                             
        WRITE ( CNUNIT ) AX, AY, AZ                                             
        WRITE ( CNUNIT ) BX, BY, BZ                                             
        WRITE ( CNUNIT ) CX, CY, CZ                                             
                                                                                
        CLOSE ( UNIT = CNUNIT )                                                 
                                                                                
        RETURN                                                                  
        END                                                                     
                                                                                
                                                                                
                                                                                
        SUBROUTINE PREDIC ( DT )                                                
                                                                                
        COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ, AX, AY, AZ,                   
     :                    BX, BY, BZ, CX, CY, CZ, FX, FY, FZ                    
        COMMON / BLOCK2 / S, SV, SA, SB, SC, SF                                 
                                                                                
C    *******************************************************************        
C    ** STANDARD TAYLOR SERIES PREDICTORS                             **        
C    *******************************************************************        
                                                                                
        INTEGER     N                                                           
        PARAMETER ( N = 108 )                                                   
                                                                                
        REAL        RX(N), RY(N), RZ(N), VX(N), VY(N), VZ(N)                    
        REAL        AX(N), AY(N), AZ(N), BX(N), BY(N), BZ(N)                    
        REAL        CX(N), CY(N), CZ(N), FX(N), FY(N), FZ(N)                    
        REAL        S, SV, SA, SB, SC, SF                                       
        REAL        DT                                                          
                                                                                
        INTEGER     I                                                           
        REAL        C1, C2, C3, C4                                              
                                                                                
C    *******************************************************************        
                                                                                
        C1 = DT                                                                 
        C2 = C1 * DT / 2.0                                                      
        C3 = C2 * DT / 3.0                                                      
        C4 = C3 * DT / 4.0                                                      
                                                                                
        DO 100 I = 1, N                                                         
                                                                                
           RX(I) = RX(I) + C1*VX(I) + C2*AX(I) + C3*BX(I) + C4*CX(I)            
           RY(I) = RY(I) + C1*VY(I) + C2*AY(I) + C3*BY(I) + C4*CY(I)            
           RZ(I) = RZ(I) + C1*VZ(I) + C2*AZ(I) + C3*BZ(I) + C4*CZ(I)            
           VX(I) = VX(I) + C1*AX(I) + C2*BX(I) + C3*CX(I)                       
           VY(I) = VY(I) + C1*AY(I) + C2*BY(I) + C3*CY(I)                       
           VZ(I) = VZ(I) + C1*AZ(I) + C2*BZ(I) + C3*CZ(I)                       
           AX(I) = AX(I) + C1*BX(I) + C2*CX(I)                                  
           AY(I) = AY(I) + C1*BY(I) + C2*CY(I)                                  
           AZ(I) = AZ(I) + C1*BZ(I) + C2*CZ(I)                                  
           BX(I) = BX(I) + C1*CX(I)                                             
           BY(I) = BY(I) + C1*CY(I)                                             
           BZ(I) = BZ(I) + C1*CZ(I)                                             
                                                                                
100     CONTINUE                                                                
                                                                                
        S  = S  + C1*SV + C2*SA + C3*SB + C4*SC                                 
        SV = SV + C1*SA + C2*SB + C3*SC                                         
        SA = SA + C1*SB + C2*SC                                                 
        SB = SB + C1*SC                                                         
                                                                                
        RETURN                                                                  
        END                                                                     
                                                                                
                                                                                
                                                                                
        SUBROUTINE CORREC ( DT )                                                
                                                                                
        COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ, AX, AY, AZ,                   
     :                    BX, BY, BZ, CX, CY, CZ, FX, FY, FZ                    
        COMMON / BLOCK2 / S, SV, SA, SB, SC, SF                                 
                                                                                
C    *******************************************************************        
C    ** GEAR CORRECTOR ALGORITHM                                      **        
C    **                                                               **        
C    ** WE USE GEAR'S ALTERNATIVE COEFFICIENT GEAR0 WHICH APPLIES IN  **        
C    ** THE CASE THAT FIRST TIME DERIVATIVES APPEAR ON THE RIGHT OF   **        
C    ** THESE SECOND-ORDER DIFFERENTIAL EQUATIONS.                    **        
C    ** FOR TIMESTEP-SCALED VARIABLES THE COEFFICIENTS WOULD BE       **        
C    ** 19/90, 3/4, 1, 1/2, 1/12.                                     **        
C    **                                                               **        
C    ** REFERENCES:                                                   **        
C    **                                                               **        
C    ** GEAR, NUMERICAL INITIAL VALUE PROBLEMS IN ORDINARY            **        
C    ** DIFFERENTIAL EQUATIONS (PRENTICE-HALL, 1971).                 **        
C    ** GEAR, REPORT ANL 7126 (ARGONNE NATIONAL LAB., 1966).          **        
C    *******************************************************************        
                                                                                
        INTEGER     N                                                           
        PARAMETER ( N = 108 )                                                   
                                                                                
        REAL        RX(N), RY(N), RZ(N), VX(N), VY(N), VZ(N)                    
        REAL        AX(N), AY(N), AZ(N), BX(N), BY(N), BZ(N)                    
        REAL        CX(N), CY(N), CZ(N), FX(N), FY(N), FZ(N)                    
        REAL        S, SV, SA, SB, SC, SF                                       
        REAL        DT                                                          
                                                                                
        INTEGER     I                                                           
        REAL        C1, C2, C3, C4                                              
        REAL        CR, CV, CB, CC                                              
        REAL        CORR, CORRX, CORRY, CORRZ                                   
        REAL        AXI, AYI, AZI                                               
                                                                                
        REAL        GEAR0, GEAR1, GEAR3, GEAR4                                  
        PARAMETER ( GEAR0 = 19.0 / 90.0, GEAR1 = 3.0 /4.0,                      
     :              GEAR3 = 1.0 / 2.0,   GEAR4 = 1.0 /12.0 )                    
                                                                                
C    *******************************************************************        
                                                                                
        C1 = DT                                                                 
        C2 = C1 * DT / 2.0                                                      
        C3 = C2 * DT / 3.0                                                      
        C4 = C3 * DT / 4.0                                                      
                                                                                
        CR = GEAR0 * C2                                                         
        CV = GEAR1 * C2 / C1                                                    
        CB = GEAR3 * C2 / C3                                                    
        CC = GEAR4 * C2 / C4                                                    
                                                                                
        DO 400 I = 1, N                                                         
                                                                                
           AXI = FX(I) / ( S ** 2 ) - 2.0 * SV * VX(I) / S                      
           AYI = FY(I) / ( S ** 2 ) - 2.0 * SV * VY(I) / S                      
           AZI = FZ(I) / ( S ** 2 ) - 2.0 * SV * VZ(I) / S                      
           CORRX = AXI - AX(I)                                                  
           CORRY = AYI - AY(I)                                                  
           CORRZ = AZI - AZ(I)                                                  
                                                                                
           RX(I) = RX(I) + CR * CORRX                                           
           RY(I) = RY(I) + CR * CORRY                                           
           RZ(I) = RZ(I) + CR * CORRZ                                           
           VX(I) = VX(I) + CV * CORRX                                           
           VY(I) = VY(I) + CV * CORRY                                           
           VZ(I) = VZ(I) + CV * CORRZ                                           
           AX(I) = AXI                                                          
           AY(I) = AYI                                                          
           AZ(I) = AZI                                                          
           BX(I) = BX(I) + CB * CORRX                                           
           BY(I) = BY(I) + CB * CORRY                                           
           BZ(I) = BZ(I) + CB * CORRZ                                           
           CX(I) = CX(I) + CC * CORRX                                           
           CY(I) = CY(I) + CC * CORRY                                           
           CZ(I) = CZ(I) + CC * CORRZ                                           
                                                                                
400     CONTINUE                                                                
                                                                                
        CORR = SF - SA                                                          
        S  = S  + CR * CORR                                                     
        SV = SV + CV * CORR                                                     
        SA = SF                                                                 
        SB = SB + CB * CORR                                                     
        SC = SC + CC * CORR                                                     
                                                                                
        RETURN                                                                  
        END                                                                     
