CCL Home Page
Up Directory CCL f.32
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.32                                                    **
C    ** ROUTINES TO IMPLEMENT CELL LINKED-LISTS IN SHEARED BOUNDARIES.**
C    *******************************************************************

C    *******************************************************************
C    ** ROUTINES TO IMPLEMENT CELL LINKED-LISTS IN SHEARED BOUNDARIES.**
C    **                                                               **
C    ** ROUTINES PROVIDED:                                            **
C    **                                                               **
C    ** SUBROUTINE MAPS                                               **
C    **    SETS UP MAP OF CELL NEIGHBOURS FOR BULK OF SIMULATION BOX  **
C    ** SUBROUTINE TOPMAP ( STRAIN )                                  **
C    **    SETS UP MAP OF NEIGHBOURS FOR TOP LAYER OF CELLS           **
C    ** SUBROUTINE LINKS ( RCUT )                                     **
C    **    CONSTRUCTS LINK-LIST GIVEN MAP OF CELL NEIGHBOURS          **
C    ** SUBROUTINE FORCE ( SIGMA, RCUT, STRAIN, V, W, WXY )           **
C    **    CALCULATES FORCES, POTENTIAL, VIRIAL ETC. USING LIST       **
C    **                                                               **
C    ** PRINCIPAL VARIABLES:                                          **
C    **                                                               **
C    ** INTEGER N                 NUMBER OF ATOMS                     **
C    ** REAL    RX(N),RY(N),RZ(N) ATOMIC POSITIONS                    **
C    ** REAL    VX(N),VY(N),VZ(N) ATOMIC VELOCITIES                   **
C    ** REAL    FX(N),FY(N),FZ(N) ATOMIC FORCES                       **
C    **                                                               **
C    ** USAGE:                                                        **
C    **                                                               **
C    ** SUBROUTINE MAPS IS CALLED ONCE AT THE START OF THE SIMULATION **
C    ** TO DEFINE CELL NEIGHBOURS FOR ALL BUT THE TOP LAYER OF CELLS. **
C    ** AT EACH TIME STEP, SUBROUTINE TOPMAP IS CALLED FOR THE TOP    **
C    ** LAYER, SUBROUTINE LINKS TO ESTABLISH THE ATOM NEIGHBOUR LIST, **
C    ** AND THEN THE FORCE SUBROUTINE.                                **
C    **                                                               **
C    ** UNITS:                                                        **
C    **                                                               **
C    ** THE PROGRAM ASSUMES A BOX OF UNIT LENGTH AND TAKES THE        **
C    ** LENNARD-JONES POTENTIAL WITH UNIT WELL-DEPTH.                 **
C    ** SUMMARY FOR BOX LENGTH L, ATOMIC MASS M, AND LENNARD-JONES    **
C    ** POTENTIAL PARAMETERS SIGMA AND EPSILON:                       **
C    **                                                               **
C    **                OUR PROGRAM            LENNARD-JONES SYSTEM    **
C    ** LENGTH         L                      SIGMA                   **
C    ** MASS           M                      M                       **
C    ** ENERGY         EPSILON                EPSILON                 **
C    ** TIME           SQRT(M*L**2/EPSILON)   SQRT(M*SIGMA**2/EPSILON)**
C    ** VELOCITY       SQRT(EPSILON/M)        SQRT(EPSILON/M)         **
C    ** PRESSURE       EPSILON/L**3           EPSILON/SIGMA**3        **
C    *******************************************************************



        SUBROUTINE MAPS

        COMMON / BLOCK2 / LIST, HEAD, MAP

C    *******************************************************************
C    ** CONSTRUCTS MAP OF CELL NEIGHBOURS.                            **
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    ** HOWEVER THE TOP LAYER (IY = M) IS TACKLED SEPARATELY BECAUSE  **
C    ** OF THE SHEARED BOUNDARY CONDITIONS, IN SUBROUTINE TOPMAP.     **
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    ** 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, M3
        PARAMETER ( N = 1372 )
        PARAMETER ( M = 5, NCELL = M * M * M, MAPSIZ = 16 * NCELL )
        PARAMETER ( M3 = M * 3 )

        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 + M3, M )
     :                           + MOD ( IY - 1 + M3, M ) * M
     :                           + MOD ( IZ - 1 + M3, M ) * M * M

C    ** FIND HALF THE NEAREST NEIGHBOURS OF EACH CELL **

        DO 50 IZ = 1, M

           DO 40 IY = 1, M - 1

              DO 30 IX = 1, M

                 IMAP = ( ICELL ( IX, IY, IZ ) - 1 ) * 16

                 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 )
                 MAP( IMAP + 14 ) = 0
                 MAP( IMAP + 15 ) = 0
                 MAP( IMAP + 16 ) = 0

30            CONTINUE

40         CONTINUE

50      CONTINUE

        RETURN
        END



        SUBROUTINE TOPMAP ( STRAIN )

        COMMON / BLOCK2 / LIST, HEAD, MAP

C    *******************************************************************
C    ** CALCULATES CELL NEIGHBOUR MAP FOR TOP LAYER.                  **
C    **                                                               **
C    ** THIS SUBROUTINE SUPPLEMENTS THE LIST OF NEIGHBOURING CELLS    **
C    ** FOR THE TOP LAYER (IY = M) WITH SHEARED BOUNDARY CONDITIONS   **
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    ** REAL    STRAIN             THE X-DISPLACEMENT OF NEXT BOX UP  **
C    **                                                               **
C    ** USAGE:                                                        **
C    **                                                               **
C    ** THE SUBROUTINE IS CALLED AT EVERY TIMESTEP IN THE SIMULATION  **
C    ** JUST BEFORE THE FORCE SUBROUTINE                              **
C    *******************************************************************

        INTEGER     N, M, NCELL, MAPSIZ, M3
        PARAMETER ( N = 1372 )
        PARAMETER ( M = 5, NCELL = M * M * M, MAPSIZ = 16 * NCELL )
        PARAMETER ( M3 = M * 3 )

        REAL        STRAIN

        INTEGER     LIST(N), HEAD(NCELL), MAP(MAPSIZ)
        INTEGER     IX, IY, IZ, IMAP, ICELL, IIX

C    *******************************************************************

C    ** STATEMENT FUNCTION TO GIVE CELL INDEX **

        ICELL ( IX, IY, IZ ) = 1 + MOD ( IX - 1 + M3, M )
     :                           + MOD ( IY - 1 + M3, M ) * M
     :                           + MOD ( IZ - 1 + M3, M ) * M * M

C    ** CALCULATE X OFFSET IN CELL LENGTHS WHERE STRAIN **
C    ** IS BETWEEN -1/2 AND +1/2 AND BOX LENGTH = 1.0   **
C    ** ADDING 1.0 SIMPLY GUARANTEES A POSITIVE RESULT  **

        STRAIN = STRAIN - ANINT ( STRAIN )
        IIX   = INT ( ( STRAIN + 1.0 ) * REAL ( M ) )

C    ** FIND HALF THE NEAREST NEIGHBOURS OF EACH CELL **

        IY = M

        DO 50 IZ = 1, M

           DO 30 IX = 1, M

              IMAP = ( ICELL ( IX, IY, IZ ) - 1 ) * 16

              MAP( IMAP + 1  ) = ICELL ( IX + 1      , IY    , IZ     )
              MAP( IMAP + 2  ) = ICELL ( IX + 1 - IIX, IY + 1, IZ     )
              MAP( IMAP + 3  ) = ICELL ( IX     - IIX, IY + 1, IZ     )
              MAP( IMAP + 4  ) = ICELL ( IX - 1 - IIX, IY + 1, IZ     )
              MAP( IMAP + 5  ) = ICELL ( IX + 1      , IY    , IZ - 1 )
              MAP( IMAP + 6  ) = ICELL ( IX + 1 - IIX, IY + 1, IZ - 1 )
              MAP( IMAP + 7  ) = ICELL ( IX     - IIX, IY + 1, IZ - 1 )
              MAP( IMAP + 8  ) = ICELL ( IX - 1 - IIX, IY + 1, IZ - 1 )
              MAP( IMAP + 9  ) = ICELL ( IX + 1      , IY    , IZ + 1 )
              MAP( IMAP + 10 ) = ICELL ( IX + 1 - IIX, IY + 1, IZ + 1 )
              MAP( IMAP + 11 ) = ICELL ( IX     - IIX, IY + 1, IZ + 1 )
              MAP( IMAP + 12 ) = ICELL ( IX - 1 - IIX, IY + 1, IZ + 1 )
              MAP( IMAP + 13 ) = ICELL ( IX          , IY    , IZ + 1 )
              MAP( IMAP + 14 ) = ICELL ( IX - 2 - IIX, IY + 1, IZ     )
              MAP( IMAP + 15 ) = ICELL ( IX - 2 - IIX, IY + 1, IZ - 1 )
              MAP( IMAP + 16 ) = ICELL ( IX - 2 - IIX, IY + 1, IZ + 1 )

30         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    ** SUBROUTINE TO SET UP LINKED LIST AND THE HEAD OF CHAIN ARRAYS **
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    **                                                               **
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    ** THE ROUTINE IS CALLED EVERY TIMESTEP BEFORE THE FORCE ROUTINE.**
C    *******************************************************************

        INTEGER     N, M, NCELL, MAPSIZ
        PARAMETER ( N = 1372 )
        PARAMETER ( M = 5, NCELL = M * M * M, MAPSIZ = 16 * 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

           WRITE(*,'('' CELL SIZE TOO SMALL FOR CUTOFF '')')
           STOP

        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, STRAIN, V, W, WXY )

        COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ, FX, FY, FZ
        COMMON / BLOCK2 / LIST, HEAD, MAP

C    *******************************************************************
C    ** COMPUTES FORCES, ETC. USING A LINK LIST IN SHEARED BOUNDARIES.**
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    STRAIN             X OFFSET OF SUCCESSIVE BOXES       **
C    ** REAL    V                  THE POTENTIAL ENERGY               **
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 SUBROUTINES MAPS AND TOPMAP.                            **
C    *******************************************************************

        INTEGER     N, M, NCELL, MAPSIZ
        PARAMETER ( N = 1372 )
        PARAMETER ( M = 5, NCELL = M * M * M, MAPSIZ = 16 * 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, STRAIN, V, W, WXY

        REAL        RXI, RYI, RZI, FXIJ, FYIJ, FZIJ, RCUTSQ
        REAL        VIJ, WIJ, FIJ
        REAL        SIGSQ, FXI, FYI, FZI, SR2, SR6, SR12
        REAL        RIJSQ, RXIJ, RYIJ, RZIJ, CORY
        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
        WXY  = 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)
                 CORY  = ANINT ( RYIJ )
                 RXIJ  = RXIJ - CORY * STRAIN
                 RXIJ  = RXIJ - ANINT ( RXIJ )
                 RYIJ  = RYIJ - CORY
                 RZIJ  = RZIJ - ANINT ( RZIJ )
                 RIJSQ = RXIJ * RXIJ + RYIJ * RYIJ + RZIJ * RZIJ

                 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
                    WXY   = WXY + RXIJ * FYIJ

                 ENDIF

                 J = LIST(J)
                 GOTO 2000

              ENDIF

C          ** LOOP OVER NEIGHBOURING CELLS **

              JCELL0 = 16 * ( ICELL - 1 )

              DO 4000 NABOR = 1, 16

                 JCELL = MAP( JCELL0 + NABOR )

                 IF ( JCELL .GT. 0 ) THEN

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)
                       CORY  = ANINT ( RYIJ )
                       RXIJ  = RXIJ - CORY * STRAIN
                       RXIJ  = RXIJ - ANINT ( RXIJ )
                       RYIJ  = RYIJ - CORY
                       RZIJ  = RZIJ - ANINT ( RZIJ )
                       RIJSQ = RXIJ * RXIJ + RYIJ * RYIJ + RZIJ * RZIJ

                       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
                          WXY   = WXY + RXIJ * FYIJ

                       ENDIF

                       J = LIST(J)
                       GOTO 3000

                    ENDIF

                 ENDIF

4000          CONTINUE

              FX(I) = FXI
              FY(I) = FYI
              FZ(I) = FZI

              I = LIST(I)
              GOTO 1000

           ENDIF

5000    CONTINUE

C    ** INCORPORATE ENERGY FACTORS **

        DO 6000 I = 1, N

           FX(I) = FX(I) * 24.0
           FY(I) = FY(I) * 24.0
           FZ(I) = FZ(I) * 24.0

6000    CONTINUE

        V   = V   * 4.0
        W   = W   * 24.0 / 3.0
        WXY = WXY * 24.0

        RETURN
        END


Modified: Fri Oct 15 16:00:00 1993 GMT
Page accessed 5570 times since Sat Aug 26 23:03:08 2000 GMT