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    **           LOCAL LISTS.                                        **
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    ** MKLIST    CREATES LOCAL LISTS OF COORDINATES.
C    ** MKVECT    CREATES VECTOR OF FORCES.
C    ** SCALET    SCALES THE VELOCITIES TO ENABLE EQUILIBRATION TO A  **
C    **           PARTICULAR TEMPERATURE.                             **
C    **                                                               **
C    ** PRINCIPAL VARIABLES:                                          **
C    **                                                               **
C    ** INTEGER VARIABLES:                                            **
C    **                                                               **
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    ** NEQUIL                       NUMBER OF EQUILIBRATION TIMESTEPS**
C    **                              TO PERFORM.                      **
C    ** NSTEP                        NUMBER OF TIMESTEPS TO PERFORM.  **
C    **                                                               **
C    ** REAL VARIABLES:                                               **
C    **                                                               **
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    ** 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'

      INTEGER MAXCL
      DOUBLE PRECISION, ALLOCATABLE :: RX (:), RY (:), RZ (:)
      DOUBLE PRECISION, ALLOCATABLE :: VX (:), VY (:), VZ (:)
      DOUBLE PRECISION, ALLOCATABLE :: FX (:), FY (:), FZ (:)
      INTEGER, ALLOCATABLE :: NATMCL (:, :, :), HEAD (:, :, :),
     1        LIST (:, :, :, :), IDL (:, :, :, :)
      DOUBLE PRECISION, ALLOCATABLE ::
     1  FXL (:, :, :, :), FYL (:, :, :, :), FZL (:, :, :, :),
     2  RXL (:, :, :, :), RYL (:, :, :, :), RZL (:, :, :, :)

      DOUBLE PRECISION DT
      DOUBLE PRECISION DENS, DENLJ, EQTEMP
      DOUBLE PRECISION MASS, FREE, 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 MX,MY,MZ
      INTEGER NC,M
      INTEGER STEP, NSTEP, IPRINT, NORM, NATM
      INTEGER NEQUIL
      INTEGER I

      NATM = 0

      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 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     WRITE INPUT DATA

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

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

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

      NATM = 4 * NC * NC * NC
      ALLOCATE (RX (NATM), RY (NATM), RZ (NATM),
     1          VX (NATM), VY (NATM), VZ (NATM),
     2          FX (NATM), FY (NATM), FZ (NATM))

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

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

      CALL HEADER (NW1, NC, NATM)

      CALL PARAM (NW1, NSTEP, NEQUIL, IPRINT, 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     CALCULATE DIMENSIONS OF SIMULATION CELL STRUCTURE

      MX = IDINT (1.0D0 / RCUT)
      MY = IDINT (1.0D0 / RCUT)
      MZ = IDINT (1.0D0 / RCUT)

      ALLOCATE (NATMCL (MX, MY, MZ), HEAD (MX, MY, MZ),
     1          LIST (LMAX, MX, MY, MZ), IDL (LMAX, MX, MY, MZ))
      ALLOCATE (FXL (LMAX, MX, MY, MZ), FYL (LMAX, MX, MY, MZ),
     1          FZL (LMAX, MX, MY, MZ))
      ALLOCATE (RXL (LMAX, MX, MY, MZ), RYL (LMAX, MX, MY, MZ),
     1          RZL (LMAX, MX, MY, MZ))

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 MAP OF NEIGHBOURING CELLS THAT WE WISH TO CONSIDER
C     IN FORCE LOOP

      CALL MAPS (MX, MY, MZ)

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

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

C     START OF MD CYCLE

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

      DO STEP = 1, NSTEP

        IF (STEP .EQ. NEQUIL + 1) THEN

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

          CALL TIMER(TSTART)
        ENDIF

        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 LOCAL LISTS OF COORDINATES FROM VECTOR

        CALL MKLIST (MAXCL, MX * MY * MZ, MX, MY, MZ,
     1               NATMCL, IDL, RXL, RYL, RZL,
     2               NATM, RX, RY, RZ)

C       CALCULATE FORCES AND POTENTIAL.  MORE THAN 95% OF CPU TIME

        CALL FORCE (MX, MY, MZ, NATMCL, MAXCL, FXL, FYL, FZL, V, W,
     1              RXL, RYL, RZL)

C       CREATE VECTOR OF FORCES FROM LOCAL LISTS

        CALL MKVECT (NATM, FX, FY, FZ, MAXCL, MX * MY * MZ, MX, MY, MZ,
     1               NATMCL, IDL, FXL, FYL, FZL)

C       PERFORM SECOND PART OF VELOCITY VERLET INTEGRATION.

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

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

        IF (STEP .GT. NEQUIL) THEN
          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
        ENDIF

C       OPTIONALLY PRINT INFORMATION

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

        IF (STEP .LE. NEQUIL) THEN

C         IF STILL EQUILIBRATING, SCALE VELOCITIES

          CALL SCALET (NATM, VX, VY, VZ, K, EQTEMP)
        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
      DOUBLE PRECISION RX (NATM), RY (NATM), RZ (NATM)
      DOUBLE PRECISION VX (NATM), VY (NATM), VZ (NATM)
      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 I, LX, LY, LZ
        DOUBLE PRECISION A

        A = 1.0 / DBLE (NC)

C       SEED FOR RANDOM NUMBER GENERATOR

        IS = 2215

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

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

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

              I = I + 1
              RX (I) = BASEX + A *  DBLE (LX)          + EPS
              RY (I) = BASEY + A * (DBLE (LY) + 0.5D0) - EPS
              RZ (I) = BASEZ + A * (DBLE (LZ) + 0.5D0) - EPS
              VX (I) = 0.0D0 - VX (I - 1)
              VY (I) = 0.0D0 - VY (I - 1)
              VZ (I) = 0.0D0 - VZ (I - 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 MAPS (MX, MY, MZ)
      IMPLICIT NONE
      INTEGER MX, MY,MZ

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

        INCLUDE 'neighb.inc'

        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


C      SUBROUTINE MKLIST (MAXCL, MX, MY, MZ,
C     1                   NATMCL, IDL, RXL, RYL, RZL,
C     2                   NATM, RX, RY, RZ)
C      IMPLICIT NONE
C      INTEGER MAXCL
C      INTEGER MX, MY, MZ
C      INCLUDE 'md1.inc'
C      INTEGER NATMCL (MX, MY, MZ), IDL (LMAX, MX, MY, MZ)
C      DOUBLE PRECISION RXL (LMAX, MX, MY, MZ),
C     1                 RYL (LMAX, MX, MY, MZ),
C     2                 RZL (LMAX, MX, MY, MZ)
C      INTEGER NATM
C      DOUBLE PRECISION RX (NATM), RY (NATM), RZ (NATM)
C
CC     ************************************************************
CC     ** BUILDS LOCAL LISTS FROM VECTOR OF PARTICLES            **
CC     ************************************************************
C
C        DOUBLE PRECISION BASEX, BASEY, BASEZ
C        PARAMETER (BASEX = -0.5, BASEY = -0.5, BASEZ = -0.5)
C
C        INTEGER I, IX, IY, IZ, FREE
C
CC       ZERO LISTS
C
C        DO IX = 1, MX
C          DO IY = 1, MY
C            DO IZ = 1, MZ
C              NATMCL (IX, IY, IZ) = 0
C            ENDDO
C          ENDDO
C        ENDDO
C
C        MAXCL  = 0
C
CC       ALLOCATE ATOMS TO LISTS.
C
C        DO I = 1, NATM
C          IX = 1 + IDINT ((RX (I) - BASEX) * MX)
C          IY = 1 + IDINT ((RY (I) - BASEY) * MY)
C          IZ = 1 + IDINT ((RZ (I) - BASEZ) * MZ)
C
C          FREE = NATMCL (IX, IY, IZ) + 1
C
C          IDL (FREE, IX, IY, IZ) = I
C
C          RXL (FREE, IX, IY, IZ) = RX (I)
C          RYL (FREE, IX, IY, IZ) = RY (I)
C          RZL (FREE, IX, IY, IZ) = RZ (I)
C
C          NATMCL (IX, IY, IZ) = FREE
C          IF (FREE .GT. MAXCL) THEN
C            MAXCL = FREE
C          ENDIF
C        ENDDO
C      END

      SUBROUTINE MKLIST (MAXCL, NCELLS, MX, MY, MZ,
     1                   NATMCL, IDL, RXL, RYL, RZL,
     2                   NATM, RX, RY, RZ)
      IMPLICIT NONE
      INTEGER MAXCL, NCELLS
      INTEGER MX, MY, MZ
      INCLUDE 'md1.inc'
      INTEGER NATMCL (MX, MY, MZ), IDL (LMAX, MX, MY, MZ)
      DOUBLE PRECISION RXL (LMAX, MX, MY, MZ),
     1                 RYL (LMAX, MX, MY, MZ),
     2                 RZL (LMAX, MX, MY, MZ)
      INTEGER NATM
      DOUBLE PRECISION RX (NATM), RY (NATM), RZ (NATM)

C     ************************************************************
C     ** BUILDS LOCAL LISTS OF PARTICLE COORDINATES             **
C     ** (THIS IS A DIFFICULT COMBINING SCATTER OPERATION,      **
C     ** FROM THE ATOMIC COORDINATE VECTOR TO THE LOCAL LISTS)  **
C     ************************************************************

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

        INTEGER I
        INTEGER, DIMENSION (NATM) :: IX, IY, IZ
        INTEGER J

        INTEGER NATMCV (NCELLS), IDV (NCELLS)
        DOUBLE PRECISION RXV (NCELLS), RYV (NCELLS), RZV (NCELLS)

C       ZERO LISTS

        NATMCL = 0
        MAXCL  = 0

C       SCATTER ATOMIC IDS TO LOCAL LISTS, IDL.

        IX = 1 + IDINT ((RX - BASEX) * MX)
        IY = 1 + IDINT ((RY - BASEY) * MY)
        IZ = 1 + IDINT ((RZ - BASEZ) * MZ)

C        LOGICAL ORPHAN (NATM)
C        INTEGER ID (NATM)
C        LOGICAL UPDATE (MX, MY, MZ)
C        INTEGER SCATID (MX, MY, MZ)
C
C        ID = (/ (I = 1, NATM) /)
C        NATMCL = 0
C        ORPHAN = .TRUE.
C        DO WHILE (ANY (ORPHAN))
C          SCATID = -1
C          UPDATE = .FALSE.
C          SCATID = COPY_SCATTER (ID, SCATID,
C     1                           IX, IY, IZ, MASK = ORPHAN)
C          FORALL (I = 1 : NATM,
C     1            SCATID (IX (I), IY (I), IZ (I)) .EQ. I)
C            ORPHAN (I) = .FALSE.
C            UPDATE (IX (I), IY (I), IZ (I)) = .TRUE.
C          ENDFORALL
C
C          FORALL (X = 1 : MX, Y = 1 : MY, Z = 1 : MZ,
C     1            UPDATE (X, Y, Z))
C            NATMCL (X, Y, Z) = NATMCL (X, Y, Z) + 1
C            IDL (NATMCL (X, Y, Z), X, Y, Z) = SCATID (X, Y, Z)
C          ENDFORALL
C        ENDDO

C       **** FULL HPF VERSION ABOVE REPLACED BY FOLLOWING ``DO LOOP'' ****

        DO I = 1, NATM
          NATMCL (IX (I), IY (I), IZ (I)) =
     1      NATMCL (IX (I), IY (I), IZ (I)) + 1
          IDL (NATMCL (IX (I), IY (I), IZ (I)), IX (I), IY (I), IZ (I))
     1      = I
        ENDDO

        MAXCL = MAXVAL (NATMCL)

C       GIVEN LOCAL LIST OF ATOMIC IDS, GATHER ATOMIC COORDINATES.

        NATMCV = RESHAPE (NATMCL, (/ NCELLS /))
        DO I = 1, MAXCL
          IDV = RESHAPE (IDL (I, :, :, :), (/ NCELLS /))
C          WHERE (I .LE. NATMCV)
C            RXV = RX (IDV)
C            RYV = RY (IDV)
C            RZV = RZ (IDV)
C          ENDWHERE

C         ****  ``WHERE'' REPLACED BY FOLLOWING ``DO LOOP'', TO WORK ROUND
C                 APPARENT BUG IN NAG COMPILER ****

          DO J = 1, NCELLS
            IF (I .LE. NATMCV (J)) THEN
              RXV (J) = RX (IDV (J))
              RYV (J) = RY (IDV (J))
              RZV (J) = RZ (IDV (J))
            ENDIF
          ENDDO
          RXL (I, :, :, :) = RESHAPE (RXV, (/ MX, MY, MZ /))
          RYL (I, :, :, :) = RESHAPE (RYV, (/ MX, MY, MZ /))
          RZL (I, :, :, :) = RESHAPE (RZV, (/ MX, MY, MZ /))
        ENDDO
      END

      SUBROUTINE MKVECT (NATM, FX, FY, FZ, MAXCL, NCELLS, MX, MY, MZ,
     1                   NATMCL, IDL, FXL, FYL, FZL)
      IMPLICIT NONE
      INCLUDE 'md1.inc'
      INTEGER NATM
      DOUBLE PRECISION FX (NATM), FY (NATM), FZ (NATM)
      INTEGER NCELLS, MAXCL
      INTEGER MX, MY, MZ
      INTEGER NATMCL (MX, MY, MZ), IDL (LMAX, MX, MY, MZ)
      DOUBLE PRECISION FXL (LMAX, MX, MY, MZ),
     1                 FYL (LMAX, MX, MY, MZ),
     2                 FZL (LMAX, MX, MY, MZ)

C     ************************************************************
C     ** SCATTER LOCAL LISTS OF FORCES INTO ATOMIC FORCE VECTOR **
C     ************************************************************
        INTEGER I
        INTEGER NATMCV (NCELLS), IDV (NCELLS)
        DOUBLE PRECISION FXV (NCELLS), FYV (NCELLS), FZV (NCELLS)
        INTEGER J

        NATMCV = RESHAPE (NATMCL, (/ NCELLS /))
        DO I = 1, MAXCL
          IDV = RESHAPE (IDL (I, :, :, :), (/ NCELLS /))
          FXV = RESHAPE (FXL (I, :, :, :), (/ NCELLS /))
          FYV = RESHAPE (FYL (I, :, :, :), (/ NCELLS /))
          FZV = RESHAPE (FZL (I, :, :, :), (/ NCELLS /))
C          WHERE (I .LE. NATMCV)
C            FX (IDV) = FXV
C            FY (IDV) = FYV
C            FZ (IDV) = FZV
C          ENDWHERE

C         ****  ``WHERE'' REPLACED BY FOLLOWING ``DO LOOP'', TO WORK ROUND
C                 APPARENT BUG IN NAG COMPILER ****

          DO J = 1, NCELLS
            IF (I .LE. NATMCV (J)) THEN
              FX (IDV (J)) = FXV (J)
              FY (IDV (J)) = FYV (J)
              FZ (IDV (J)) = FZV (J)
            ENDIF
          ENDDO
        ENDDO
      END


      SUBROUTINE MOVEA (NATM, RX, RY, RZ, VX, VY, VZ, FX, FY, FZ, DT)
      IMPLICIT NONE
      INTEGER NATM
      DOUBLE PRECISION RX (NATM), RY (NATM), RZ (NATM)
      DOUBLE PRECISION VX (NATM), VY (NATM), VZ (NATM)
      DOUBLE PRECISION FX (NATM), FY (NATM), FZ (NATM)
      DOUBLE PRECISION DT

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

        DOUBLE PRECISION DT2, DTSQ2

        DT2   = 0.5D0 * DT
        DTSQ2 = DT * DT2

C       INTEGRATE COORDINATES A COMPLETE TIMESTEP

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

C       PUT PARTICLE BACK IN BOX (ASSUMING IT DOESN'T GO FURTHER
C       THAN 1 UNIT OUT ON NEGATIVE SIDE!)

        RX = MOD (1.5D0 + RX, 1.0D0) - 0.5D0
        RY = MOD (1.5D0 + RY, 1.0D0) - 0.5D0
        RZ = MOD (1.5D0 + RZ, 1.0D0) - 0.5D0

C       INTEGRATE VELOCITIES HALF A TIMESTEP

        VX = VX + DT2 * FX
        VY = VY + DT2 * FY
        VZ = VZ + DT2 * FZ
      END


      SUBROUTINE MOVEB (NATM, VX, VY, VZ, K, FX, FY, FZ, DT)
      IMPLICIT NONE
      INTEGER NATM
      DOUBLE PRECISION VX (NATM), VY (NATM), VZ (NATM)
      DOUBLE PRECISION K
      DOUBLE PRECISION FX (NATM), FY (NATM), FZ (NATM)
      DOUBLE PRECISION DT

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

        DOUBLE PRECISION DT2

        DT2 = 0.5D0 * DT
      
C       INTEGRATE VELOCITIES HALF A TIMESTEP

        VX = VX + DT2 * FX
        VY = VY + DT2 * FY
        VZ = VZ + DT2 * FZ

C       SUM KINETIC ENERGY

        K = 0.5D0 * SUM ((VX * VX) + (VY * VY) + (VZ * VZ))
      END


      SUBROUTINE FORCE (MX, MY, MZ, NATMCL, MAXCL, FXL, FYL, FZL, V, W,
     1                  RXL, RYL, RZL)
      IMPLICIT NONE
      INTEGER MX, MY, MZ
      INTEGER NATMCL (MX, MY, MZ)
      INTEGER MAXCL
      INCLUDE 'md1.inc'
      DOUBLE PRECISION FXL (LMAX, MX, MY, MZ),
     1                 FYL (LMAX, MX, MY, MZ),
     2                 FZL (LMAX, MX, MY, MZ)
      DOUBLE PRECISION V, W
      DOUBLE PRECISION RXL (LMAX, MX, MY, MZ),
     1                 RYL (LMAX, MX, MY, MZ),
     2                 RZL (LMAX, MX, MY, MZ)

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

        INCLUDE 'neighb.inc'

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

        DOUBLE PRECISION RCUTSQ, SIGSQ
        DOUBLE PRECISION VRCUT, DVRCUT, DVRC12
        DOUBLE PRECISION SRC6

        LOGICAL, DIMENSION (MX, MY, MZ) :: LIMASK, LJMASK, RCMASK
        DOUBLE PRECISION, DIMENSION (MX, MY, MZ) ::
     1    RXIJ, RYIJ, RZIJ, RIJSQ, RIJ, SR2, SR6, VIJ, WIJ,
     2    FIJ, FXIJ, FYIJ, FZIJ, RXJ, RYJ, RZJ,
     3    DSPX, DSPY, DSPZ
        INTEGER SX, SY, SZ
        INTEGER NATMCN (MX, MY, MZ)
        DOUBLE PRECISION, DIMENSION (LMAX, MX, MY, MZ) ::
     1    RXN, RYN, RZN, FXN, FYN, FZN

        INTEGER I, J, NABOR

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

C       CALCULATE VALUES FOR SHIFTED FORCE POTENTIAL

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

C       ZERO FORCES AND POTENTIAL

        FXL (1 : MAXCL, :, :, :) = 0.0D0
        FYL (1 : MAXCL, :, :, :) = 0.0D0
        FZL (1 : MAXCL, :, :, :) = 0.0D0

        V = 0.0D0
        W = 0.0D0

C       INTERACTIONS INTERNAL TO CELL

        DO I = 1, MAXCL
          LIMASK = I .LE. NATMCL

C         LOOP OVER ALL MOLECULES BELOW I IN THE CURRENT CELL

          DO J = I + 1, MAXCL
            LJMASK = LIMASK .AND. (J .LE. NATMCL)
            WHERE (LJMASK)
              RXIJ  = RXL (I, :, :, :) - RXL (J, :, :, :)
              RYIJ  = RYL (I, :, :, :) - RYL (J, :, :, :)
              RZIJ  = RZL (I, :, :, :) - RZL (J, :, :, :)
              RIJSQ = RXIJ**2 + RYIJ**2 + RZIJ**2
            ENDWHERE
            RCMASK = LJMASK .AND. (RIJSQ .LT. RCUTSQ)
            WHERE (RCMASK)
              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
              FIJ   = WIJ / RIJSQ
              FXIJ  = FIJ * RXIJ
              FYIJ  = FIJ * RYIJ
              FZIJ  = FIJ * RZIJ
              FXL (I, :, :, :) = FXL (I, :, :, :) + FXIJ
              FYL (I, :, :, :) = FYL (I, :, :, :) + FYIJ
              FZL (I, :, :, :) = FZL (I, :, :, :) + FZIJ

              FXL (J, :, :, :) = FXL (J, :, :, :) - FXIJ
              FYL (J, :, :, :) = FYL (J, :, :, :) - FYIJ
              FZL (J, :, :, :) = FZL (J, :, :, :) - FZIJ
            ENDWHERE
            V = V + SUM (VIJ, MASK = RCMASK)
            W = W + SUM (WIJ, MASK = RCMASK)
          ENDDO
        ENDDO

C       LOOP OVER ALL NEIGHBOURING CELLS

        DO NABOR = 1, 13

          SX = SXAR (NABOR)
          SY = SYAR (NABOR)
          SZ = SZAR (NABOR)
          NATMCN (:, :, :) = 
     1           CSHIFT (CSHIFT (CSHIFT (NATMCL (:, :, :),
     2                  SX, DIM = 1), SY, DIM = 2), SZ, DIM = 3)
          RXN (1 : MAXCL, :, :, :) =
     1           CSHIFT (CSHIFT (CSHIFT (RXL (1 : MAXCL, :, :, :),
     2                  SX, DIM = 2), SY, DIM = 3), SZ, DIM = 4)
          RYN (1 : MAXCL, :, :, :) =
     1           CSHIFT (CSHIFT (CSHIFT (RYL (1 : MAXCL, :, :, :),
     2                  SX, DIM = 2), SY, DIM = 3), SZ, DIM = 4)
          RZN (1 : MAXCL, :, :, :) =
     1           CSHIFT (CSHIFT (CSHIFT (RZL (1 : MAXCL, :, :, :),
     2                  SX, DIM = 2), SY, DIM = 3), SZ, DIM = 4)
          FXN (1 : MAXCL, :, :, :) = 0.0D0
          FYN (1 : MAXCL, :, :, :) = 0.0D0
          FZN (1 : MAXCL, :, :, :) = 0.0D0

          DSPX = SPREAD (SPREAD (DSPXAR (NABOR, 1 : MX), 2, MY), 3, MZ)
          DSPY = SPREAD (SPREAD (DSPYAR (NABOR, 1 : MY), 1, MX), 3, MZ)
          DSPZ = SPREAD (SPREAD (DSPZAR (NABOR, 1 : MZ), 1, MX), 2, MY)

          DO J = 1, MAXCL
            LJMASK = (J .LE. NATMCN)
            WHERE (LJMASK)
              RXJ = RXN (J, :, :, :) + DSPX
              RYJ = RYN (J, :, :, :) + DSPY
              RZJ = RZN (J, :, :, :) + DSPZ
            ENDWHERE

            DO I = 1, MAXCL
              LIMASK = LJMASK .AND. (I .LE. NATMCL)
              WHERE (LIMASK)
                RXIJ  = RXL (I, :, :, :) - RXJ
                RYIJ  = RYL (I, :, :, :) - RYJ
                RZIJ  = RZL (I, :, :, :) - RZJ
                RIJSQ = RXIJ**2 + RYIJ**2 + RZIJ**2
              ENDWHERE
              RCMASK = LIMASK .AND. (RIJSQ .LT. RCUTSQ)
              WHERE (RCMASK)
                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

                FXL (I, :, :, :) = FXL (I, :, :, :) + FXIJ
                FYL (I, :, :, :) = FYL (I, :, :, :) + FYIJ
                FZL (I, :, :, :) = FZL (I, :, :, :) + FZIJ

                FXN (J, :, :, :) = FXN (J, :, :, :) - FXIJ
                FYN (J, :, :, :) = FYN (J, :, :, :) - FYIJ
                FZN (J, :, :, :) = FZN (J, :, :, :) - FZIJ
              ENDWHERE
              V = V + 2 * SUM (VIJ, MASK = RCMASK)
              W = W + 2 * SUM (WIJ, MASK = RCMASK)
            ENDDO
          ENDDO

          FXL (1 : MAXCL, :, :, :) = FXL (1 : MAXCL, :, :, :) +
     1           CSHIFT (CSHIFT (CSHIFT (FXN (1 : MAXCL, :, :, :),
     2                  -SX, DIM = 2), -SY, DIM = 3), -SZ, DIM = 4)
          FYL (1 : MAXCL, :, :, :) = FYL (1 : MAXCL, :, :, :) +
     1           CSHIFT (CSHIFT (CSHIFT (FYN (1 : MAXCL, :, :, :),
     2                  -SX, DIM = 2), -SY, DIM = 3), -SZ, DIM = 4)
          FZL (1 : MAXCL, :, :, :) = FZL (1 : MAXCL, :, :, :) +
     1           CSHIFT (CSHIFT (CSHIFT (FZN (1 : MAXCL, :, :, :),
     2                  -SX, DIM = 2), -SY, DIM = 3), -SZ, DIM = 4)

        ENDDO

C       INCORPORATE ENERGY FACTORS

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

        FXL (1 : MAXCL, :, :, :) = 48.0D0 * FXL (1 : MAXCL, :, :, :)
        FYL (1 : MAXCL, :, :, :) = 48.0D0 * FYL (1 : MAXCL, :, :, :)
        FZL (1 : MAXCL, :, :, :) = 48.0D0 * FZL (1 : MAXCL, :, :, :)
      END


      SUBROUTINE SCALET (NATM, VX, VY, VZ, K, EQTEMP)
      IMPLICIT NONE
      INTEGER NATM
      DOUBLE PRECISION VX (NATM), VY (NATM), VZ (NATM)
      DOUBLE PRECISION EQTEMP, K

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

        DOUBLE PRECISION SCALEF

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

C       SCALE VELOCITIES

        VX = VX * SCALEF
        VY = VY * SCALEF
        VZ = VZ * SCALEF
      END

