C    *******************************************************************
C    **                                                               **
C    ** MOLECULAR DYNAMICS FOR A LENNARD JONES FLUID                  **
C    ** 3D VERSION PASSING EDGES IN X,Y,Z DIRECTIONS                  **
C    **                                                               **
C    ** AS FAR AS WE ARE AWARE THIS PROGRAM WORKS CORRECTLY, BUT WE   **
C    ** CAN ACCEPT NO RESPONSIBILITY FOR THE CONSEQUNCES OF ANY       **
C    ** ERRORS, AND WOULD BE GRATEFUL TO HEAR FROM YOU IF YOU FIND    **
C    ** ANY. YOU SHOULD ALWAYS CHECK A PROGRAM OR SUBROUTINE FOR YOUR **
C    ** PARTICULAR APPLICATION.                                       **
C    **                                                               **
C    ** THIS PROGRAM WILL SIMULATE A THREE DIMENSIONAL LENNARD-JONES  **
C    ** FLUID. THE PROGRAM USES A THREE DIMENSIONAL SPATIAL           **
C    ** DECOMPOSITION OF THE SIMULATION SPACE INTO SUB-CUBES.         **
C    **                                                               **
C    ** SUBROUTINES CALLED:                                           **
C    **                                                               **
C    ** FORCE     CALCULATES THE  FORCES AND POTENTIAL USING THE      **
C    **           LINKED LIST.                                        **
C    ** MAPS      CREATES NEIGHBOUR LIST OF CELLS USED BY FORCE.      **
C    ** MOVEA     PERFORMS FIRST PART OF VELOCITY VERLET ALGORITHM.   **
C    ** MOVEB     PERFORMS SECOND PART OF VELOCITY VERLET ALGORITHM.  **
C    ** MOVOUT    CREATES LINKED LIST.
C    ** SCALET    SCALES THE VELOCITIES TO ENABLE EQUILIBRATION TO A  **
C    **           PARTICULAR TEMPERATURE.                             **
C    **                                                               **
C    ** PRINCIPAL VARIABLES:                                          **
C    **                                                               **
C    ** INTEGER VARIABLES:                                            **
C    **                                                               **
C    ** HEAD(NCELL)                  ARRAY HOLDING INITIAL POINTER FOR**
C    **                              EACH LINK-CELL.                  **
C    ** ISAVE                        NUMBER OF TIMESTEPS BETWEEN SAVES**
C    **                              OF CONFIGURATIONAL INFORMATION.  **
C    ** LIST(NMAX)                   ARRAY HOLDING LINKED LIST.       **
C    ** MAP(MAPSIZ)                  ARRAY CONTAINING NEIGHBOUR LIST  **
C    **                              OF CELLS CONSIDERED IN FORCE     **
C    **                              CALCULATION.                     **
C    ** MAPSIZ                       SIZE OF ARRAY MAP.               **
C    ** MMX,MMY,MMZ                  MAXIMUM NUMBER OF LINK-CELLS IN  **
C    **                              EACH COORDINATE DIRECTION.       **
C    ** MX,MY,MZ                     NUMBER OF LINK CELLS IN EACH     **
C    **                              COORDINATE DIRECTION.            **
C    ** NATM                         NUMBER OF REAL ATOMS CURRENTLY   **
C    **                              CONTROLLED BY THIS PROCESSOR.    **
C    ** NCELL                        MAXIMUM TOTAL NUMBER OF LINK-CELLS*
C    ** NMAX                         SIZE OF COORDINATE ARRAYS AND    **
C    **                              LIST ARRAY. N.B. SHOULD BE LARGE **
C    **                              ENOUGH TO HOLD ALL REAL AND EDGE **
C    **                              COORDINATES.                     **
C    ** NEQUIL                       NUMBER OF EQUILIBRATION TIMESTEPS**
C    **                              TO PERFORM.                      **
C    ** NSTEP                        NUMBER OF TIMESTEPS TO PERFORM.  **
C    **                                                               **
C    ** REAL VARIABLES:                                               **
C    **                                                               **
C    ** CELLIX,CELLIY,CELLIZ         REAL VALUES OF MX, MY, AND MZ    **
C    ** DT                           TIMESTEP.                        **
C    ** EQTEMP                       TEMPERATURE TO WHICH THE SYSTEM  **
C    **                              IS BEING EQUILIBRATED.           **
C    ** RX(NMAX),RY(NMAX),RZ(NMAX)   ARRAYS WHICH HOLD COORDINATES.   **
C    ** RCUT                         CUTOFF                           **
C    ** SFX,SFY,SFZ                  SCALE FACTORS TO ENABLE LINKED   **
C    **                              LIST CALCULATION FROM PROCESSOR  **
C    **                              COORDINATE SYSTEM.               **
C    ** SIGMA                        SIGMA FOR LENNARD-JONES POTENTIAL**
C    ** VX(NMAX),VY(NMAX),VZ(NMAX)    ARRAYS WHICH HOLD VELOCITIES.   **
C    ** V,K,W                        TOTAL POTENTIAL AND KINETIC ENERGY*
C    **                              AND VIRIAL FOR ATOMS CONTROLLED  **
C    **                              BY THIS PROCESSOR.               **
C    **                                                               **
C    *******************************************************************

      PROGRAM MD1

      IMPLICIT CHARACTER ( A-Z )

      DOUBLE PRECISION SIGMA, RCUT
      COMMON /POTNTL/ SIGMA, RCUT

      INTEGER NW1
      PARAMETER (NW1 = 11)

      INCLUDE 'md1.inc'

      DOUBLE PRECISION RX (NMAX), RY (NMAX), RZ (NMAX)
      DOUBLE PRECISION VX (NMAX), VY (NMAX), VZ (NMAX)
      DOUBLE PRECISION FX (NMAX), FY (NMAX), FZ (NMAX)
      DOUBLE PRECISION DT
      DOUBLE PRECISION DENS, DENLJ, EQTEMP
      DOUBLE PRECISION MASS, FREE, PI, TEMP
      DOUBLE PRECISION E, V, K, W
      DOUBLE PRECISION EN, VN, KN, PRES
      DOUBLE PRECISION ACE, ACV, ACK, ACP
      DOUBLE PRECISION ACESQ, ACVSQ, ACKSQ, ACPSQ
      DOUBLE PRECISION AVE, AVV, AVK, AVP, AVT
      DOUBLE PRECISION FLE, FLV, FLK, FLP, FLT

      DOUBLE PRECISION TSTART, TEND, TIME

      INTEGER HEAD (MMX, MMY, MMZ), LIST (NMAX)
      INTEGER MX,MY,MZ
      INTEGER NC,M
      INTEGER STEP, NSTEP, ISAVE, IPRINT, NORM, NATM
      INTEGER NEQUIL, ISCALE, NPROSX, NPROSY, NPROSZ
      INTEGER I

      DO I = 1, NMAX
        FX (I) = 0.0D0 
        FY (I) = 0.0D0
        FZ (I) = 0.0D0
      ENDDO

      NATM = 0

      PI = ACOS(-1.D0)

C    ** READ INPUT DATA **

      WRITE(*,'(1H1,'' **** PROGRAM MDLJ ****                   '')')
      WRITE(*,'('' MOLECULAR DYNAMICS SIMULATION                '')')
      WRITE(*,'('' WITH LINKED LIST                             '')')

C     NUMBER OF TIME STEPS

      NSTEP = 600

C     NUMBER OF EQUILIBRATION TIME STEPS

      NEQUIL = 500

C     NUMBER OF STEPS BETWEEN DATA SAVES

      ISAVE = 2*NSTEP

C     NUMBER OF STEPS BETWEEN OUTPUT

      IPRINT = 100

C     THE FOLLOWING IN LENNARD-JONES UNITS...

C     THE TEMPERATURE

      EQTEMP = 0.8D0

C     THE DENSITY

      DENS = 0.75D0

C     THE POTENTIAL CUTOFF DISTANCE

      RCUT = 2.0D0**(1.0D0/6.0D0)

C     THE TIMESTEP

      DT = 0.01D0

C     REDUNDANT VARIABLES...

      ISCALE = 0
      NPROSX = 1
      NPROSY = 1
      NPROSZ = 1

C     WRITE INPUT DATA

      OPEN (NW1,FILE='result')
      REWIND(NW1)

C     READ NUMBER OF FCC UNIT-CELLS TO USED IN ALL THREE
C     COORDINATE DIRECTIONS

      WRITE(*,'('' ENTER NC '')')
      READ(*,*) NC

      IF (4*NC*NC*NC .GT. NMAX) STOP ' Too many atoms...'

      CALL FCC (NATM, RX, RY, RZ, VX, VY, VZ, NC)

      WRITE(*,'('' NUMBER OF ATOMS BEING USED '',I8)')NATM

      CALL HEADER (NW1, NPROSX, NPROSY, NPROSZ, NC, NATM)

      CALL PARAM(NW1,NSTEP,NEQUIL,IPRINT,ISCALE,EQTEMP,DENS,RCUT,DT)

C     CONVERT INPUT DATA TO PROGRAM UNITS

      SIGMA  = (DENS/DBLE(NATM))**(1.0D0/3.0D0)
      RCUT   = RCUT*SIGMA
      M      = IDINT(1.D0/RCUT)
      MASS   = 1.0D0
      DENLJ  = DENS
      DENS   = DENS/(SIGMA**3)
      DT     = DT*SIGMA
      FREE   = DBLE(3*(NATM-1))

C     CHECK TO SEE WHETHER SYSTEM CAN BE SIMULATED

      IF (M.LT.1)  STOP 'SYTEM TOO SMALL FOR ARRAY'

C     ZERO ACCUMULATORS

      ACE    = 0.0D0
      ACV    = 0.0D0
      ACK    = 0.0D0
      ACP    = 0.0D0

      ACESQ  = 0.0D0
      ACVSQ  = 0.0D0
      ACKSQ  = 0.0D0
      ACPSQ  = 0.0D0

      FLE    = 0.0D0
      FLV    = 0.0D0
      FLK    = 0.0D0
      FLP    = 0.0D0
      FLT    = 0.0D0

C     WRITE OUT SOME USEFUL INFORMATION

      WRITE(*,'('' SIGMA/BOX              =  '',F10.4)')  SIGMA
      WRITE(*,'('' RCUT/BOX               =  '',F10.4)')  RCUT
      WRITE(*,'('' DT                     =  '',F10.4)')  DT

C     CALCULATE DIMENSIONS OF PROCESSORS LINK-CELL STRUCTURE

      MX = MIN (MMX, IDINT (1.0D0 / RCUT))
      MY = MIN (MMY, IDINT (1.0D0 / RCUT))
      MZ = MIN (MMZ, IDINT (1.0D0 / RCUT))

C     CALL SUBROUTINE TO CALCULATE NEIGHBOUR LIST OF NEIGHBOURING
C     CELLS THAT WE WISH TO CONSIDER IN FORCE LOOP

      CALL MAPS (MX, MY, MZ)

C     ******* DEBUG *******
      WRITE (6, *) 'MX, MY, MZ =', MX, MY, MZ

C     RECEIVE CONFIGURATIONAL INFORMATION FROM MASTER
C     LOOPS OVER ALL STEPS

      WRITE(*,'(/'' ** MOLECULAR DYNAMICS BEGINS ** ''/ )')
      WRITE(*,10001)
      WRITE(NW1,'(/'' ** MOLECULAR DYNAMICS BEGINS ** ''/ )')
      WRITE(NW1,10001)

C     START OF MD CYCLE

      DO STEP=1,MIN(NSTEP,NEQUIL)
        V = 0
        K = 0
        W = 0
        PRES = 0

C       PERFORM FIRST PART OF VELOCITY VERLET INTEGRATION

        CALL MOVEA (NATM, RX, RY, RZ, VX, VY, VZ, FX, FY, FZ, DT)

C       CREATE LINKED LISTS

        CALL MOVOUT (NATM, LIST, MX, MY, MZ, HEAD, RX, RY, RZ)

C       CALCULATE FORCES AND POTENTIAL USING LINKED SCALAR ROUTINE
C       MORE THAN 95% OF CPU TIME

        CALL FORCE (NATM, MX, MY, MZ, FX, FY, FZ, V, W,
     1              RX, RY, RZ, LIST, HEAD)

C       SECOND PART OF VELOCITY VERLET INTEGRATION

        CALL MOVEB (NATM, VX, VY, VZ, K, FX, FY, FZ, DT)

        E    = K+V
        EN   = E/DBLE(NATM)
        VN   = V/DBLE(NATM)
        KN   = K/DBLE(NATM)
        TEMP = 2.0D0*K/FREE
        PRES = DENS*TEMP+W
        PRES = PRES*SIGMA**3

C       OPTIONALLY PRINT INFORMATION

        IF ( MOD ( STEP, IPRINT ) .EQ. 0 ) THEN
             WRITE(*,'(1X,I8,5(2X,F10.6))')
     1                 STEP, EN, KN, VN, PRES, TEMP
             WRITE(NW1,'(1X,I8,5(2X,F10.6))')
     1                 STEP, EN, KN, VN, PRES, TEMP
          ENDIF


C       IF STILL EQUILIBRATING CALL SUBROUTINE TO SCALE VELOCITIES

        CALL SCALET (NATM, VX, VY, VZ, K, EQTEMP)

      ENDDO

C     START OF PRODUCTION MD CYCLE ** TIME THIS LOOP **

      CALL TIMER(TSTART)

      DO STEP=NEQUIL+1,NSTEP

        V = 0
        K = 0
        W = 0
        PRES = 0

C       CALL SUBROUTINE TO PERFORM FIRST
C       PART OF VELOCITY VERLET INTEGRATION
C       VECTOR OPERATIONS

        CALL MOVEA (NATM, RX, RY, RZ, VX, VY, VZ, FX, FY, FZ, DT)

C       CALL SUBROUTINE TO CREATE LINKED LIST AND PASS EDGES
C       BETWEEN PROCESSORS
C       THIS ROUTINE DOES ALL THE COMMUNICATIONS

        CALL MOVOUT (NATM, LIST, MX, MY, MZ, HEAD, RX, RY, RZ)

C       CALL SUBROUTINE TO CALCULATE FORCES AND POTENTIAL
C       USING LINKED
C       SCALAR ROUTINE MORE THAN 95% OF CPU TIME

        CALL FORCE (NATM, MX, MY, MZ, FX, FY, FZ, V, W,
     1              RX, RY, RZ, LIST, HEAD)

C       CALL SUBROUTINE TO PERFORM SECOND PART
C       OF VELOCITY VERLET INTEGRATION
C       VECTOR OPERATIONS

        CALL MOVEB (NATM, VX, VY, VZ, K, FX, FY, FZ, DT)

C       CALL SUBROUTINE TO PASS THERMODYNAMIC VARIABLES TO MASTER

C       CALCULATE AND ACCUMULATE GLOBAL DATA

        E    = K+V
        EN   = E/DBLE(NATM)
        VN   = V/DBLE(NATM)
        KN   = K/DBLE(NATM)
        TEMP = 2.0D0*K/FREE
        PRES = DENS*TEMP+W
        PRES = PRES*SIGMA**3

        ACE    = ACE  +EN
        ACV    = ACV  +VN
        ACK    = ACK  +KN
        ACP    = ACP  +PRES
        ACESQ  = ACESQ+EN*EN
        ACVSQ  = ACVSQ+VN*VN
        ACKSQ  = ACKSQ+KN*KN
        ACPSQ  = ACPSQ+PRES*PRES

C       OPTIONALLY PRINT INFORMATION

        IF ( MOD ( STEP, IPRINT ) .EQ. 0 ) THEN
             WRITE(*,'(1X,I8,5(2X,F10.6))')
     :                 STEP, EN, KN, VN, PRES, TEMP
             WRITE(NW1,'(1X,I8,5(2X,F10.6))')
     :                 STEP, EN, KN, VN, PRES, TEMP
          ENDIF

      ENDDO

      CALL TIMER(TEND)

      TIME = TEND-TSTART

C    *******************************************************************
C    ** ENDS THE LOOP OVER CYCLES                                     **
C    *******************************************************************

      WRITE(*,'(/'' ** MOLECULAR DYNAMICS ENDS  ** ''///)')

C     WRITE OUT FINAL AVERAGES

      NORM = DBLE(NSTEP-NEQUIL)
      AVE  = ACE/NORM
      AVK  = ACK/NORM
      AVV  = ACV/NORM
      AVP  = ACP/NORM

      ACESQ = (ACESQ/NORM)-AVE**2
      ACKSQ = (ACKSQ/NORM)-AVK**2
      ACVSQ = (ACVSQ/NORM)-AVV**2
      ACPSQ = (ACPSQ/NORM)-AVP**2

      IF (ACESQ.GT.0.0D0) FLE = DSQRT(ACESQ)
      IF (ACKSQ.GT.0.0D0) FLK = DSQRT(ACKSQ)
      IF (ACVSQ.GT.0.0D0) FLV = DSQRT(ACVSQ)
      IF (ACPSQ.GT.0.0D0) FLP = DSQRT(ACPSQ)

      AVT = AVK/1.5D0
      FLT = FLK/1.5D0

      WRITE(*,'('' AVERAGES'',5(2X,F10.5))')
     :             AVE, AVK, AVV, AVP, AVT
      WRITE(*,'('' FLUCTS  '',5(2X,F10.5))')
     :             FLE, FLK, FLV, FLP, FLT

      CALL FLOPS (NW1, NC, NSTEP - NEQUIL, RCUT, TIME,  MX, MY, MZ)

      CLOSE(NW1)

      STOP
10001 FORMAT(/1X,'TIMESTEP  ..ENERGY..  ..KINETIC..',
     :          '  ..POTENT..  .PRESSURE.  ..TEMPER..  ')

      END


      SUBROUTINE FCC (NATM, RX, RY, RZ, VX, VY, VZ, NC)
      IMPLICIT NONE
      INTEGER NATM
      INCLUDE 'md1.inc'
      DOUBLE PRECISION RX (NMAX), RY (NMAX), RZ (NMAX)
      DOUBLE PRECISION VX (NMAX), VY (NMAX), VZ (NMAX)
      INTEGER NC

C     *********************************************************************
C     ** ROUTINE TO CREATE AN FCC CONFIGURATION.                         **
C     **                                                                 **
C     ** THE INITIAL CONFIGURATION IS A FACE-CENTERED CUBIC ARRANGEMENT, **
C     ** CONSISTING OF NC * NC * NC UNIT CELLS, EACH CONTAINING 4 ATOMS. **
C     ** NOTE WELL THAT THESE UNIT CELLS HAVE NO CONNECTION WITH THE     **
C     ** CELLS USED DURING UPDATING!  THE UNIT CELL SIZE IS A * A * A.   **
C     ** THE FOUR ATOMS WITHIN THIS CELL ARE DISPLACED BY                **
C     **     (0, 0, 0), (A/2, 0, A/2), (A/2, A/2, 0), (0, A/2, A/2)      **
C     ** FROM THE CORNER.  THE `EPS' SHIFTS DON'T APPEAR TO SERVE ANY    **
C     ** USEFUL PURPOSE, BUT THEY WERE IN THE ORIGINAL CODE.             **
C     ** VELOCITIES ARE INITIALISED RANDOMLY, BUT ALTERNATE ATOMS HAVE   **
C     ** OPPOSITE VELOCITIES (SO THE TOTAL MOMENTUM IS ZERO).            **
C     *********************************************************************

        DOUBLE PRECISION EPS
        PARAMETER (EPS = 1.0D-6)
        DOUBLE PRECISION BASEX, BASEY, BASEZ
        PARAMETER (BASEX = -0.5, BASEY = -0.5, BASEZ = -0.5)

        INTEGER IS
        INTEGER LX, LY, LZ
        DOUBLE PRECISION A

        A = 1.0 / DBLE (NC)

C       SEED FOR RANDOM NUMBER GENERATOR

        IS = 2215

        NATM = 0
        DO LZ = 0, NC - 1
          DO LY = 0, NC - 1
            DO LX = 0, NC - 1
              NATM = NATM + 1
              RX (NATM) = BASEX + A * DBLE (LX) + EPS
              RY (NATM) = BASEY + A * DBLE (LY) + EPS
              RZ (NATM) = BASEZ + A * DBLE (LZ) + EPS
              CALL XSRAND (IS, VX (NATM))
              CALL XSRAND (IS, VY (NATM))
              CALL XSRAND (IS, VZ (NATM))

              NATM = NATM + 1
              RX (NATM) = BASEX + A * (DBLE (LX) + 0.5D0) - EPS
              RY (NATM) = BASEY + A *  DBLE (LY)          + EPS
              RZ (NATM) = BASEZ + A * (DBLE (LZ) + 0.5D0) - EPS
              VX (NATM) = 0.0D0 - VX (NATM - 1)
              VY (NATM) = 0.0D0 - VY (NATM - 1)
              VZ (NATM) = 0.0D0 - VZ (NATM - 1)

              NATM = NATM + 1
              RX (NATM) = BASEX + A * (DBLE (LX) + 0.5D0) - EPS
              RY (NATM) = BASEY + A * (DBLE (LY) + 0.5D0) - EPS
              RZ (NATM) = BASEZ + A *  DBLE (LZ)          + EPS
              CALL XSRAND (IS, VX (NATM))
              CALL XSRAND (IS, VY (NATM))
              CALL XSRAND (IS, VZ (NATM))

              NATM = NATM + 1
              RX (NATM) = BASEX + A *  DBLE (LX)          + EPS
              RY (NATM) = BASEY + A * (DBLE (LY) + 0.5D0) - EPS
              RZ (NATM) = BASEZ + A * (DBLE (LZ) + 0.5D0) - EPS
              VX (NATM) = 0.0D0 - VX (NATM - 1)
              VY (NATM) = 0.0D0 - VY (NATM - 1)
              VZ (NATM) = 0.0D0 - VZ (NATM - 1)
            ENDDO
          ENDDO
        ENDDO
      END


      SUBROUTINE XSRAND(IS,C)
      IMPLICIT NONE
      INTEGER IS
      DOUBLE PRECISION C

C     *********************************************************************
C     **   ROUTINE TO GENERATE RANDOM NUMBERS.                           **
C     **   THIS GENERATOR HAS NOT BEEN TESTED AND SHOULD NOT BE RELIED   **
C     **   UPON TO PRODUCE GENUINLY RANDOM NUMBERS!!                     **
C     *********************************************************************

        INTEGER IMULT, IMOD, IS1, IS2, ISS2
        DOUBLE PRECISION SCALE

        DATA IMULT /16807/, IMOD /2147483647/, SCALE /4.656612875D-10/

        IF (IS .LE. 0) IS = 1
        IS2  = MOD (IS, 32768)
        IS1  = (IS - IS2) / 32768
        ISS2 = IS * IMULT
        IS2  = MOD (ISS2, 32768)
        IS1  = MOD (IS1 * IMULT + (ISS2 - IS2) / 32768, 65536)
        IS   = MOD (IS1 * 32768 + IS2, IMOD)
        C = 2.0D0 * (SCALE * DBLE (IS) - 0.5D0) * 1.0D-6
      END


      SUBROUTINE MOVOUT (NATM, LIST, MX, MY, MZ, HEAD, RX, RY, RZ)
      IMPLICIT NONE
      INCLUDE 'md1.inc'
      INTEGER NATM
      INTEGER HEAD (MMX, MMY, MMZ), LIST (NMAX)
      DOUBLE PRECISION RX (NMAX), RY (NMAX), RZ (NMAX)

C     ************************************************************
C     ** BUILDS LINKED LISTS FROM VECTOR OF PARTICLES           **
C     ************************************************************

        INTEGER MX, MY, MZ

        DOUBLE PRECISION BASEX, BASEY, BASEZ
        PARAMETER (BASEX = -0.5, BASEY = -0.5, BASEZ = -0.5)

        INTEGER I, IX, IY, IZ

C       ZERO LINKED LISTS

        DO IX = 1, MX
          DO IY = 1, MY
            DO IZ = 1, MZ
              HEAD (IX, IY, IZ) = 0
            ENDDO
          ENDDO
        ENDDO

        DO I = 1, NATM
          IX = 1 + IDINT ((RX (I) - BASEX) * MX)
          IY = 1 + IDINT ((RY (I) - BASEY) * MY)
          IZ = 1 + IDINT ((RZ (I) - BASEZ) * MZ)
          LIST (I) = HEAD (IX, IY, IZ)
          HEAD (IX, IY, IZ) = I
        ENDDO
      END


      SUBROUTINE MAPS (MX, MY, MZ)
      IMPLICIT NONE
      INTEGER MX, MY,MZ

C     *******************************************************************
C     ** INITIALISE TABLES FOR ACCESSING NEIGHBOURS IN `FORCE'         **
C     *******************************************************************

        INCLUDE 'md1.inc'
        INTEGER SXAR (13), SYAR (13), SZAR (13)
        INTEGER JXAR (13, MMX), JYAR (13, MMY), JZAR (13, MMZ)
        DOUBLE PRECISION DSPXAR (13, MMX), DSPYAR (13, MMY),
     1                   DSPZAR (13, MMZ)
        COMMON /NEIGHB/ DSPXAR, DSPYAR, DSPZAR,
     1                  SXAR, SYAR, SZAR, JXAR, JYAR, JZAR

        INTEGER NABOR, S
        INTEGER IX, IY, IZ
        INTEGER JX, JY, JZ
        INTEGER SX, SY, SZ

        DO NABOR = 1, 13

C         `SX', `SY', `SZ' RUN OVER DISPLACEMENTS OF 13 NEIGHBOURS
C         WITH LEXICOGRAPHICALLY GREATER (IGNORING BC'S) COORDS.

          S = 13 + NABOR

          SX = MOD (S / 3, 3) - 1
          SY =      S / 9     - 1
          SZ = MOD (S,     3) - 1

          SXAR (NABOR) = SX
          SYAR (NABOR) = SY
          SZAR (NABOR) = SZ

          DO IX = 1, MX
            JX = 1 + MOD (IX + SX + MX - 1, MX)
            JXAR   (NABOR, IX) = JX
            DSPXAR (NABOR, IX) = (IX + SX - JX) / MX
          ENDDO

          DO IY = 1, MY
            JY = 1 + MOD (IY + SY + MY - 1, MY)
            JYAR   (NABOR, IY) = JY
            DSPYAR (NABOR, IY) = (IY + SY - JY) / MY
          ENDDO

          DO IZ = 1, MZ
            JZ = 1 + MOD (IZ + SZ + MZ - 1, MZ)
            JZAR   (NABOR, IZ) = JZ
            DSPZAR (NABOR, IZ) = (IZ + SZ - JZ) / MZ
          ENDDO
        ENDDO
      END


      SUBROUTINE MOVEA (NATM, RX, RY, RZ, VX, VY, VZ, FX, FY, FZ, DT)
      IMPLICIT NONE
      INTEGER NATM
      INCLUDE 'md1.inc'
      DOUBLE PRECISION RX (NMAX), RY (NMAX), RZ (NMAX)
      DOUBLE PRECISION VX (NMAX), VY (NMAX), VZ (NMAX)
      DOUBLE PRECISION FX (NMAX), FY (NMAX), FZ (NMAX)
      DOUBLE PRECISION DT

C     *******************************************************************
C     ** FIRST PART OF THE VELOCITY VERLET ALGORITHM                   **
C     *******************************************************************

        INTEGER I
        DOUBLE PRECISION DT2, DTSQ2

        DT2   = 0.5D0 * DT
        DTSQ2 = DT * DT2

        DO I = 1, NATM

C         INTEGRATE COORDINATES A COMPLETE TIMESTEP

          RX (I) = RX (I) + DT * VX (I) + DTSQ2 * FX (I)
          RY (I) = RY (I) + DT * VY (I) + DTSQ2 * FY (I)
          RZ (I) = RZ (I) + DT * VZ (I) + DTSQ2 * FZ (I)

          RX (I) = RX (I) - DINT (RX (I) + RX (I))
          RY (I) = RY (I) - DINT (RY (I) + RY (I))
          RZ (I) = RZ (I) - DINT (RZ (I) + RZ (I))

C         INTEGRATE VELOCITIES HALF A TIMESTEP

          VX (I) = VX (I) + DT2 * FX (I)
          VY (I) = VY (I) + DT2 * FY (I)
          VZ (I) = VZ (I) + DT2 * FZ (I)
        ENDDO
      END


      SUBROUTINE MOVEB (NATM, VX, VY, VZ, K, FX, FY, FZ, DT)
      IMPLICIT NONE
      INTEGER NATM
      INCLUDE 'md1.inc'
      DOUBLE PRECISION VX (NMAX), VY (NMAX), VZ (NMAX)
      DOUBLE PRECISION K
      DOUBLE PRECISION FX (NMAX), FY (NMAX), FZ (NMAX)
      DOUBLE PRECISION DT

C     *******************************************************************
C     ** SECOND PART OF VELOCITY VERLET ALGORITHM                      **
C     *******************************************************************

        DOUBLE PRECISION DT2
        INTEGER I

        DT2 = 0.5D0 * DT
        K   = 0.0D0
      
        DO I = 1, NATM

C         INTEGRATE VELOCITIES HALF A TIMESTEP

          VX (I) = VX (I) + DT2 * FX (I)
          VY (I) = VY (I) + DT2 * FY (I)
          VZ (I) = VZ (I) + DT2 * FZ (I)

C         SUM KINETIC ENERGY

          K = K + VX (I) * VX (I) + VY (I) * VY (I) + VZ (I) * VZ (I)
        ENDDO

        K = 0.5D0 * K
      END


      SUBROUTINE FORCE (NATM, MX, MY, MZ, FX, FY, FZ, V, W,
     1                  RX, RY, RZ, LIST, HEAD)
      IMPLICIT NONE
      INTEGER MX, MY, MZ
      INTEGER NATM
      INCLUDE 'md1.inc'
      DOUBLE PRECISION FX (NMAX), FY (NMAX), FZ (NMAX)
      DOUBLE PRECISION V, W
      DOUBLE PRECISION RX (NMAX), RY (NMAX), RZ (NMAX)
      INTEGER HEAD (MMX, MMY, MMZ), LIST (NMAX)

C     *******************************************************************
C     ** ROUTINE TO COMPUTE FORCES AND POTENTIAL USING A LINK LIST     **
C     *******************************************************************

        INTEGER SXAR (13), SYAR (13), SZAR (13)
        INTEGER JXAR (13, MMX), JYAR (13, MMY), JZAR (13, MMZ)
        DOUBLE PRECISION DSPXAR (13, MMX), DSPYAR (13, MMY),
     1                   DSPZAR (13, MMZ)
        COMMON /NEIGHB/ DSPXAR, DSPYAR, DSPZAR,
     1                  SXAR, SYAR, SZAR, JXAR, JYAR, JZAR

        DOUBLE PRECISION SIGMA, RCUT
        COMMON /POTNTL/ SIGMA, RCUT

        DOUBLE PRECISION RCUTSQ, SIGSQ
        DOUBLE PRECISION VRCUT, DVRCUT, DVRC12
        DOUBLE PRECISION RXJ, RYJ, RZJ, RXIJ, RYIJ, RZIJ
        DOUBLE PRECISION RIJ, RIJSQ, SR2, SR6
        DOUBLE PRECISION FXIJ, FYIJ, FZIJ
        DOUBLE PRECISION FIJ, VIJ, WIJ

        INTEGER IX, IY, IZ
        INTEGER I, J, ILOC, NABOR

        SIGSQ  = SIGMA**2
        RCUTSQ = RCUT**2

C       CALCULATE VALUES FOR SHIFTED FORCE POTENTIAL

        SR6    = (SIGMA / RCUT)**6
        VRCUT  = SR6 * (SR6 - 1.0D0)
        DVRCUT = -1.0D0 * SR6 * (SR6 - 0.5D0) / RCUT
        DVRC12 = 12.0D0 * DVRCUT

C       ZERO FORCES AND POTENTIAL

        DO I = 1, NATM
          FX(I) = 0.0D0
          FY(I) = 0.0D0
          FZ(I) = 0.0D0
        ENDDO

        V = 0.0D0
        W = 0.0D0

C       LOOP OVER ALL CELLS

        DO IX = 1, MX
          DO IY = 1, MY
            DO IZ = 1, MZ
        
C             INTERACTIONS INTERNAL TO CELL

              ILOC = HEAD (IX, IY, IZ)

              I = ILOC
              DO WHILE (I .GT. 0)

C               LOOP OVER ALL MOLECULES BELOW I IN THE CURRENT CELL

                J = LIST (I)
                DO WHILE (J .GT. 0)
                  RXIJ  = RX (I) - RX (J)
                  RYIJ  = RY (I) - RY (J)
                  RZIJ  = RZ (I) - RZ (J)
                  RIJSQ = RXIJ * RXIJ + RYIJ * RYIJ + RZIJ * RZIJ
                  IF (RIJSQ .LT. RCUTSQ) THEN
                    RIJ   = DSQRT (RIJSQ)
                    SR2   = SIGSQ / RIJSQ
                    SR6   = SR2 * SR2 * SR2
                    VIJ   = SR6 * (SR6 - 1.0) -
     1                      VRCUT - DVRC12 * (RIJ - RCUT)
                    WIJ   = SR6 * (SR6 - 0.5) + DVRCUT * RIJ
                    V     = V + VIJ
                    W     = W + WIJ
                    FIJ   = WIJ / RIJSQ
                    FXIJ  = FIJ * RXIJ
                    FYIJ  = FIJ * RYIJ
                    FZIJ  = FIJ * RZIJ
                    FX (I) = FX (I) + FXIJ
                    FY (I) = FY (I) + FYIJ
                    FZ (I) = FZ (I) + FZIJ

                    FX (J) = FX (J) - FXIJ
                    FY (J) = FY (J) - FYIJ
                    FZ (J) = FZ (J) - FZIJ
                  ENDIF
                  J = LIST (J)
                ENDDO

                I = LIST (I)
              ENDDO

C             LOOP OVER ALL NEIGHBOURING CELLS

              DO NABOR = 1, 13
                J = HEAD (JXAR (NABOR, IX), JYAR (NABOR, IY),
     1                    JZAR (NABOR, IZ))

                DO WHILE (J .NE. 0)
                  RXJ   = RX (J) + DSPXAR (NABOR, IX)
                  RYJ   = RY (J) + DSPYAR (NABOR, IY)
                  RZJ   = RZ (J) + DSPZAR (NABOR, IZ)

                  I = ILOC
                  DO WHILE (I .NE. 0)
                    RXIJ  = RX (I) - RXJ
                    RYIJ  = RY (I) - RYJ
                    RZIJ  = RZ (I) - RZJ
                    RIJSQ = RXIJ * RXIJ + RYIJ * RYIJ + RZIJ * RZIJ
                    IF (RIJSQ .LT. RCUTSQ) THEN
                      RIJ   = DSQRT (RIJSQ)
                      SR2   = SIGSQ / RIJSQ
                      SR6   = SR2 * SR2 * SR2
                      VIJ   = (SR6 * (SR6 - 1.0) -
     1                         VRCUT - DVRC12 * (RIJ - RCUT)) * 0.5
                      WIJ   = SR6 * (SR6 - 0.5) + DVRCUT * RIJ
                      FIJ   = WIJ / RIJSQ
                      WIJ   = 0.5D0 * WIJ
                      FXIJ  = FIJ * RXIJ
                      FYIJ  = FIJ * RYIJ
                      FZIJ  = FIJ * RZIJ

C                     SUM FORCES AND POTENTIAL

                      V      = V + VIJ
                      W      = W + WIJ
                      FX (I) = FX (I) + FXIJ
                      FY (I) = FY (I) + FYIJ
                      FZ (I) = FZ (I) + FZIJ

C                     SUM FORCES AND POTENTIAL

                      V      = V + VIJ
                      W      = W + WIJ
                      FX (J) = FX (J) - FXIJ
                      FY (J) = FY (J) - FYIJ
                      FZ (J) = FZ (J) - FZIJ
                    ENDIF
                    I = LIST (I)
                  ENDDO

                  J = LIST (J)
                ENDDO
              ENDDO
            ENDDO
          ENDDO
        ENDDO

C       INCORPORATE ENERGY FACTORS

        V =  4.0D0 * V
        W = 16.0D0 * W

        DO I = 1, NATM
          FX (I) = 48.0D0 * FX (I)
          FY (I) = 48.0D0 * FY (I)
          FZ (I) = 48.0D0 * FZ (I)
        ENDDO

      END


      SUBROUTINE SCALET (NATM, VX, VY, VZ, K, EQTEMP)
      IMPLICIT NONE
      INCLUDE 'md1.inc'
      INTEGER NATM
      DOUBLE PRECISION VX (NMAX), VY (NMAX), VZ (NMAX)
      DOUBLE PRECISION EQTEMP, K

C     *******************************************************************
C     ** SCALE VELOCITIES FOR EQUILIBRATION                            **
C     *******************************************************************

        DOUBLE PRECISION SCALEF
        INTEGER I

        SCALEF =
     1    DSQRT (EQTEMP / (2.0D0 * K / (3.0D0 * DBLE (NATM - 1))))

C       SCALE VELOCITIES

        DO I = 1, NATM
          VX (I) = VX (I) * SCALEF
          VY (I) = VY (I) * SCALEF
          VZ (I) = VZ (I) * SCALEF
        ENDDO

      END

