Forces

The forces routine, which is the most CPU intensive routine, was originally written for arbitrary dimension. It was found that the routine could not be vectorized properly without coding it explicitly in three dimensions. This is the only part of the code that is dimension dependent, apart form the initialization routine SPM_quotFCC". It is also desirable to perform potential and pressure calculations at this stage.

C Compute the forces on all the particles due to all the other particles        
C as a function of the position in phase space GAMMA,                           
C and store the result in an array F(i)                                         
                                                                                
       SUBROUTINE FORCES( X )                                                   
       INCLUDE (PARAM)                                                          
                                                                                
C R = displacement between two particles (vector quantity)                      
C RIJSQ = |R|**2                                                                
C FIJ = SCALF(RIJSQ) = scalar part of force term                                
C FKIJ = Kth component of the force between two particles                       
                                                                                
C Three dimensional version                                                     
                                                                                
       DOUBLE PRECISION F(NP*D), X(DIMPS)                                       
       DOUBLE PRECISION FIJ, FKIJ, SCALF, RIJSQ                                 
       DOUBLE PRECISION R1,R2,R3, RSI, R6, V                                    
C F$ in front of the index variables are used to get PI to interact             
C correctly with *VOCL                                                          
       INTEGER F$I,F$J                                                          
       COMMON /FORCE/ F, V                                                      
                                                                                
       DOUBLE PRECISION STEP,Z                                                  
       STEP(Z) = 0.5 + SIGN(0.5D0,Z)                                            
                                                                                
       DO 6 F$I=1,D*NP                                                          
         F(F$I)=0                                                               
6        CONTINUE                                                               
                                                                                
       V = 0                                                                    
                                                                                
*VOCL LOOP,TEMP(R1,R2,R3,RIJSQ,FIJ,FKIJ)                                        
*VOCL LOOP,VI(F$I)                                                              
*VOCL LOOP,NOVREC(F)                                                            
*VOCL LOOP,F$I.LE.NP                                                            
*VOCL LOOP,F$J.LT.NP                                                            
*VOCL LOOP,F$J.LT.F$I      
                                                     
       DO 1 F$I=2,NP                                                            
        DO 2 F$J=1,F$I-1                                                       
         R1 = X(F$I) - X(F$J)                                                   
         R2 = X(F$I+NP) - X(F$J+NP)                                             
         R3 = X(F$I+2*NP) - X(F$J+2*NP)         
                                         
C Perform image particle correction
*VOCL STMT,IF(60)                                                               
         IF (ABS(R1).GT.CUBSIZ/2) R1 = R1 - SIGN(CUBSIZ,R1)                     
*VOCL STMT,IF(60)                                                               
         IF (ABS(R2).GT.CUBSIZ/2) R2 = R2 - SIGN(CUBSIZ,R2)                     
*VOCL STMT,IF(60)                                                               
         IF (ABS(R3).GT.CUBSIZ/2) R3 = R3 - SIGN(CUBSIZ,R3)                     

         RIJSQ = R1 ** 2 + R2 ** 2 + R3 ** 2                                    
                                                                                
         RSI = 1/RIJSQ                                                          
         R6 = RSI ** 3                                                          
         V = V + 4 * (R6 * (R6-1) - RCUT ** (-6) * (RCUT**(-6)-1))              
     &    * STEP(RCUT**2-RIJSQ)                                                 
         FIJ = 4 * 6 * RSI * R6 *(2*R6-1) * STEP(RCUT**2-RIJSQ)                 
                                                                                
         FKIJ = R1 * FIJ                                                        
         F(F$I) = F(F$I) + FKIJ                                                 
         F(F$J) = F(F$J) - FKIJ                                                 
         FKIJ = R2 * FIJ                                                        
         F(F$I+NP) = F(F$I+NP) + FKIJ                                           
         F(F$J+NP) = F(F$J+NP) - FKIJ                                           
         FKIJ = R3 * FIJ                                                        
         F(F$I+2*NP) = F(F$I+2*NP) + FKIJ                                       
         F(F$J+2*NP) = F(F$J+2*NP) - FKIJ                                       
2        CONTINUE                                                               
1       CONTINUE                                                                 
       RETURN                                                                   
       END