C ******************************************************************* C ** FICHE F.17 ** C ** FORCE CALCULATION FOR LENNARD-JONES ATOMS. ** C ******************************************************************* SUBROUTINE FORCE ( EPSLON, SIGMA, RCUT, BOX, V, VC, W ) COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ, FX, FY, FZ C ******************************************************************* C ** FORCE CALCULATION FOR LENNARD-JONES ATOMS. ** C ** ** C ** IN THIS WE AIM TO SHOW HOW THE FORCES, POTENTIAL ENERGY AND ** C ** VIRIAL FUNCTION ARE CALCULATED IN A FAIRLY EFFICIENT WAY. ** C ** UNDOUBTEDLY FURTHER IMPROVEMENT WOULD BE POSSIBLE ON SPECIFIC ** C ** MACHINES. ** C ** THE POTENTIAL IS V(R) = 4*EPSLON*((SIGMA/R)**12-(SIGMA/R)**6) ** C ** WE INCLUDE SPHERICAL CUTOFF AND MINIMUM IMAGING IN CUBIC BOX. ** C ** THE BOX LENGTH IS BOX. THE CUTOFF IS RCUT. ** C ** THE ROUTINE ACTUALLY RETURNS TWO DIFFERENT POTENTIAL ENERGIES.** C ** V IS CALCULATED USING THE LENNARD-JONES POTENTIAL TO BE USED ** C ** FOR CALCULATING THE THERMODYNAMIC INTERNAL ENERGY. ** C ** LONG-RANGE CORRECTIONS SHOULD BE APPLIED TO THIS OUTSIDE THE ** C ** ROUTINE, IN THE FORM ** C ** SR3 = ( SIGMA / RCUT ) ** 3 ** C ** SR9 = SR3 ** 3 ** C ** DENS = REAL(N) * ( SIGMA / BOX ) ** 3 ** C ** VLRC = ( 8.0 /9.0 ) * PI * EPSLON * DENS * REAL ( N ) ** C ** : * ( SR9 - 3.0 * SR3 ) ** C ** WLRC = ( 16.0 / 9.0 ) * PI * EPSLON * DENS * REAL( N )** C ** : * ( 2.0 * SR9 - 3.0 * SR3 ) ** C ** V = V + VLRC ** C ** W = W + WLRC ** C ** VC IS CALCULATED USING THE SHIFTED LENNARD-JONES POTENTIAL, ** C ** WITH NO DISCONTINUITY AT THE CUTOFF, TO BE USED IN ASSESSING ** C ** ENERGY CONSERVATION. ** C ** NO REDUCED UNITS ARE USED: FOR THIS POTENTIAL WE COULD SET ** C ** EPSLON = 1 AND EITHER SIGMA = 1 OR BOX = 1 TO IMPROVE SPEED. ** C ** ** C ** PRINCIPAL VARIABLES: ** C ** ** C ** INTEGER N NUMBER OF MOLECULES ** C ** REAL RX(N),RY(N),RZ(N) MOLECULAR POSITIONS ** C ** REAL VX(N),VY(N),VZ(N) MOLECULAR VELOCITIES (NOT USED) ** C ** REAL FX(N),FY(N),FZ(N) MOLECULAR FORCES ** C ** REAL SIGMA PAIR POTENTIAL LENGTH PARAMETER ** C ** REAL EPSLON PAIR POTENTIAL ENERGY PARAMETER ** C ** REAL RCUT PAIR POTENTIAL CUTOFF ** C ** REAL BOX SIMULATION BOX LENGTH ** 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 SIGMA, EPSLON, RCUT, BOX, V, VC, W REAL RX(N), RY(N), RZ(N) REAL VX(N), VY(N), VZ(N) REAL FX(N), FY(N), FZ(N) INTEGER I, J, NCUT REAL BOXINV, RCUTSQ, SIGSQ, EPS4, EPS24 REAL RXI, RYI, RZI, FXI, FYI, FZI REAL RXIJ, RYIJ, RZIJ, RIJSQ, FXIJ, FYIJ, FZIJ REAL SR2, SR6, SR12, VIJ, WIJ, FIJ C ******************************************************************* C ** CALCULATE USEFUL QUANTITIES ** BOXINV = 1.0 / BOX RCUTSQ = RCUT ** 2 SIGSQ = SIGMA ** 2 EPS4 = EPSLON * 4.0 EPS24 = EPSLON * 24.0 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 = SIGSQ / RIJSQ SR6 = SR2 * SR2 * SR2 SR12 = SR6 ** 2 VIJ = SR12 - SR6 V = V + VIJ WIJ = VIJ + SR12 W = W + WIJ FIJ = WIJ / RIJSQ 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 = SIGSQ / RCUTSQ SR6 = SR2 * SR2 * SR2 SR12 = SR6 * SR6 VIJ = SR12 - SR6 VC = V - REAL ( NCUT ) * VIJ C ** MULTIPLY RESULTS BY ENERGY FACTORS ** DO 300 I = 1, N FX(I) = FX(I) * EPS24 FY(I) = FY(I) * EPS24 FZ(I) = FZ(I) * EPS24 300 CONTINUE V = V * EPS4 VC = VC * EPS4 W = W * EPS24 / 3.0 RETURN END