Equations of Motion

  As an example of how equations of motion are coded as a function, here is the function implementing the current-statting dynamics in equation (). 

C Equations of motions for the Norton ensemble, with Nose-Hoover                
C thermostat and current statting mechanisms in the form                        
C (d/dt) GAMMA(I) = NORTON(I).                                                  
                                                                                
       DOUBLE PRECISION FUNCTION NORTON( I,GAMMA)                               
       INCLUDE (PARAM)                                                          
       DOUBLE PRECISION GAMMA(DIMPS+2), F(D*NP)                                 
       DOUBLE PRECISION V, K, CALJ                                              
       COMMON /FORCE/ F,V                                                       
       INTEGER I , ALPHA, LAMBDA, P                                             
C Offsets into GAMMA                                                            
       PARAMETER (ALPHA=DIMPS+1,LAMBDA=DIMPS+2,P=DIMPS/2)                       
       INCLUDE (STFUN)                                                          
                                                                                
*VOCL STMT,IF(0)                                                                
       IF (I.LE.DIMPS/2) THEN                                                   
         NORTON = GAMMA(P+I)                                                    
*VOCL STMT,IF(100)                                                              
       ELSE IF (THERMO) THEN                                                    
*VOCL STMT,IF(33)                                                               
         IF (COLOUR.AND.I.LE.P+NP) THEN                                         
           NORTON = F(I-P) - GAMMA(ALPHA) * GAMMA(I)  -                         
     &         GAMMA(LAMBDA) * E(I)                                             
         ELSE IF (I.EQ.ALPHA) THEN                                              
C equation of motion for alpha                                                  
           NORTON = ( K( GAMMA( P+1 ) ) / (1.5*NP*TEMP) - 1 ) /                 
     &           TAU ** 2                                                       
         ELSE IF (COLOUR.AND.I.EQ.LAMBDA) THEN                                  
C equation of motion for lambda                                                 
           NORTON = NP/QL * (CALJ( GAMMA( P+1 )) - J0)                          
         ELSE                                                                   
            NORTON = F(I-P) - GAMMA(ALPHA) * GAMMA(I)                           
         ENDIF                                                                  
       ELSE                                                                     
         NORTON = F( I - P)                                                     
         ENDIF                                                                  
       RETURN                                                                   
       END

There are two things to note in this listing. The first is that there is generally some preprocessing to be done independently of the component number SPM_quotI". If Fortran allowed vector valued functions, this would be no problem, as we could do this preprocessing , and then return all components of SPM_quotF" at once. Instead, we must use a kludge in which a call to SPM_quotF" with SPM_quotI" set to 0 is performed before the individual components are returned. The second thing to note is that Fortran does not allow the equivalencing of subroutine parameters. Instead, to make thing a little more readable, I have defined the constants SPM_quotALPHA", SPM_quotLAMBDA" and SPM_quotP" as indexes to SPM_quotGAMMA", so for example SPM_quotGAMMA(ALPHA)" is the thermostatting multiplier, and SPM_quotGAMMA(P+1)" is .