C    *******************************************************************        
C    ** FICHE F.3                                                     **        
C    ** TWO SEPARATE PARTS: FORTRAN AND BASIC VERSIONS.               **        
C    *******************************************************************        
                                                                                
                                                                                
                                                                                
C    *******************************************************************        
C    ** FICHE F.3 - PART A                                            **        
C    ** FORTRAN PROGRAM USING THE LEAPFROG ALGORITHM.                 **        
C    *******************************************************************        
                                                                                
                                                                                
                                                                                
        PROGRAM FROGGY                                                          
                                                                                
        COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ                                
                                                                                
C    *******************************************************************        
C    ** FORTRAN PROGRAM TO CONDUCT MOLECULAR DYNAMICS OF ATOMS.       **        
C    **                                                               **        
C    ** A SPECIAL LOW-STORAGE FORM OF THE LEAPFROG VERLET ALGORITHM   **        
C    ** IS USED AND WE TAKE THE LENNARD-JONES POTENTIAL AS AN EXAMPLE.**        
C    **                                                               **        
C    ** REFERENCE:                                                    **        
C    **                                                               **        
C    ** FINCHAM AND HEYES, CCP5 QUARTERLY, 6, 4, 1982.                **        
C    **                                                               **        
C    ** ROUTINES REFERENCED:                                          **        
C    **                                                               **        
C    ** SUBROUTINE READCN ( CNFILE )                                  **        
C    **    READS IN CONFIGURATION                                     **        
C    ** SUBROUTINE FORCE ( DT, SIGMA, RCUT, NEWV, NEWVC, NEWW )       **        
C    **    CALCULATES THE ACCELERATIONS AND ADVANCES THE VELOCITIES   **        
C    **    FROM T - DT/2 TO T + DT/2. ALSO CALCULATES POTENTIAL       **        
C    **    ENERGY AND VIRIAL AT TIMESTEP T.                           **        
C    ** SUBROUTINE MOVE ( DT )                                        **        
C    **    ADVANCES POSITIONS FROM T TO T + DT.                       **        
C    ** SUBROUTINE KINET ( NEWK )                                     **        
C    **    CALCULATES KINETIC ENERGY.                                 **        
C    ** SUBROUTINE WRITCN ( CNFILE )                                  **        
C    **    WRITES OUT CONFIGURATION                                   **        
C    **                                                               **        
C    ** PRINCIPAL VARIABLES:                                          **        
C    **                                                               **        
C    ** INTEGER N                 NUMBER OF ATOMS                     **        
C    ** REAL    DT                TIMESTEP                            **        
C    ** REAL    RX(N),RY(N),RZ(N) ATOMIC POSITIONS                    **        
C    ** REAL    VX(N),VY(N),VZ(N) ATOMIC VELOCITIES                   **        
C    ** REAL    ACV,ACK ETC.      AVERAGE VALUE ACCUMULATORS          **        
C    ** REAL    AVV,AVK ETC.      AVERAGE VALUES                      **        
C    ** REAL    ACVSQ, ACKSQ ETC. AVERAGE SQUARED VALUE ACCUMULATORS  **        
C    ** REAL    FLV,FLK ETC.      FLUCTUATION AVERAGES                **        
C    **                                                               **        
C    ** USAGE:                                                        **        
C    **                                                               **        
C    ** THE LEAPFROG ALGORITHM IS OF THE FORM                         **        
C    ** VX(T + DT/2) = VX(T - DT/2) + DT * AX(T)  (SIMILARLY FOR Y,Z) **        
C    ** RX(T + DT)   = RX(T) + DT * VX(T + DT/2)  (SIMILARLY FOR Y,Z) **        
C    ** TO SAVE STORAGE IN THIS PROGRAM WE ACCUMULATE VALUES AX(T)    **        
C    ** DIRECTLY ONTO THE VELOCITIES VX(T - DT/2) IN THE FORCE LOOP.  **        
C    ** THIS MEANS THAT AN APPROXIMATION MUST BE USED FOR THE KINETIC **        
C    ** ENERGY AT EACH TIME STEP:                                     **        
C    ** K = ( OLDK + NEWK ) / 2 + ( OLDV - 2 * V + NEWV ) / 8         **        
C    ** WHERE K    = K( T        ), V    = V( T )                     **        
C    **       OLDK = K( T - DT/2 ), OLDV = V( T - DT )                **        
C    **       NEWK = K( T + DT/2 ), NEWV = V( T + DT )                **        
C    ** AT THE START OF A STEP THE FOLLOWING VARIABLES ARE STORED:    **        
C    ** R    : R(STEP)         V    : V(STEP+1/2)                     **        
C    ** OLDV : V(STEP-2)       V    : V(STEP-1)        NEWV : V(STEP) **        
C    ** OLDK : K(STEP-1/2)     NEWK : K(STEP+1/2)                     **        
C    ** THIS PROGRAM USES UNIT 10 FOR CONFIGURATION INPUT AND OUTPUT  **        
C    **                                                               **        
C    ** UNITS:                                                        **        
C    **                                                               **        
C    ** THE PROGRAM USES LENNARD-JONES REDUCED UNITS FOR USER INPUT   **        
C    ** AND OUTPUT BUT CONDUCTS SIMULATION IN A BOX OF UNIT LENGTH.   **        
C    ** SUMMARY FOR BOX LENGTH L, ATOMIC MASS M, AND LENNARD-JONES    **        
C    ** POTENTIAL PARAMETERS SIGMA AND EPSILON:                       **        
C    **                OUR PROGRAM           LENNARD-JONES SYSTEM     **        
C    ** LENGTH         L                     SIGMA                    **        
C    ** MASS           M                     M                        **        
C    ** ENERGY         EPSILON               EPSILON                  **        
C    ** TIME           SQRT(M*L**2/EPSILON)  SQRT(M*SIGMA**2/EPSILON) **        
C    ** VELOCITY       SQRT(EPSILON/M)       SQRT(EPSILON/M)          **        
C    ** PRESSURE       EPSILON/L**3          EPSILON/SIGMA**3         **        
C    *******************************************************************        
                                                                                
        INTEGER     N                                                           
        PARAMETER ( N = 108 )                                                   
                                                                                
        REAL        RX(N), RY(N), RZ(N)                                         
        REAL        VX(N), VY(N), VZ(N)                                         
        INTEGER     STEP, NSTEP, IPRINT                                         
        REAL        DENS, NORM, RCUT, DT, SIGMA                                 
        REAL        V, K, E, W                                                  
        REAL        OLDK, NEWK, OLDV, NEWV, NEWW                                
        REAL        OLDVC, NEWVC, VC, KC, EC                                    
        REAL        VN, KN, EN, ECN, PRES, TEMP                                 
        REAL        ACV, ACK, ACE, ACEC, ACP, ACT                               
        REAL        AVV, AVK, AVE, AVEC, AVP, AVT                               
        REAL        ACVSQ, ACKSQ, ACESQ, ACECSQ, ACPSQ, ACTSQ                   
        REAL        FLV, FLK, FLE, FLEC, FLP, FLT                               
        REAL        SR3, SR9, VLRC, WLRC, PI, SIGCUB                            
        CHARACTER   CNFILE*30, TITLE*80                                         
        REAL        FREE                                                        
                                                                                
        PARAMETER ( FREE = 3.0 )                                                
        PARAMETER ( PI = 3.1415927 )                                            
                                                                                
C    *******************************************************************        
                                                                                
C    ** READ IN INITIAL PARAMETERS **                                           
                                                                                
        WRITE(*,'(1H1,'' **** PROGRAM FROGGY ****                  '')')        
        WRITE(*,'(//, '' MOLECULAR DYNAMICS OF LENNARD-JONES ATOMS '')')        
        WRITE(*,'(    '' LEAPFROG ALGORITHM WITH MINIMAL STORAGE   '')')        
                                                                                
        WRITE(*,'('' ENTER RUN TITLE                               '')')        
        READ (*,'(A)') TITLE                                                    
        WRITE(*,'('' ENTER NUMBER OF STEPS                         '')')        
        READ (*,*) NSTEP                                                        
        WRITE(*,'('' ENTER NUMBER OF STEPS BETWEEN OUTPUT LINES    '')')        
        READ (*,*) IPRINT                                                       
        WRITE(*,'('' ENTER CONFIGURATION FILENAME                  '')')        
        READ (*,'(A)') CNFILE                                                   
        WRITE(*,'('' ENTER THE FOLLOWING IN LENNARD-JONES UNITS    '')')        
        WRITE(*,'('' ENTER DENSITY                                 '')')        
        READ (*,*) DENS                                                         
        WRITE(*,'('' ENTER POTENTIAL CUTOFF DISTANCE               '')')        
        READ (*,*) RCUT                                                         
        WRITE(*,'('' ENTER TIME STEP                               '')')        
        READ (*,*) DT                                                           
                                                                                
        WRITE(*,'(//1X,A)') TITLE                                               
        WRITE(*,'('' NUMBER OF ATOMS  = '',I10  )') N                           
        WRITE(*,'('' NUMBER OF STEPS  = '',I10  )') NSTEP                       
        WRITE(*,'('' OUTPUT FREQUENCY = '',I10  )') IPRINT                      
        WRITE(*,'('' POTENTIAL CUTOFF = '',F10.4)') RCUT                        
        WRITE(*,'('' DENSITY          = '',F10.4)') DENS                        
        WRITE(*,'('' TIME STEP        = '',F10.6)') DT                          
                                                                                
C    ** READ CONFIGURATION INTO COMMON / BLOCK1 / VARIABLES **                  
                                                                                
        CALL READCN ( CNFILE )                                                  
                                                                                
C    ** CONVERT INPUT DATA TO PROGRAM UNITS **                                  
                                                                                
        SIGMA = ( DENS / REAL ( N ) ) ** ( 1.0 / 3.0 )                          
        RCUT  = RCUT * SIGMA                                                    
        DT    = DT * SIGMA                                                      
        DENS  = DENS / ( SIGMA ** 3 )                                           
                                                                                
C    ** CALCULATE LONG-RANGE CORRECTIONS **                                     
C    ** NOTE: SPECIFIC TO LENNARD-JONES  **                                     
                                                                                
        SR3    = ( SIGMA / RCUT ) ** 3                                          
        SR9    = SR3 ** 3                                                       
        SIGCUB = SIGMA ** 3                                                     
        VLRC = ( 8.0 /9.0 ) * PI * DENS * SIGCUB * REAL ( N )                   
     :           * ( SR9 - 3.0 * SR3 )                                          
        WLRC = ( 16.0 / 9.0 ) * PI * DENS * SIGCUB * 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                                                              
                                                                                
        ACVSQ  = 0.0                                                            
        ACKSQ  = 0.0                                                            
        ACESQ  = 0.0                                                            
        ACECSQ = 0.0                                                            
        ACPSQ  = 0.0                                                            
        ACTSQ  = 0.0                                                            
                                                                                
        FLV  = 0.0                                                              
        FLK  = 0.0                                                              
        FLE  = 0.0                                                              
        FLEC = 0.0                                                              
        FLP  = 0.0                                                              
        FLT  = 0.0                                                              
                                                                                
C    ** CALCULATE INITIAL VALUES **                                             
                                                                                
        CALL FORCE ( -DT, SIGMA, RCUT, NEWV, NEWVC, NEWW )                      
        CALL MOVE  ( -DT )                                                      
        CALL FORCE ( -DT, SIGMA, RCUT, V, VC, W )                               
        CALL FORCE (  DT, SIGMA, RCUT, V, VC, W )                               
        CALL KINET ( OLDK )                                                     
        CALL MOVE  (  DT )                                                      
        CALL FORCE (  DT, SIGMA, RCUT, NEWV, NEWVC, NEWW )                      
        CALL KINET ( NEWK )                                                     
                                                                                
C    ** INCLUDE LONG-RANGE CORRECTIONS **                                       
                                                                                
        V    = V + VLRC                                                         
        W    = W + WLRC                                                         
        NEWV = NEWV + VLRC                                                      
        NEWW = NEWW + WLRC                                                      
                                                                                
        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 MOVE ( DT )                                                     
                                                                                
           OLDV  = V                                                            
           V     = NEWV                                                         
           OLDVC = VC                                                           
           VC    = NEWVC                                                        
           W     = NEWW                                                         
                                                                                
           CALL FORCE ( DT, SIGMA, RCUT, NEWV, NEWVC, NEWW )                    
                                                                                
C       ** INCLUDE LONG-RANGE CORRECTIONS **                                    
                                                                                
           NEWV = NEWV + VLRC                                                   
           NEWW = NEWW + WLRC                                                   
                                                                                
C       ** CALCULATE KINETIC ENERGY AT CURRENT STEP **                          
                                                                                
           K  = ( NEWK + OLDK ) / 2.0                                           
     :        + ( NEWV  - 2.0 * V  + OLDV  ) / 8.0                              
           KC = ( NEWK + OLDK ) / 2.0                                           
     :        + ( NEWVC - 2.0 * VC + OLDVC ) / 8.0                              
           OLDK = NEWK                                                          
                                                                                
           CALL KINET ( NEWK )                                                  
                                                                                
C       ** CALCULATE INSTANTANEOUS VALUES **                                    
                                                                                
           E    = K + V                                                         
           EC   = KC + VC                                                       
           VN   = V  / REAL ( N )                                               
           KN   = K  / REAL ( N )                                               
           EN   = E  / REAL ( N )                                               
           ECN  = EC / REAL ( N )                                               
           TEMP = 2.0 * KN / FREE                                               
           PRES = DENS * TEMP + W                                               
                                                                                
C       ** CONVERT TO LENNARD-JONES UNITS **                                    
                                                                                
           PRES = PRES * SIGMA ** 3                                             
                                                                                
C       ** INCREMENT ACCUMULATORS **                                            
                                                                                
           ACE  = ACE  + EN                                                     
           ACEC = ACEC + ECN                                                    
           ACK  = ACK  + KN                                                     
           ACV  = ACV  + VN                                                     
           ACP  = ACP  + PRES                                                   
                                                                                
           ACESQ  = ACESQ  + EN  ** 2                                           
           ACECSQ = ACECSQ + ECN ** 2                                           
           ACKSQ  = ACKSQ  + KN  ** 2                                           
           ACVSQ  = ACVSQ  + VN  ** 2                                           
           ACPSQ  = ACPSQ  + PRES ** 2                                          
                                                                                
C       ** OPTIONALLY PRINT INFORMATION **                                      
                                                                                
           IF ( MOD( STEP, IPRINT ) .EQ. 0 ) THEN                               
                                                                                
              WRITE(*,'(1X,I8,6(2X,F10.4))')                                    
     :                 STEP, EN, ECN, KN, VN, PRES, TEMP                        
           ENDIF                                                                
                                                                                
1000    CONTINUE                                                                
                                                                                
C    *******************************************************************        
C    ** MAIN LOOP ENDS                                                **        
C    *******************************************************************        
                                                                                
        WRITE(*,'(/1X,''**** END OF DYNAMICS **** ''//)')                       
                                                                                
C    ** OUTPUT RUN AVERAGES **                                                  
                                                                                
        NORM   = REAL ( NSTEP )                                                 
                                                                                
        AVE  = ACE  / NORM                                                      
        AVEC = ACEC / NORM                                                      
        AVK  = ACK  / NORM                                                      
        AVV  = ACV  / NORM                                                      
        AVP  = ACP  / 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                                  
                                                                                
        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  )                           
                                                                                
        AVT = AVK * 2.0 / FREE                                                  
        FLT = FLK * 2.0 / FREE                                                  
                                                                                
        WRITE(*,'('' AVERAGES'',6(2X,F10.5))')                                  
     :            AVE, AVEC, AVK, AVV, AVP, AVT                                 
        WRITE(*,'('' FLUCTS  '',6(2X,F10.5))')                                  
     :            FLE, FLEC, FLK, FLV, FLP, FLT                                 
                                                                                
C    ** WRITE OUT FINAL CONFIGURATION **                                        
                                                                                
        CALL WRITCN ( CNFILE )                                                  
                                                                                
        STOP                                                                    
                                                                                
10001   FORMAT(//1X,'TIMESTEP  ..ENERGY..  CUTENERGY.'                          
     :              '  ..KINETIC.  ..POTENT..',                                 
     :              '  .PRESSURE.  ..TEMPER..'/)                                
        END                                                                     
                                                                                
                                                                                
                                                                                
        SUBROUTINE FORCE ( DT, SIGMA, RCUT, V, VC, W )                          
                                                                                
        COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ                                
                                                                                
C    *******************************************************************        
C    ** ROUTINE TO COMPUTE FORCES WITH CUTOFF & MINIMUM IMAGING.      **        
C    **                                                               **        
C    ** THE ROUTINE ACTUALLY RETURNS TWO DIFFERENT POTENTIAL ENERGIES.**        
C    ** V IS CALCULATED USING THE UNSHIFTED LENNARD-JONES POTENTIAL.  **        
C    ** WHEN LONG-RANGE TAIL CORRECTIONS ARE ADDED, THIS MAY BE USED  **        
C    ** TO CALCULATE THERMODYNAMIC INTERNAL ENERGY ETC.               **        
C    ** VC IS CALCULATED USING THE SHIFTED LENNARD-JONES POTENTIAL    **        
C    ** WITH NO DISCONTINUITY AT THE CUTOFF.  THIS MAY BE USED TO     **        
C    ** CHECK ENERGY CONSERVATION.                                    **        
C    ** PROGRAM UNITS MAKE EPSILON = MASS = 1, BUT SIGMA IS NOT UNITY **        
C    ** SINCE WE TAKE A BOX OF UNIT LENGTH.                           **        
C    ** THE ROUTINE IS COMBINED WITH THE LEAPFROG ALGORITHM FOR       **        
C    ** ADVANCING VELOCITIES, I.E. FORCES ARE ACCUMULATED DIRECTLY    **        
C    ** ONTO VELOCITIES.                                              **        
C    *******************************************************************        
                                                                                
        INTEGER     N                                                           
        PARAMETER ( N = 108 )                                                   
        REAL        RX(N), RY(N), RZ(N)                                         
        REAL        VX(N), VY(N), VZ(N)                                         
        REAL        V, RCUT, DT, SIGMA, W, VC                                   
                                                                                
        INTEGER     I, J, NCUT                                                  
        REAL        RCUTSQ, VCUT, SIGSQ                                         
        REAL        RXI, RYI, RZI, VXI, VYI, VZI                                
        REAL        RXIJ, RYIJ, RZIJ, RIJSQ                                     
        REAL        SR2, SR6, SR12                                              
        REAL        VIJ, WIJ, VELIJ, DVX, DVY, DVZ                              
                                                                                
C    *******************************************************************        
                                                                                
C    ** CALCULATE SQUARED DISTANCES FOR INNER LOOP **                           
                                                                                
        SIGSQ  = SIGMA ** 2                                                     
        RCUTSQ = RCUT ** 2                                                      
                                                                                
C    ** CONDUCT OUTER LOOP OVER ATOMS **                                        
                                                                                
        V    = 0.0                                                              
        W    = 0.0                                                              
        NCUT = 0                                                                
                                                                                
        DO 1000 I = 1, N - 1                                                    
                                                                                
           RXI = RX(I)                                                          
           RYI = RY(I)                                                          
           RZI = RZ(I)                                                          
           VXI = VX(I)                                                          
           VYI = VY(I)                                                          
           VZI = VZ(I)                                                          
                                                                                
C       ** CONDUCT INNER LOOP OVER ATOMS **                                     
                                                                                
           DO 900 J = I + 1, N                                                  
                                                                                
              RXIJ = RXI - RX(J)                                                
              RXIJ = RXIJ - ANINT ( RXIJ )                                      
                                                                                
              IF ( ABS ( RXIJ ) .LT. RCUT ) THEN                                
                                                                                
                 RYIJ = RYI - RY(J)                                             
                 RYIJ = RYIJ - ANINT ( RYIJ )                                   
                 RIJSQ = RXIJ ** 2 + RYIJ ** 2                                  
                                                                                
                 IF ( RIJSQ .LT. RCUTSQ ) THEN                                  
                                                                                
                    RZIJ = RZI - RZ(J)                                          
                    RZIJ = RZIJ - ANINT ( RZIJ )                                
                    RIJSQ = RIJSQ + RZIJ ** 2                                   
                                                                                
                    IF ( RIJSQ .LT. RCUTSQ ) THEN                               
                                                                                
                       SR2   = SIGSQ / RIJSQ                                    
                       SR6   = SR2 ** 3                                         
                       SR12  = SR6 ** 2                                         
                       VIJ   = 4.0 * ( SR12 - SR6 )                             
                       WIJ   = 24.0 * ( 2.0 * SR12 - SR6 )                      
                       VELIJ = WIJ * DT / RIJSQ                                 
                       DVX   = VELIJ * RXIJ                                     
                       DVY   = VELIJ * RYIJ                                     
                       DVZ   = VELIJ * RZIJ                                     
                       VXI   = VXI + DVX                                        
                       VYI   = VYI + DVY                                        
                       VZI   = VZI + DVZ                                        
                       VX(J) = VX(J) - DVX                                      
                       VY(J) = VY(J) - DVY                                      
                       VZ(J) = VZ(J) - DVZ                                      
                       V     = V + VIJ                                          
                       W     = W + WIJ                                          
                       NCUT  = NCUT + 1                                         
                                                                                
                    ENDIF                                                       
                                                                                
                 ENDIF                                                          
                                                                                
              ENDIF                                                             
                                                                                
900        CONTINUE                                                             
                                                                                
C       ** END OF INNER LOOP **                                                 
                                                                                
           VX(I) = VXI                                                          
           VY(I) = VYI                                                          
           VZ(I) = VZI                                                          
                                                                                
1000    CONTINUE                                                                
                                                                                
C    ** END OF OUTER LOOP **                                                    
                                                                                
C    ** CALCULATE POTENTIAL SHIFT **                                            
                                                                                
        SR2  = SIGSQ / RCUTSQ                                                   
        SR6  = SR2 * SR2 * SR2                                                  
        SR12 = SR6 * SR6                                                        
        VIJ  = 4.0 * ( SR12 - SR6 )                                             
        VC   = V - REAL ( NCUT ) * VIJ                                          
                                                                                
        W    = W / 3.0                                                          
                                                                                
        RETURN                                                                  
        END                                                                     
                                                                                
                                                                                
                                                                                
        SUBROUTINE READCN ( CNFILE )                                            
                                                                                
        COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ                                
                                                                                
C    *******************************************************************        
C    ** SUBROUTINE TO READ IN INITIAL CONFIGURATION FROM UNIT 10      **        
C    *******************************************************************        
                                                                                
        INTEGER     N                                                           
        PARAMETER ( N = 108 )                                                   
        CHARACTER   CNFILE*(*)                                                  
        REAL        RX(N), RY(N), RZ(N)                                         
        REAL        VX(N), VY(N), VZ(N)                                         
                                                                                
        INTEGER     CNUNIT, NN                                                  
        PARAMETER ( CNUNIT = 10 )                                               
                                                                                
C     ******************************************************************        
                                                                                
        OPEN ( UNIT = CNUNIT, FILE = CNFILE, STATUS = 'OLD',                    
     :         FORM = 'UNFORMATTED' )                                           
                                                                                
        READ ( CNUNIT ) NN                                                      
        IF ( NN .NE. N ) STOP ' INCORRECT NUMBER OF ATOMS '                     
        READ ( CNUNIT ) RX, RY, RZ                                              
        READ ( CNUNIT ) VX, VY, VZ                                              
                                                                                
        CLOSE ( UNIT = CNUNIT )                                                 
                                                                                
        RETURN                                                                  
        END                                                                     
                                                                                
                                                                                
                                                                                
        SUBROUTINE WRITCN ( CNFILE )                                            
                                                                                
        COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ                                
                                                                                
C    *******************************************************************        
C    ** ROUTINE TO WRITE OUT FINAL CONFIGURATION TO UNIT 10           **        
C    *******************************************************************        
                                                                                
        INTEGER     N                                                           
        PARAMETER ( N = 108 )                                                   
        CHARACTER   CNFILE*(*)                                                  
        REAL        RX(N), RY(N), RZ(N)                                         
        REAL        VX(N), VY(N), VZ(N)                                         
                                                                                
        INTEGER     CNUNIT                                                      
        PARAMETER ( CNUNIT = 10 )                                               
                                                                                
C    ****************************************************************           
                                                                                
        OPEN  ( UNIT = CNUNIT, FILE = CNFILE, STATUS = 'UNKNOWN',               
     :          FORM = 'UNFORMATTED' )                                          
                                                                                
        WRITE ( CNUNIT ) N                                                      
        WRITE ( CNUNIT ) RX, RY, RZ                                             
        WRITE ( CNUNIT ) VX, VY, VZ                                             
                                                                                
        CLOSE ( UNIT = CNUNIT )                                                 
                                                                                
        RETURN                                                                  
        END                                                                     
                                                                                
                                                                                
                                                                                
        SUBROUTINE MOVE ( DT )                                                  
                                                                                
        COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ                                
                                                                                
C    *******************************************************************        
C    ** LEAPFROG VERLET ALGORITHM FOR ADVANCING POSITIONS             **        
C    *******************************************************************        
                                                                                
        INTEGER     N                                                           
        PARAMETER ( N = 108 )                                                   
        REAL        RX(N), RY(N), RZ(N)                                         
        REAL        VX(N), VY(N), VZ(N)                                         
        REAL        DT                                                          
                                                                                
        INTEGER     I                                                           
                                                                                
C    *******************************************************************        
                                                                                
        DO 1000 I = 1, N                                                        
                                                                                
           RX(I) = RX(I) + VX(I) * DT                                           
           RY(I) = RY(I) + VY(I) * DT                                           
           RZ(I) = RZ(I) + VZ(I) * DT                                           
                                                                                
1000    CONTINUE                                                                
                                                                                
        RETURN                                                                  
        END                                                                     
                                                                                
                                                                                
                                                                                
        SUBROUTINE KINET ( K )                                                  
                                                                                
        COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ                                
                                                                                
C    *******************************************************************        
C    ** COMPUTES KINETIC ENERGY                                       **        
C    *******************************************************************        
                                                                                
        INTEGER     N                                                           
        PARAMETER ( N = 108 )                                                   
        REAL        RX(N), RY(N), RZ(N)                                         
        REAL        VX(N), VY(N), VZ(N)                                         
        REAL        K                                                           
                                                                                
        INTEGER     I                                                           
                                                                                
C    *******************************************************************        
                                                                                
        K  = 0.0                                                                
                                                                                
        DO 100 I = 1, N                                                         
                                                                                
           K = K + VX(I) ** 2 + VY(I) ** 2 + VZ(I) ** 2                         
                                                                                
100     CONTINUE                                                                
                                                                                
        K = K * 0.5                                                             
                                                                                
        RETURN                                                                  
        END                                                                     
                                                                                
                                                                                
                                                                                
C    *******************************************************************        
C    ** FICHE F.3 - PART B                                            **        
C    ** BASIC PROGRAM USING THE LEAPFROG VERLET ALGORITHM.            **        
C    *******************************************************************        
                                                                                
C    *******************************************************************        
C    ** THE FOLLOWING PROGRAMS WERE WRITTEN ON AN ACORN MODEL B MICRO,**        
C    ** USING BBC BASIC AND SOME MACHINE CODE FOR THE GRAPHICS.       **        
C    ** CONSEQUENTLY MINOR ALTERATIONS MAY BE NEEDED FOR OTHER MICROS.**        
C    ** THE FIRST PROGRAM RUNS MOLECULAR DYNAMICS FROM AN INITIAL     **        
C    ** CONFIGURATION WHICH MAY BE GENERATED USING THE SECOND PROGRAM.**        
C    ** DEFAULT INPUT PARAMETERS ARE PROVIDED IN THE MD PROGRAM:      **        
C    ** THESE CAN BE ALTERED TO SUIT THE USER.                        **        
C    *******************************************************************        
                                                                                
                                                                                
                                                                                
       MOLECULAR DYNAMICS PROGRAM                                               
                                                                                
  10 IF PAGE<>&1700 THEN PAGE=&1700:CHAIN"2.FROG1"                              
   20 MODE1                                                                     
   30 CLS:PRINT''''TAB(12);"Program Froggy"                                     
   40PRINT'TAB(10);"Molecular Dynamics"'                                        
   50PRINTTAB(7);" of Lennard-Jones atoms"                                      
   60 PRINT'TAB(1);"Leapfrog algorithm with minimal storage"                    
   70 PRINT'TAB(7);"Repulsive WCA potential"'                                   
   80 PRINTTAB(9);"in two dimensions"                                           
   90 PRINT'TAB(15);"Fiche 3b"                                                  
  100 INPUTTAB(10,31)"( press <RETURN> )"A$                                     
  110 CLS                                                                       
  120*KEY0"100"                                                                 
  130*KEY1"2.DATA"                                                              
  140*KEY2"0.3"                                                                 
  150*KEY3"0.02"                                                                
  160 N=10:FREE=2:x1=200:y1=200:x2=512:y2=512                                   
  170 HIMEM=&2E00                                                               
  180DIM X(12),Y(12),RX(N),RY(N),VX(N),VY(N)                                    
  190HIMEM=&2E00:GCOL 3,4:VDU 19,2,1,0,0,0                                      
  200 *FX138,0,128                                                              
  210NSTEP%= FNgetnum( "Enter number of steps")                                 
  220 *FX138,0,129                                                              
  230PRINT"Enter configuration filename ";:INPUTCNFILE$                         
  240PRINT'TAB(5);"IN Lennard-Jones units,"                                     
  250 *FX138,0,130                                                              
  260DENS= FNgetnum("Enter density")                                            
  270 *FX138,0,131                                                              
  280DT= FNgetnum("Enter time step")                                            
  290PROCreadcn(CNFILE$):SIGMA=(DENS/N)^(1/2):RCUT=SIGMA*(2.0^(1.0/6.0))        
  300DT=DT*SIGMA                                                                
  310CLS                                                                        
  320R=(DENS*(x2-x1)*(y2-y1)/N)^0.5/2                                           
  330PROCcircle(R)                                                              
  340PROCcode                                                                   
  350PROCforce( -DT,SIGMA,RCUT )                                                
  360PROCmove( -DT )                                                            
  370PROCforce( -DT,SIGMA,RCUT )                                                
  380PROCforce(  DT,SIGMA,RCUT )                                                
  390PROCmove( DT )                                                             
  400PROCforce(DT,SIGMA,RCUT )                                                  
  410GCOL 3,1                                                                   
  420PROCbox(x1-2*R,y1,x1+x2-2*R,y1+y2)                                         
  430GCOL 3,3                                                                   
  440PROCgraphics(&2F00,FALSE)                                                  
  450PROCerase                                                                  
  460FOR STP%=1 TO NSTEP%                                                       
  470PROCmove( DT )                                                             
  480PROCforce(DT,SIGMA,RCUT )                                                  
  490PROCgraphics(&2F04,TRUE)                                                   
  500NEXT STP%                                                                  
  510END                                                                        
  520DEF FNgetnum(S$)                                                           
  530PRINT S$;" ";:INPUT A                                                      
  540=A                                                                         
  550DEF PROCreadcn( FILE$ )                                                    
  560LOCAL PORT                                                                 
  570PORT=OPENUP( FILE$ )                                                       
  580IF PORT=0 THEN PRINT''CHR$(34);FILE$;CHR$(34);" - File not found"':END     
  590FOR I%= 1 TO N:INPUT#PORT,RX(I%),RY(I%):NEXT I%                            
  600FOR I%= 1 TO N:INPUT#PORT,VX(I%),VY(I%):NEXT I%                            
  610CLOSE #PORT                                                                
  620DEF PROCmove( DT )                                                         
  630FOR I%= 1 TO N                                                             
  640RX(I%)= RX(I%)+VX(I%)*DT                                                   
  650IF RX(I%)>0.5 REPEAT:RX(I%)=RX(I%)-1:UNTIL RX(I%)<0.5                      
  660IF RX(I%)<-0.5 REPEAT:RX(I%)=RX(I%)+1:UNTIL RX(I%)>-0.5                    
  670RY(I%)= RY(I%)+VY(I%)*DT                                                   
  680IF RY(I%)>0.5 REPEAT:RY(I%)=RY(I%)-1:UNTIL RY(I%)<0.5                      
  690IF RY(I%)<-0.5 REPEAT:RY(I%)=RY(I%)+1:UNTIL RY(I%)>-0.5                    
  700NEXT I%                                                                    
  710ENDPROC                                                                    
  720DEF PROCforce( DT,SIGMA,RCUT )                                             
  730A=SIGMA*SIGMA:B=RCUT*RCUT                                                  
  740FOR I%=1 TO N-1                                                            
  750C=RX(I%):D=RY(I%):H=VX(I%):I=VY(I%)                                        
  760FOR J%= I%+1 TO N                                                          
  770E=C-RX(J%):E=E-INT(E+0.5)                                                  
  780IF ABS(E)>=RCUT THEN 830                                                   
  790F=D-RY(J%):F=F-INT(F+0.5):G=E*E+F*F                                        
  800IF G>=B THEN 830                                                           
  810S=A/G:T=S*S*S:U=T*T:W=24*(U+U-S):J=W*DT/G:M=J*E:O=J*F:H=H+M:I=I+O          
  820VX(J%)=VX(J%)-M:VY(J%)=VY(J%)-O                                            
  830NEXT J%                                                                    
  840VX(I%)=H:VY(I%)=I                                                          
  850NEXT I%                                                                    
  860ENDPROC                                                                    
  870DEF PROCcircle(R)                                                          
  880D=100:T=0:X(1)=R*SIN(T):Y(1)=SQR(R^2-X(1)^2)                               
  890FOR A=2 TO 12                                                              
  900X(A)=X(A-1)*COS(D)+Y(A-1)*SIN(D):Y(A)=Y(A-1)*COS(D)-X(A-1)*SIN(D)          
  910NEXT                                                                       
  920B%=&2E00                                                                   
  930FOR A%=1 TO 12                                                             
  940?B%=25:B%?1=1:B%?2=X(A%) MOD 256:B%?4=Y(A%) MOD 256                        
  950IF X(A%)<0 B%?3=255-(X(A%) DIV 256) ELSE B%?3=X(A%) DIV 256                
  960IF Y(A%)<0 B%?5=255-(Y(A%) DIV 256) ELSE B%?5=Y(A%) DIV 256                
  970B%=B%+6                                                                    
  980NEXT                                                                       
  990?&2E17=0                                                                   
 1000ENDPROC                                                                    
 1010DEF PROCcode                                                               
 1020P%=&1100                                                                   
 1030FOR I%=0 TO 2 STEP2                                                        
 1040\OPT I%                                                                    
 1050 .circle LDY#0                                                             
 1060.loop LDA &2E00,Y                                                          
 1070JSR &FFEE                                                                  
 1080INY                                                                        
 1090CPY #72                                                                    
 1100BNE loop                                                                   
 1110RTS                                                                        
 1120.code LDX #20                                                              
 1130 LDA #&00                                                                  
 1140 STA &70                                                                   
 1150 LDA #&2F                                                                  
 1160 STA &71                                                                   
 1170 .l2 JSR getnum                                                            
 1180 JSR pl2                                                                   
 1190 JSR image                                                                 
 1200 INC &70                                                                   
 1210 INC &70                                                                   
 1220 INC &70                                                                   
 1230 INC &70                                                                   
 1240 DEX                                                                       
 1250 BNE l2                                                                    
 1260 RTS                                                                       
 1270.getnum LDY#3                                                              
 1280.getnum1 LDA(&70),Y                                                        
 1290 STA &72,Y                                                                 
 1300 DEY                                                                       
 1310 BPL getnum1                                                               
 1320 RTS                                                                       
 1330.pl2 LDA #25                                                               
 1340 JSR &FFEE                                                                 
 1350 LDA #4                                                                    
 1360 JSR &FFEE                                                                 
 1370 LDY #0                                                                    
 1380.pl21 LDA &72,Y                                                            
 1390 JSR &FFEE                                                                 
 1400 INY                                                                       
 1410 CPY #4                                                                    
 1420 BNE pl21                                                                  
 1430 JSR circle                                                                
 1440 RTS                                                                       
 1450                                                                           
 1460.image                                                                     
 1470 LDA &73                                                                   
 1480 ADC #1                                                                    
 1490 STA &73                                                                   
 1500 CMP #3                                                                    
 1510 BPL P%+5                                                                  
 1520 JSR pl2                                                                   
 1530 JSR getnum                                                                
 1540 .n7 LDA &75                                                               
 1550 ADC #1                                                                    
 1560 STA &75                                                                   
 1570 CMP #3                                                                    
 1580 BPL P%+5                                                                  
 1590 JSR pl2                                                                   
 1600JSR getnum                                                                 
 1610LDA &75                                                                    
 1620SEC                                                                        
 1630SBC #2                                                                     
 1640STA &75                                                                    
 1650CMP #0                                                                     
 1660BNE P%+11                                                                  
 1670 LDA &74                                                                   
 1680 CMP #&80                                                                  
 1690 BMI P%+5                                                                  
 1700JSR pl2                                                                    
 1710JSR getnum                                                                 
 1720LDA &73                                                                    
 1730SEC                                                                        
 1740SBC #2                                                                     
 1750STA &73                                                                    
 1760CMP #0                                                                     
 1770BNE P%+11                                                                  
 1780  LDA &72                                                                  
 1790 CMP #&80                                                                  
 1800 BMI P%+5                                                                  
 1810 JSR pl2                                                                   
 1820 RTS                                                                       
 1830.go JSR code                                                               
 1840 LDY #0                                                                    
 1850 .l3 LDA &2F04,Y                                                           
 1860 STA &2F00,Y                                                               
 1870 INY                                                                       
 1880 CPY #81                                                                   
 1890 BNE l3                                                                    
 1900 RTS                                                                       
 1910|                                                                          
 1920NEXT                                                                       
 1930ENDPROC                                                                    
 1940DEF PROCgraphics(B%,run)                                                   
 1950FOR A%=1 TO 10                                                             
 1960 D%=x1+(RX(A%)+0.5)*x2:C%=y1+(RY(A%)+0.5)*y2:?(B%+(A%-1)*8)=D% MOD 256     
 1970?(B%+(A%-1)*8+1)=D% DIV 256:?(B%+(A%-1)*8+2)=C% MOD 256:                   
 1980?(B%+(A%-1)*8+3)=C% DIV 256                                                
 1990 NEXT                                                                      
 2000 IF run CALL go                                                            
 2010ENDPROC                                                                    
 2020DEF PROCerase                                                              
 2030FOR A%=1 TO 10                                                             
 2040 X%=x1+(RX(A%)+0.5)*x2                                                     
 2050 Y%=y1+(RY(A%)+0.5)*y2                                                     
 2060MOVE x1+(RX(A%)+0.5)*x2,y1+(RY(A%)+0.5)*y2                                 
 2070CALL circle                                                                
 2080 IF X%+512<256*3 MOVEX%+512,Y%:CALL circle                                 
 2090 IF Y%+512<256*3 MOVEX%,Y%+512:CALL circle                                 
 2100 IF (Y%-512<256) AND (Y%-512)>128 MOVEX%,Y%-512:CALL circle                
 2110 IF X%-512<256 AND X%-512>128 MOVEX%-512,Y%:CALL circle                    
 2120NEXT                                                                       
 2130ENDPROC                                                                    
 2140DEF PROCbox(A%,B%,C%,D%)                                                   
 2150MOVE A%,B%:DRAW A%,D%:DRAW C%,D%:DRAW C%,B%:DRAW A%,B%                     
 2160ENDPROC                                                                    
                                                                                
                                                                                
                                                                                
      CONFIGURATION GENERATOR                                                   
                                                                                
   10 MODE 1                                                                    
   20REM PROGRAM SETUP                                                          
   30REM                                                                        
   40N= 10: REM CONSTANT                                                        
   50DIM RX(N),RY(N)                                                            
   60DIM VX(N),VY(N)                                                            
   70 A1= 3.949846138: A3= 0.252408784: A5= 0.076542912                         
   80 A7= 0.008355968: A9= 0.029899776                                          
   90PROCfcc                                                                    
  100 *KEY0"1.75"                                                               
  110 *KEY1"2.DATA"                                                             
  120                                                                           
  130 *FX138,0,128                                                              
  140PRINT'" Enter Temperature in LJ units ";:INPUT TEMP                        
  150PROCcomvel( TEMP )                                                         
  160                                                                           
  170 *FX138,0,129                                                              
  180PRINT'"Enter Output Data File ";:INPUT CNFILE$                             
  190                                                                           
  200PORT= OPENOUT( CNFILE$ )                                                   
  210IF PORT=0:PRINT''"Error in opening ";CHR$(34);CNFILE$;CHR$(34)'':END       
  220                                                                           
  230FOR I%=1 TO N:PRINT #PORT,RX(I%),RY(I%):NEXT I%                            
  240FOR I%=1 TO N:PRINT #PORT,VX(I%),VY(I%):NEXT I%                            
  250CLOSE #PORT                                                                
  260END                                                                        
  270                                                                           
  280 DEF PROCfcc                                                               
  290 LOCAL X,Y,M,NC,KL,KX                                                      
  300 NC=INT((N/4)^(1/2) +0.5)                                                  
  310 X=1/NC:Y=1                                                                
  320RX(1)=0:RY(1)=0                                                            
  330RX(2)=2*Y/6:RY(2)=0                                                        
  340RX(3)=4*Y/6:RY(3)=0                                                        
  350RX(4)=Y/6:RY(4)=Y/5                                                        
  360RX(5)=3*Y/6:RY(5)=Y/5                                                      
  370RX(6)=0:RY(6)=2*Y/5                                                        
  380RX(7)=2*Y/6:RY(7)=2*Y/5                                                    
  390RX(8)=4*Y/6:RY(8)=2*Y/5                                                    
  400RX(9)=Y/6:RY(9)=3*Y/5                                                      
  410RX(10)=3*Y/6:RY(10)=3*Y/5                                                  
  420 FOR A%=1 TO 10:RX(A%)=RX(A%)-0.5:RY(A%)=RY(A%)-0.5:NEXT                   
  430 ENDPROC                                                                   
  440DEF PROCcomvel( TEMP )                                                     
  450                                                                           
  460LOCAL RTEMP,SUMX,SUMY,GAUSS,DUMMY                                          
  470                                                                           
  480RTEMP= SQR( TEMP )                                                         
  490FOR I%= 1 TO N                                                             
  500VX(I%)= RTEMP*FNgauss                                                      
  510VY(I%)= RTEMP*FNgauss                                                      
  520NEXT I%                                                                    
  530SUX= 0:SUMY= 0                                                             
  540FOR I%=1 TO N                                                              
  550SUMX= SUMX+VX(I%)                                                          
  560SUMY= SUMY+VY(I%)                                                          
  570NEXT I%                                                                    
  580SUMX= SUMX/N: SUMY= SUMY/N                                                 
  590FOR I%=1 TO N                                                              
  600VX(I%)= VX(I%)-SUMX                                                        
  610VY(I%)= VY(I%)-SUMY                                                        
  620NEXT I%                                                                    
  630ENDPROC                                                                    
  640DEF FNgauss                                                                
  650LOCAL SUM,R,R2                                                             
  660FOR L%=1 TO 12                                                             
  670SUM= SUM+RND(1)                                                            
  680NEXT L%                                                                    
  690R= (SUM-6)/4: R2= R*R                                                      
  700= ((((A9*R2+A7)*R2+A5)*R2+A3)*R2+A1)*R                                     
  710END                                                                        
