CCL Home Page
Up Directory CCL f.20
C    *******************************************************************
C    ** THIS FORTRAN CODE IS INTENDED TO ILLUSTRATE POINTS MADE IN    **
C    ** THE TEXT. TO OUR KNOWLEDGE IT WORKS CORRECTLY. HOWEVER IT IS  **
C    ** THE RESPONSIBILITY OF THE USER TO TEST IT, IF IT IS USED IN A **
C    ** RESEARCH APPLICATION.                                         **
C    *******************************************************************

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


Modified: Fri Oct 15 16:00:00 1993 GMT
Page accessed 5628 times since Sat Aug 26 23:02:53 2000 GMT