C ******************************************************************* C ** FICHE F.20 ** C ** CONSTRUCTION OF CELL LINKED-LISTS AND USE IN FORCE ROUTINE. ** C ******************************************************************* C ******************************************************************* C ** CONSTRUCTION OF CELL LINKED-LISTS AND USE IN FORCE ROUTINE. ** C ** ** C ** REFERENCES: ** C ** ** C ** QUENTREC AND BROT, J. COMPUT. PHYS. 13, 430, 1975. ** C ** HOCKNEY AND EASTWOOD, COMPUTER SIMULATION USING PARTICLES, ** C ** MCGRAW HILL, 1981. ** C ** ** C ** ROUTINES SUPPLIED: ** C ** ** C ** SUBROUTINE MAPS ** C ** SETS UP MAP OF CELL STRUCTURE FOR USE IN FORCE ** C ** SUBROUTINE LINKS ( RCUT ) ** C ** SETS UP HEAD OF CHAIN ARRAY AND LINKED LIST ** C ** SUBROUTINE FORCE ( SIGMA, RCUT, V, W ) ** C ** CALCULATES FORCES USING A LINKED LIST ** C ** ** C ** USAGE: ** C ** ** C ** SUBROUTINE MAPS IS CALLED ONCE AT THE START OF A SIMULATION ** C ** TO ESTABLISH CELL NEIGHBOUR IDENTITIES. AT EACH TIMESTEP, ** C ** SUBROUTINE LINKS IS CALLED TO SET UP THE LINKED LIST AND THIS ** C ** IS IMMEDIATELY USED BY SUBROUTINE FORCE. ** C ******************************************************************* SUBROUTINE MAPS COMMON / BLOCK2 / LIST, HEAD, MAP C ******************************************************************* C ** ROUTINE TO SET UP A LIST OF NEIGHBOURING CELLS ** C ** ** C ** PRINCIPAL VARIABLES: ** C ** ** C ** INTEGER M NUMBER OF CELLS IN EACH DIRECTION ** C ** INTEGER MAPSIZ SIZE OF CELL-CELL MAP ** C ** INTEGER MAP(MAPSIZ) LIST OF NEIGHBOURING CELLS ** C ** ** C ** USAGE: ** C ** ** C ** THIS SUBROUTINE SETS UP A LIST OF THE THIRTEEN NEIGHBOURING ** C ** CELLS OF EACH OF THE SMALL CELLS IN THE CENTRAL BOX. THE ** C ** EFFECTS OF THE PERIODIC BOUNDARY CONDITIONS ARE INCLUDED. ** C ** THE SUBROUTINE IS CALLED ONCE AT THE BEGINNING OF THE ** C ** SIMULATION AND THE MAP IS USED IN THE FORCE SUBROUTINE ** C ******************************************************************* INTEGER N, M, NCELL, MAPSIZ PARAMETER ( N = 1372, M = 5, NCELL = M * M * M ) PARAMETER ( MAPSIZ = 13 * NCELL ) INTEGER LIST(N), HEAD(NCELL), MAP(MAPSIZ) INTEGER IX, IY, IZ, IMAP, ICELL C ******************************************************************* C ** STATEMENT FUNCTION TO GIVE CELL INDEX ** ICELL ( IX, IY, IZ) = 1 + MOD ( IX - 1 + M, M ) : + MOD ( IY - 1 + M, M ) * M : + MOD ( IZ - 1 + M, M ) * M * M C ** FIND HALF THE NEAREST NEIGHBOURS OF EACH CELL ** DO 50 IZ = 1, M DO 40 IY = 1, M DO 30 IX = 1, M IMAP = ( ICELL ( IX, IY, IZ ) - 1 ) * 13 MAP( IMAP + 1 ) = ICELL( IX + 1, IY , IZ ) MAP( IMAP + 2 ) = ICELL( IX + 1, IY + 1, IZ ) MAP( IMAP + 3 ) = ICELL( IX , IY + 1, IZ ) MAP( IMAP + 4 ) = ICELL( IX - 1, IY + 1, IZ ) MAP( IMAP + 5 ) = ICELL( IX + 1, IY , IZ - 1 ) MAP( IMAP + 6 ) = ICELL( IX + 1, IY + 1, IZ - 1 ) MAP( IMAP + 7 ) = ICELL( IX , IY + 1, IZ - 1 ) MAP( IMAP + 8 ) = ICELL( IX - 1, IY + 1, IZ - 1 ) MAP( IMAP + 9 ) = ICELL( IX + 1, IY , IZ + 1 ) MAP( IMAP + 10 ) = ICELL( IX + 1, IY + 1, IZ + 1 ) MAP( IMAP + 11 ) = ICELL( IX , IY + 1, IZ + 1 ) MAP( IMAP + 12 ) = ICELL( IX - 1, IY + 1, IZ + 1 ) MAP( IMAP + 13 ) = ICELL( IX , IY , IZ + 1 ) 30 CONTINUE 40 CONTINUE 50 CONTINUE RETURN END SUBROUTINE LINKS ( RCUT ) COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ, FX, FY, FZ COMMON / BLOCK2 / LIST, HEAD, MAP C ******************************************************************* C ** ROUTINE TO SET UP LINKED LIST AND THE HEAD OF CHAIN ARRAYS ** C ** ** C ** PRINCIPAL VARIABLES: ** C ** ** C ** INTEGER N NUMBER OF ATOMS ** C ** INTEGER M NUMBER OF CELLS IN EACH DIRECTION ** C ** INTEGER NCELL TOTAL NUMBER OF CELLS (M**3) ** C ** INTEGER LIST(N) LINKED LIST OF ATOMS ** C ** INTEGER HEAD(NCELL) HEAD OF CHAIN FOR EACH CELL ** C ** REAL RX(N),RY(N),RZ(N) POSITIONS ** C ** REAL RCUT THE CUTOFF DISTANCE FOR THE FORCE ** C ** ** C ** USAGE: ** C ** ** C ** EACH ATOM IS SORTED INTO ONE OF THE M**3 SMALL CELLS. ** C ** THE FIRST ATOM IN EACH CELL IS PLACED IN THE HEAD ARRAY. ** C ** SUBSEQUENT ATOMS ARE PLACED IN THE LINKED LIST ARRAY. ** C ** ATOM COORDINATES ARE ASSUMED TO BE BETWEEN -0.5 AND +0.5. ** C ** THE ROUTINE IS CALLED EVERY TIMESTEP BEFORE THE FORCE ROUTINE.** C ******************************************************************* INTEGER N, M, NCELL, MAPSIZ PARAMETER ( N = 1372, M = 5, NCELL = M * M * M ) PARAMETER ( MAPSIZ = 13 * NCELL ) REAL RX(N), RY(N), RZ(N) REAL VX(N), VY(N), VZ(N) REAL FX(N), FY(N), FZ(N) INTEGER HEAD(NCELL), LIST(N), MAP(MAPSIZ) REAL CELLI, RCUT, CELL INTEGER ICELL, I C ******************************************************************* C ** ZERO HEAD OF CHAIN ARRAY ** DO 10 ICELL = 1, NCELL HEAD(ICELL) = 0 10 CONTINUE CELLI = REAL ( M ) CELL = 1.0 / CELLI IF ( CELL .LT. RCUT ) THEN STOP ' CELL SIZE TOO SMALL FOR CUTOFF ' ENDIF C ** SORT ALL ATOMS ** DO 20 I = 1, N ICELL = 1 + INT ( ( RX(I) + 0.5 ) * CELLI ) : + INT ( ( RY(I) + 0.5 ) * CELLI ) * M : + INT ( ( RZ(I) + 0.5 ) * CELLI ) * M * M LIST(I) = HEAD(ICELL) HEAD(ICELL) = I 20 CONTINUE RETURN END SUBROUTINE FORCE ( SIGMA, RCUT, V, W ) COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ, FX, FY, FZ COMMON / BLOCK2 / LIST, HEAD, MAP C ******************************************************************* C ** ROUTINE TO COMPUTE FORCES AND POTENTIAL USING A LINK LIST ** C ** ** C ** PRINCIPAL VARIABLES: ** C ** ** C ** INTEGER N NUMBER OF ATOMS ** C ** INTEGER M NUMBER OF CELLS IN EACH DIRECTION ** C ** INTEGER NCELL NUMBER OF SMALL CELLS (M**3) ** C ** INTEGER MAPSIZ SIZE OF CELL-CELL MAP ** C ** INTEGER LIST(N) THE LINKED LIST ** C ** INTEGER HEAD(NCELL) THE HEAD OF CHAIN ARRAY ** C ** INTEGER MAP(MAPSIZ) LIST OF NEIGHBOURING CELLS ** C ** REAL RX(N),RY(N),RZ(N) POSITIONS ** C ** REAL FX(N),FY(N),FZ(N) FORCES ** C ** REAL SIGMA THE LJ LENGTH PARAMETER ** C ** REAL RCUT THE CUT-OFF DISTANCE ** C ** REAL V THE POTENTIAL ENERGY ** C ** REAL W THE VIRIAL ** C ** ** C ** USAGE: ** C ** ** C ** FORCE IS CALLED IN AN MD PROGRAM TO CALCULATE THE FORCE ON ** C ** EACH ATOM. THE ROUTINE IS WRITTEN FOR A LIQUID OF LENNARD ** C ** JONES ATOMS. SUBROUTINE FORCE REQUIRES A LINKED LIST, SET UP ** C ** USING SUBROUTINE LINKS, AND THE MAP OF THE SMALL CELLS SET UP ** C ** USING SUBROUTINE MAPS. ** C ******************************************************************* INTEGER N, M, NCELL, MAPSIZ PARAMETER ( N = 1372, M = 5, NCELL = M * M * M ) PARAMETER ( MAPSIZ = 13 * NCELL) REAL RX(N), RY(N), RZ(N) REAL VX(N), VY(N), VZ(N) REAL FX(N), FY(N), FZ(N) INTEGER HEAD(NCELL), LIST(N), MAP(MAPSIZ) REAL RCUT, SIGMA, V, W REAL RXI, RYI, RZI, FXIJ, FYIJ, FZIJ, FIJ, RCUTSQ REAL SIGSQ, FXI, FYI, FZI, SR2, SR6, VIJ, WIJ REAL RIJSQ, RXIJ, RYIJ, RZIJ INTEGER ICELL, JCELL0, JCELL, I, J, NABOR C ******************************************************************* SIGSQ = SIGMA ** 2 RCUTSQ = RCUT ** 2 C ** ZERO FORCES AND POTENTIAL ** 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 C ** LOOP OVER ALL CELLS ** DO 5000 ICELL = 1, NCELL I = HEAD(ICELL) C ** LOOP OVER ALL MOLECULES IN THE CELL ** 1000 IF ( I .GT. 0 ) THEN RXI = RX(I) RYI = RY(I) RZI = RZ(I) FXI = FX(I) FYI = FY(I) FZI = FZ(I) C ** LOOP OVER ALL MOLECULES BELOW I IN THE CURRENT CELL ** J = LIST(I) 2000 IF ( J .GT. 0 ) THEN 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 = 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 ENDIF J = LIST(J) GO TO 2000 ENDIF C ** LOOP OVER NEIGHBOURING CELLS ** JCELL0 = 13 * (ICELL - 1) DO 4000 NABOR = 1, 13 JCELL = MAP ( JCELL0 + NABOR ) C ** LOOP OVER ALL MOLECULES IN NEIGHBOURING CELLS ** J = HEAD(JCELL) 3000 IF ( J .NE. 0 ) THEN 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 = 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 ENDIF J = LIST(J) GO TO 3000 ENDIF 4000 CONTINUE FX(I) = FXI FY(I) = FYI FZ(I) = FZI I = LIST(I) GO TO 1000 ENDIF 5000 CONTINUE C ** INCORPORATE ENERGY FACTORS ** V = 4.0 * V W = 48.0 * W / 3.0 DO 50 I = 1, N FX(I) = 48.0 * FX(I) FY(I) = 48.0 * FY(I) FZ(I) = 48.0 * FZ(I) 50 CONTINUE RETURN END