C ******************************************************************* C ** FICHE F.19 ** C ** FORCE ROUTINE USING A VERLET NEIGHBOUR LIST. ** C ******************************************************************* C ******************************************************************* C ** FORCE ROUTINE USING A VERLET NEIGHBOUR LIST. ** C ** ** C ** PRINCIPAL VARIABLES: ** C ** ** C ** INTEGER N NUMBER OF ATOMS ** C ** REAL RCUT CUTOFF DISTANCE FOR THE FORCE ** C ** REAL RLIST OUTER RADIUS OF THE VERLET LIST ** C ** REAL SIGMA LENNARD JONES LENGTH PARAMETER ** C ** REAL V POTENTIAL ENERGY ** C ** REAL W VIRIAL ** C ** REAL RX(N),RY(N),RZ(N) ATOM POSITIONS ** C ** REAL FX(N),FY(N),FZ(N) FORCE ON AN ATOM ** C ** LOGICAL UPDATE IF TRUE THE LIST IS UPDATED ** C ** ** C ** ROUTINES SUPPLIED: ** C ** ** C ** SUBROUTINE FORCE ( RCUT, RLIST, SIGMA, UPDATE, V, W ) ** C ** CALCULATE FORCES USING THE VERLET LIST, UPDATES THE LIST ** C ** SUBROUTINE CHECK ( RCUT, RLIST, UPDATE ) ** C ** SETS UPDATE TO TRUE WHEN THE LIST NEEDS TO BE UPDATED ** C ** SUBROUTINE SAVE ** C ** SAVES CURRENT CONFIGURATION FOR FUTURE CHECKING. ** C ** ** C ** USAGE: ** C ** ** C ** AT THE START OF A RUN, FORCE IS CALLED WITH UPDATE SET TRUE. ** C ** THIS SETS UP THE INITIAL VERLET LIST AND SAVES THE POSITIONS. ** C ** THEREAFTER CHECK IS CALLED JUST BEFORE EACH CALL OF FORCE TO ** C ** DECIDE WHETHER OR NOT A LIST UPDATE IS NECESSARY. ** C ** THESE ROUTINES CAN BE USED IN A CONVENTIONAL MD PROGRAM AND ** C ** CAN BE EASILY ADAPTED FOR USE IN AN MC PROGRAM. FORCES IS ** C ** SPECIFIC TO A FLUID OF LENNARD JONES ATOMS. ** C ******************************************************************* SUBROUTINE FORCE ( RCUT, RLIST, SIGMA, UPDATE, V, W ) COMMON / BLOCK1 / RX, RY, RZ, FX, FY, FZ COMMON / BLOCK2 / POINT, LIST C ******************************************************************* C ** CALCULATES THE FORCE ON AN ATOM IN A LENNARD-JONES LIQUID ** C ** ** C ** PRINCIPAL VARIABLES: ** C ** ** C ** REAL RCUT CUTOFF DISTANCE FOR THE FORCE ** C ** REAL RLIST OUTER RADIUS OF THE VERLET LIST ** C ** REAL SIGMA LENNARD JONES LENGTH PARAMETER ** C ** REAL V POTENTIAL ENERGY ** C ** REAL W VIRIAL ** C ** REAL RX(N),RY(N),RZ(N) ATOM POSITIONS ** C ** REAL FX(N),FY(N),FZ(N) FORCE ON AN ATOM ** C ** INTEGER POINT(N) INDEX TO THE NEIGHBOUR LIST ** C ** INTEGER LIST(MAXNAB) VERLET NEIGHBOUR LIST ** C ** LOGICAL UPDATE IF TRUE THE LIST IS UPDATED ** C ** ** C ** USAGE: ** C ** ** C ** FORCE IS CALLED IN TWO MODES. IF UPDATE IS TRUE, THE VERLET ** C ** LIST IS UPDATED AND THE FORCES CALCULATED AT THE SAME TIME. ** C ** THE CURRENT CONFIGURATION IS SAVED FOR FUTURE CHECKING. ** C ** IF UPDATE IS FALSE, THE VERLET LIST IS USED TO FIND THE ** C ** NEIGHBOURS OF A GIVEN ATOM AND CALCULATE THE FORCES. ** C ** THE ATOMS ARE IN A BOX OF UNIT LENGTH, CENTRED AT THE ORIGIN. ** C ******************************************************************* INTEGER N, MAXNAB PARAMETER ( N = 108, MAXNAB = N * 25 ) REAL RX(N), RY(N), RZ(N), FX(N), FY(N), FZ(N) REAL RLIST, RCUT, SIGMA, V, W INTEGER POINT(N), LIST(MAXNAB) LOGICAL UPDATE REAL RXI, RYI, RZI, FXI, FYI, FZI REAL RIJSQ, SR2, SR6, VIJ, WIJ, FIJ REAL SIGSQ, RCUTSQ, RLSTSQ REAL RXIJ, RYIJ, RZIJ, FXIJ, FYIJ, FZIJ INTEGER I, J, NLIST INTEGER JBEG, JEND, JNAB C ******************************************************************* RLSTSQ = RLIST * RLIST SIGSQ = SIGMA * SIGMA RCUTSQ = RCUT * RCUT C ** ZERO FORCES ** DO 10 I = 1, N FX(I) = 0.0 FY(I) = 0.0 FZ(I) = 0.0 10 CONTINUE V = 0.0 W = 0.0 IF ( UPDATE ) THEN C ** SAVE CURRENT CONFIGURATION, CONSTRUCT ** C ** NEIGHBOUR LIST AND CALCULATE FORCES ** CALL SAVE NLIST = 0 DO 100 I = 1, N - 1 POINT(I) = NLIST + 1 RXI = RX(I) RYI = RY(I) RZI = RZ(I) FXI = FX(I) FYI = FY(I) FZI = FZ(I) DO 99 J = I + 1, N RXIJ = RXI - RX(J) RYIJ = RYI - RY(J) RZIJ = RZI - RZ(J) RXIJ = RXIJ - ANINT ( RXIJ ) RYIJ = RYIJ - ANINT ( RYIJ ) RZIJ = RZIJ - ANINT ( RZIJ ) RIJSQ = RXIJ * RXIJ + RYIJ * RYIJ + RZIJ * RZIJ IF ( RIJSQ .LT. RLSTSQ ) THEN NLIST = NLIST + 1 LIST(NLIST) = J C ** REMOVE THIS CHECK IF MAXNAB IS APPROPRIATE ** IF ( NLIST .EQ. MAXNAB ) STOP 'LIST TOO SMALL' IF ( RIJSQ .LT. RCUTSQ ) THEN SR2 = SIGSQ / RIJSQ SR6 = SR2 * SR2 * SR2 VIJ = SR6 * ( SR6 - 1.0 ) WIJ = SR6 * ( SR6 - 0.5 ) V = V + VIJ W = W + WIJ FIJ = WIJ / RIJSQ FXIJ = RXIJ * FIJ FYIJ = RYIJ * FIJ FZIJ = RZIJ * FIJ FXI = FXI + FXIJ FYI = FYI + FYIJ FZI = FZI + FZIJ FX(J) = FX(J) - FXIJ FY(J) = FY(J) - FYIJ FZ(J) = FZ(J) - FZIJ ENDIF ENDIF 99 CONTINUE FX(I) = FXI FY(I) = FYI FZ(I) = FZI 100 CONTINUE POINT(N) = NLIST + 1 ELSE C ** USE THE LIST TO FIND THE NEIGHBOURS ** DO 200 I = 1, N - 1 JBEG = LIST(I) JEND = LIST(I+1) - 1 C ** CHECK THAT ATOM I HAS NEIGHBOURS ** IF( JBEG .LE. JEND ) THEN RXI = RX(I) RYI = RY(I) RZI = RZ(I) FXI = FX(I) FYI = FY(I) FZI = FZ(I) DO 199 JNAB = JBEG, JEND J = LIST(JNAB) RXIJ = RXI - RX(J) RYIJ = RYI - RY(J) RZIJ = RZI - RZ(J) RXIJ = RXIJ - ANINT( RXIJ ) RYIJ = RYIJ - ANINT( RYIJ ) RZIJ = RZIJ - ANINT( RZIJ ) RIJSQ = RXIJ * RXIJ + RYIJ * RYIJ + RZIJ * RZIJ IF ( RIJSQ .LT. RCUTSQ ) THEN SR2 = SIGSQ / RIJSQ SR6 = SR2 * SR2 * SR2 VIJ = SR6 * ( SR6 - 1.0 ) WIJ = SR6 * ( SR6 - 0.5 ) V = V + VIJ W = W + WIJ FIJ = WIJ / RIJSQ FXIJ = RXIJ * FIJ FYIJ = RYIJ * FIJ FZIJ = RZIJ * FIJ FXI = FXI + FXIJ FYI = FYI + FYIJ FZI = FZI + FZIJ FX(J) = FX(J) - FXIJ FY(J) = FY(J) - FYIJ FZ(J) = FZ(J) - FZIJ ENDIF 199 CONTINUE FX(I) = FXI FY(I) = FYI FZ(I) = FZI ENDIF 200 CONTINUE ENDIF V = 4.0 * V W = 48.0 * W / 3.0 DO 300 I = 1, N FX(I) = 48.0 * FX(I) FY(I) = 48.0 * FY(I) FZ(I) = 48.0 * FZ(I) 300 CONTINUE RETURN END SUBROUTINE CHECK ( RCUT, RLIST, UPDATE ) COMMON / BLOCK1 / RX, RY, RZ, FX, FY, FZ COMMON / BLOCK3 / RX0, RY0, RZ0 C ******************************************************************* C ** DECIDES WHETHER THE LIST NEEDS TO BE RECONSTRUCTED. ** C ** ** C ** PRINCIPAL VARIABLES: ** C ** ** C ** REAL RX(N),RY(N),RZ(N) ATOM POSITIONS ** C ** REAL RX0(N),RY0(N),RZ0(N) COORDINATES AT LAST UPDATE ** C ** REAL RLIST RADIUS OF VERLET LIST ** C ** REAL RCUT CUTOFF DISTANCE FOR FORCES ** C ** REAL DISPMX LARGEST DISPLACEMENT ** C ** INTEGER N NUMBER OF ATOMS ** C ** LOGICAL UPDATE IF TRUE THE LIST IS UPDATED ** C ** ** C ** USAGE: ** C ** ** C ** CHECK IS CALLED TO SET UPDATE BEFORE EVERY CALL TO FORCE. ** C ******************************************************************* INTEGER N PARAMETER ( N = 108 ) REAL RX(N), RY(N), RZ(N), FX(N), FY(N), FZ(N) REAL RX0(N), RY0(N), RZ0(N) REAL RCUT, RLIST LOGICAL UPDATE REAL DISPMX INTEGER I C ******************************************************************* C ** CALCULATE MAXIMUM DISPLACEMENT SINCE LAST UPDATE ** DISPMX = 0.0 DO 30 I = 1, N DISPMX = MAX ( ABS ( RX(I) - RX0(I) ), DISPMX ) DISPMX = MAX ( ABS ( RY(I) - RY0(I) ), DISPMX ) DISPMX = MAX ( ABS ( RZ(I) - RZ0(I) ), DISPMX ) 30 CONTINUE C ** A CONSERVATIVE TEST OF THE LIST SKIN CROSSING ** DISPMX = 2.0 * SQRT ( 3.0 * DISPMX ** 2 ) UPDATE = ( DISPMX .GT. ( RLIST - RCUT ) ) RETURN END SUBROUTINE SAVE COMMON / BLOCK1 / RX, RY, RZ, FX, FY, FZ COMMON / BLOCK3 / RX0, RY0, RZ0 C ******************************************************************* C ** SAVES CURRENT CONFIGURATION FOR FUTURE CHECKING. ** C ** ** C ** PRINCIPAL VARIABLES: ** C ** ** C ** REAL RX(N),RY(N),RZ(N) ATOM POSITIONS ** C ** REAL RX0(N),RY0(N),RZ0(N) STORAGE LOCATIONS FOR SAVE ** C ** INTEGER N NUMBER OF ATOMS ** C ** ** C ** USAGE: ** C ** ** C ** SAVE IS CALLED WHENEVER THE NEW VERLET LIST IS CONSTRUCTED. ** C ******************************************************************* INTEGER N PARAMETER ( N = 108 ) REAL RX(N), RY(N), RZ(N), FX(N), FY(N), FZ(N) REAL RX0(N), RY0(N), RZ0(N) INTEGER I C ******************************************************************* DO 100 I = 1, N RX0(I) = RX(I) RY0(I) = RY(I) RZ0(I) = RZ(I) 100 CONTINUE RETURN END