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    ** COUNT     SENDS THERMODYNAMIC DATA BACK TO MASTER PROGRAM FOR **
C    **           ACCUMULATION.                                       **
C    ** FORCE     CALCULATES THE  FORCES AND POTENTIAL USING THE      **
C    **           LINKED LIST.                                        **
C    ** INPARS    RECEIVES SIMULATION DATA FROM MASTER PROGRAM AND    **
C    **           INITIALISES VARIOUS VARIBLES                        **
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. SENDS ATOMS THAT HAVE MOVED    **
C    **           OUT OF A PROCESSORS SIMULATION SPACE TO THE         **
C    **           RELEVANT NEW HOME PROCESSOR. PERFORMS THE           **
C    **           COMMUNICATIONS REQUIRED TO CREATE THE EDGE CELLS    **
C    **           NEEDED TO COMPLETE THE FORCE CALCULATIONS.          **
C    ** OUTCON    SENDS THE PROCESSORS CURRENT CONFIGURATIONAL        **
C    **           INFORMATION BACK TO THE MASTER PROGRAM.             **
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    ** PASSING CONVENTIONS:                                          **
C    ** EAST      DIRECTION OF INCREASING X COORDINATE                **
C    ** WEST      DIRECTION OF DECREASING X COORDINATE                **
C    ** NORTH     DIRECTION OF INCREASING Y COORDINATE                **
C    ** SOUTH     DIRECTION OF DECREASING Y COORDINATE                **
C    ** UP        DIRECTION OF INCREASING Z COORDINATE                **
C    ** DOWN      DIRECTION OF DECREASING Z COORDINATE                **
C    **                                                               **
C    *******************************************************************

      PROGRAM MD1

      IMPLICIT CHARACTER ( A-Z )

      COMMON/COORS/RX,RY,RZ
      COMMON/VELOC/VX,VY,VZ
      COMMON/STEPS/DT
      COMMON/LJPAR/SIGMA
      COMMON/ATOMS/NATM
      COMMON/CUTOFF/RCUT
      COMMON/TEMPER/EQTEMP
      COMMON/ENERGY/V,K,W
      COMMON/ESHIFT/VRCUT,DVRCUT,DVRC12
      COMMON/LSTDAT/LIST,HEAD,MAP
      COMMON/LSTDT2/SFX,SFY,SFZ,CELLIX,CELLIY,CELLIZ,MX,MY,MZ
      COMMON/CONTRL/NSTEP,NEQUIL,ISAVE
      COMMON/FORCES/FX,FY,FZ
      COMMON/PROCS/NPROSX, NPROSY, NPROSZ

      INTEGER     NMAX, MMX, MMY, MMZ, NCELL, MAPSIZ
      INTEGER     NW1

      INCLUDE 'md1.inc'

      PARAMETER(NW1 = 11)

      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        RCUT, DT, SIGMA
      DOUBLE PRECISION        VRCUT,DVRCUT,DVRC12
      DOUBLE PRECISION        SFX,SFY,SFZ,CELLIX,CELLIY,CELLIZ
      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(NCELL), LIST(NMAX), MAP(MAPSIZ),MX,MY,MZ
      INTEGER     NC,M
      INTEGER     STEP, NSTEP, ISAVE, IPRINT, NORM, NATM
      INTEGER     NEQUIL, ISCALE, NPROSX, NPROSY, NPROSZ

      NATM = 0

      PI = ACOS(-1.D0)

C    ** READ INPUT DATA **

      WRITE(*,'(1H1,'' **** PROGRAM MDLJ ****                   '')')
      WRITE(*,'('' MOLECULAR DYNAMICS SIMULATION                '')')
      WRITE(*,'('' WITH LINKED LIST                             '')')
C      WRITE(*,'('' ENTER NUMBER OF PROCESSORS IN X DIRECTION    '')')
C      READ (*,*) NPROSX
C      WRITE(*,'('' ENTER NUMBER OF PROCESSORS IN Y DIRECTION    '')')
C      READ (*,*) NPROSY
C      WRITE(*,'('' ENTER NUMBER OF PROCESSORS IN Z DIRECTION    '')')
C      READ (*,*) NPROSZ
C  ENTER NUMBER OF TIME STEPS
      NSTEP = 600
C  ENTER NUMBER OF EQUILIBRATION TIME STEPS
      NEQUIL = 500
C  ENTER NUMBER OF STEPS BETWEEN DATA SAVES
      ISAVE = 2*NSTEP
C  ENTER NUMBER OF STEPS BETWEEN OUTPUT
      IPRINT = 100
C  ENTER THE FOLLOWING IN LENNARD-JONES UNITS
C  ENTER THE TEMPERATURE
      EQTEMP = 0.8D0
C  ENTER THE DENSITY
      DENS = 0.75D0
C  ENTER THE POTENTIAL CUTOFF DISTANCE
      RCUT = 2.0D0**(1.0D0/6.0D0)
C  ENTER 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 SUBLATTICES TO USE IN ALL THREE
C**   COORDINATE DIRECTIONS
C**   AT THIS POINT THE NUMBER OF ATOMS WOULD BE READ FROM A FILE
C**   IF SUBROUTINE INCON OR INCONB WERE BEING USED

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

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

      CALL FCC(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**   CALL SUBROUTINE TO RECEIVE SIMULATION PARAMETERS FROM MASTER

      CALL INPARS

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

C**   CALL SUBROUTINE TO 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 1000 STEP=1,MIN(NSTEP,NEQUIL)
        WRITE (6, *) 'STEP =', STEP
        V = 0
        K = 0
        W = 0
        PRES = 0

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

        CALL MOVEA

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

        CALL MOVOUT

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

        CALL FORCE

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

        CALL MOVEB

        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))')
     :                 STEP, EN, KN, VN, PRES, TEMP
             WRITE(NW1,'(1X,I8,5(2X,F10.6))')
     :                 STEP, EN, KN, VN, PRES, TEMP
          ENDIF


C**     IF STILL EQUILIBRATING CALL SUBROUTINE TO SCALE VELOCITIES

        CALL SCALET

1000    CONTINUE

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

      CALL TIMER(TSTART)

      DO 2000 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

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

        CALL MOVOUT

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

        CALL FORCE

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

        CALL MOVEB

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

2000  CONTINUE

      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)

      CLOSE(NW1)

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

      END




      SUBROUTINE INPARS

      IMPLICIT CHARACTER ( A-Z )

      COMMON/ATOMS/NATM
      COMMON/STEPS/DT
      COMMON/CUTOFF/RCUT
      COMMON/LJPAR/SIGMA
      COMMON/TEMPER/EQTEMP
      COMMON/ESHIFT/VRCUT,DVRCUT,DVRC12
      COMMON/LSTDAT/LIST,HEAD,MAP
      COMMON/LSTDT2/SFX,SFY,SFZ,CELLIX,CELLIY,CELLIZ,MX,MY,MZ
      COMMON/CONTRL/NSTEP,NEQUIL,ISAVE

C    *******************************************************************
C    ** SUBROUTINE TO READ IN L-J PARAMETERS, INITIALISE COUNTERS ETC **
C    *******************************************************************

      INTEGER     NMAX, MMX, MMY, MMZ, NCELL, MAPSIZ

      INCLUDE 'md1.inc'

      DOUBLE PRECISION         DT, RCUT, SIGMA, EQTEMP
      DOUBLE PRECISION         SR6,VRCUT,DVRCUT,DVRC12
      DOUBLE PRECISION         SFX,SFY,SFZ,CELLIX,CELLIY,CELLIZ

      INTEGER      HEAD(NCELL), LIST(NMAX), MAP(MAPSIZ),MX,MY,MZ
      INTEGER      NATM
      INTEGER      NSTEP,NEQUIL,ISAVE

C**   CALCULATE DIMENSIONS OF PROCESSORS LINK-CELL STRUCTURE

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

C**   CALCULATE SCALE FACTORS TO ENSURE EDGE CELLS ARE INCLUDED
C**   IN LINKED LIST

      SFX = (DBLE(MX)-2.D0)/DBLE(MX)
      SFY = (DBLE(MY)-2.D0)/DBLE(MY)
      SFZ = (DBLE(MZ)-2.D0)/DBLE(MZ)

C**   REAL VALUES OF NUMBER OF CELLS IN EACH COORDINATE DIRECTION

      CELLIX = DBLE(MX)
      CELLIY = DBLE(MY)
      CELLIZ = DBLE(MZ)

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

      CALL MAPS

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

      END



      SUBROUTINE FCC(NC)

C*********************************************************************
C**   ROUTINE TO CREATE AN FCC CONFIGURATION 
C*********************************************************************

      IMPLICIT CHARACTER ( A - Z )

      INTEGER     NMAX, MMX, MMY, MMZ, NCELL, MAPSIZ

      INCLUDE 'md1.inc'

      COMMON/ATOMS/NATM
      COMMON/COORS/RX,RY,RZ
      COMMON/VELOC/VX,VY,VZ

      DOUBLE PRECISION        RX(NMAX), RY(NMAX), RZ(NMAX)
      DOUBLE PRECISION        VX(NMAX), VY(NMAX), VZ(NMAX)

      INTEGER     IS,NC
      INTEGER     NATM
      INTEGER     LX,LY,LZ

      NATM = 0
C**   SEED FOR RANDOM NUMBER GENERATOR

      IS = 2215

C**   TRIPLE LOOP TO CREATE LATTICE

      DO 40 LZ=0,NC-1
        DO 50 LY=0,NC-1
          DO 60 LX=0,NC-1

C**   COORDINATES S.T. THEY WITHIN SIMULATION BOX
C**   RANDOM VELICITIES

            NATM = NATM+1
            RX(NATM) = DBLE(LX)/DBLE(NC)+1.0D-6-0.5D0
            RY(NATM) = DBLE(LY)/DBLE(NC)+1.0D-6-0.5D0
            RZ(NATM) = DBLE(LZ)/DBLE(NC)+1.0D-6-0.5D0

            CALL XSRAND(IS,VX(NATM))
            CALL XSRAND(IS,VY(NATM))
            CALL XSRAND(IS,VZ(NATM))

            NATM = NATM+1
            RX(NATM) = (DBLE(LX)+0.5D0)/DBLE(NC)-1.0D-6-0.5D0
            RY(NATM) = DBLE(LY)/DBLE(NC)+1.0D-6-0.5D0
            RZ(NATM) = (DBLE(LZ)+0.5D0)/DBLE(NC)-1.0D-6-0.5D0
            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) = (DBLE(LX)+0.5D0)/DBLE(NC)-1.0D-6-0.5D0
            RY(NATM) = (DBLE(LY)+0.5D0)/DBLE(NC)-1.0D-6-0.5D0
            RZ(NATM) = DBLE(LZ)/DBLE(NC)+1.0D-6-0.5D0
            CALL XSRAND(IS,VX(NATM))
            CALL XSRAND(IS,VY(NATM))
            CALL XSRAND(IS,VZ(NATM))

            NATM = NATM+1
            RX(NATM) = DBLE(LX)/DBLE(NC)+1.0D-6-0.5D0
            RY(NATM) = (DBLE(LY)+0.5D0)/DBLE(NC)-1.0D-6-0.5D0
            RZ(NATM) = (DBLE(LZ)+0.5D0)/DBLE(NC)-1.0D-6-0.5D0
            VX(NATM) = 0.0D0-VX(NATM-1)
            VY(NATM) = 0.0D0-VY(NATM-1)
            VZ(NATM) = 0.0D0-VZ(NATM-1)

 60         CONTINUE
 50       CONTINUE
 40     CONTINUE

      END



      SUBROUTINE XSRAND(IS,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*********************************************************************

      DOUBLE PRECISION SCALE, C

      INTEGER IS,IMULT,IMOD,IS1,IS2,ISS2

      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

C************************************************************
C**                                                        **
C** PASSES EDGE OF SIMULATION BOX TO NEIGHBOURS AND ADDS   **
C** THOSE PASSED TO IT TO THE ARRAYS OF COORDINATES.       **
C** DURING THIS PROCESS ATOMS THAT HAVE MOVED INTO OTHER   **
C** PROCESSOR'S SIMULATION SPACE ARE PASSED TO THEIR       **
C** NEW HOME PROCESSOR.                                    **
C**                                                        **
C************************************************************

      IMPLICIT CHARACTER ( A-Z )

      COMMON/COORS/RX,RY,RZ
      COMMON/VELOC/VX,VY,VZ
      COMMON/ATOMS/NATM
      COMMON/LSTDAT/LIST,HEAD,MAP
      COMMON/LSTDT2/SFX,SFY,SFZ,CELLIX,CELLIY,CELLIZ,MX,MY,MZ

      INTEGER     NMAX, MMX, MMY, MMZ,  NCELL, MAPSIZ

      INCLUDE 'md1.inc'

      DOUBLE PRECISION        RX(NMAX), RY(NMAX), RZ(NMAX)
      DOUBLE PRECISION        VX(NMAX), VY(NMAX), VZ(NMAX)
      DOUBLE PRECISION        SFX,SFY,SFZ,CELLIX,CELLIY,CELLIZ

      INTEGER     HEAD(NCELL), LIST(NMAX), MAP(MAPSIZ), MX, MY,MZ
      INTEGER     NATM, I, J, K, ICELL, NEDGE

C**   ZERO LINKED LIST

      DO 5 ICELL = 1, MX*MY*MZ
        HEAD(ICELL) = 0
 5      CONTINUE

      DO 10 I = 1, NATM
        ICELL = 1+IDINT((RX(I)*SFX+0.5D0)*CELLIX)
     &           +IDINT((RY(I)*SFY+0.5D0)*CELLIY)*MX
     &           +IDINT((RZ(I)*SFZ+0.5D0)*CELLIZ)*MX*MY
        LIST(I) = HEAD(ICELL)
        IF (I .EQ. 48) THEN
          WRITE (6, *) 'I =', I
          WRITE (6, *) 'RX, RY, RZ (I) =', RX (I), RY (I), RZ (I)
          WRITE (6, *) 'ICELL =', ICELL
        ENDIF

        HEAD(ICELL) = I
 10     CONTINUE

      NEDGE = NMAX+1

C**   CREATE WEST FACE OF BOX

      DO 20 K = 0,MZ-1
        DO 21 I= MX-1,MX*MY-1,MX
          J=HEAD(I+K*MX*MY)
 22       IF (J.NE.0) THEN
          NEDGE = NEDGE-1
          RX(NEDGE) = RX(J)-1.0D0
          RY(NEDGE) = RY(J)
          RZ(NEDGE) = RZ(J)
          ICELL = 1+IDINT((RY(NEDGE)*SFY+0.5D0)*CELLIY)*MX
     &             +IDINT((RZ(NEDGE)*SFZ+0.5D0)*CELLIZ)*MX*MY
          LIST(NEDGE) = HEAD(ICELL)
          HEAD(ICELL) = NEDGE

          J = LIST(J)
          GOTO 22
          ENDIF
 21     CONTINUE
 20   CONTINUE

C**   CREATE EAST FACE OF BOX

      DO 30 K = 0,MZ-1
        DO 31 I = 2, MX*(MY-1)+2,MX
          J=HEAD(I+K*MX*MY)
 32       IF (J.NE.0) THEN
          NEDGE = NEDGE-1
          RX(NEDGE) = RX(J)+1.0D0
          RY(NEDGE) = RY(J)
          RZ(NEDGE) = RZ(J)
          ICELL = MX+IDINT((RY(NEDGE)*SFY+0.5D0)*CELLIY)*MX
     &              +IDINT((RZ(NEDGE)*SFZ+0.5D0)*CELLIZ)*MX*MY
          LIST(NEDGE) = HEAD(ICELL)
          HEAD(ICELL) = NEDGE

          J = LIST(J)
          GOTO 32
          ENDIF
 31     CONTINUE
 30   CONTINUE

C**   CREATE SOUTH FACE OF BOX

      DO 100 K = 0,MZ-1
        DO 101 I= MX*(MY-2)+1,MX*(MY-1)
          J=HEAD(I+K*MX*MY)
 102      IF (J.NE.0) THEN
          NEDGE = NEDGE-1
          RX(NEDGE) = RX(J)
          RY(NEDGE) = RY(J)-1.0D0
          RZ(NEDGE) = RZ(J)
          ICELL = 1+IDINT((RX(NEDGE)*SFX+0.5D0)*CELLIX)
     &             +IDINT((RZ(NEDGE)*SFZ+0.5D0)*CELLIZ)*MX*MY
          LIST(NEDGE) = HEAD(ICELL)
          HEAD(ICELL) = NEDGE

          J = LIST(J)
          GOTO 102
          ENDIF
 101    CONTINUE
 100  CONTINUE

C**   CREATE NORTH FACE OF BOX

      DO 130 K = 0,MZ-1
        DO 131 I = MX+1,2*MX
          J=HEAD(I+K*MX*MY)
 132      IF (J.NE.0) THEN
          NEDGE = NEDGE-1
          RX(NEDGE) = RX(J)
          RY(NEDGE) = RY(J)+1.0D0
          RZ(NEDGE) = RZ(J)
          ICELL = 1+IDINT((RX(NEDGE)*SFX+0.5D0)*CELLIX)+MX*(MY-1)
     &             +IDINT((RZ(NEDGE)*SFZ+0.5D0)*CELLIZ)*MX*MY
          LIST(NEDGE) = HEAD(ICELL)
          HEAD(ICELL) = NEDGE

          J = LIST(J)
          GOTO 132
          ENDIF
 131    CONTINUE
 130  CONTINUE

C**   CREATE DOWN FACE OF BOX

      DO 220 K = MX*MY*(MZ-2)+1,MX*MY*(MZ-1)
        J=HEAD(K)
 221    IF (J.NE.0) THEN
        NEDGE = NEDGE-1
        RX(NEDGE) = RX(J)
        RY(NEDGE) = RY(J)
        RZ(NEDGE) = RZ(J)-1.0D0
        ICELL = 1+IDINT((RX(NEDGE)*SFX+0.5D0)*CELLIX)
     &           +IDINT((RY(NEDGE)*SFY+0.5D0)*CELLIY)*MX
        LIST(NEDGE) = HEAD(ICELL)
        HEAD(ICELL) = NEDGE

        J = LIST(J)
        GOTO 221
        ENDIF
 220  CONTINUE


C**   CREATE UP FACE OF BOX

      DO 270 K = MX*MY+1,2*MX*MY
        J=HEAD(K)
 271    IF (J.NE.0) THEN
        NEDGE = NEDGE-1
        RX(NEDGE) = RX(J)
        RY(NEDGE) = RY(J)
        RZ(NEDGE) = RZ(J)+1.0D0
        ICELL = 1+IDINT((RX(NEDGE)*SFX+0.5D0)*CELLIX)
     &           +IDINT((RY(NEDGE)*SFY+0.5D0)*CELLIY)*MX
     &           +(MZ-1)*MX*MY
        LIST(NEDGE) = HEAD(ICELL)
        HEAD(ICELL) = NEDGE

        J = LIST(J)
        GOTO 271
        ENDIF
 270  CONTINUE

C**   TEST TO SEE IF WE HAVE TOO MANY EDGE ATOMS FOR ARRAYS

      IF (NEDGE.LT.NATM) WRITE(*,*) ' LINKED LIST RECURSIVE '

      END




      SUBROUTINE MAPS

      IMPLICIT CHARACTER ( A-Z )

      COMMON/LSTDAT/LIST,HEAD,MAP
      COMMON/LSTDT2/SFX,SFY,SFZ,CELLIX,CELLIY,CELLIZ,MX,MY,MZ

C    *******************************************************************
C    ** ROUTINE TO SET UP A LIST OF NEIGHBOURING CELLS                **
C    *******************************************************************

      INTEGER     NMAX, MMX, MMY, MMZ, NCELL, MAPSIZ

      INCLUDE 'md1.inc'

      DOUBLE PRECISION        SFX,SFY,SFZ,CELLIX,CELLIY,CELLIZ

      INTEGER     LIST(NMAX), HEAD(NCELL), MAP(MAPSIZ), MX, MY,MZ
      INTEGER     IX, IY, IZ, IMAP, ICELL

C    ** STATEMENT FUNCTION TO GIVE CELL INDEX **
      ICELL(IX,IY,IZ) = 1+MOD(IX-1+MX,MX)+MOD(IY-1+MY,MY)*MX
     &                   +MOD(IZ-1+MZ,MZ)*MX*MY

C    ** FIND HALF THE NEAREST NEIGHBOURS OF EACH CELL **
      DO 40 IX = 1, MX
        DO 41 IY = 1, MY
          DO 42 IZ = 1, MZ
            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)
42        CONTINUE
41      CONTINUE
40    CONTINUE

      END



      SUBROUTINE MOVEA

      IMPLICIT CHARACTER ( A-Z )

      COMMON/COORS/RX,RY,RZ
      COMMON/VELOC/VX,VY,VZ
      COMMON/STEPS/DT
      COMMON/ATOMS/NATM
      COMMON/FORCES/FX,FY,FZ

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

      INTEGER     NMAX, MMX, MMY, MMZ, NCELL, MAPSIZ

      INCLUDE 'md1.inc'

      DOUBLE PRECISION        DT, DT2, DTSQ2
      DOUBLE PRECISION        RX(NMAX), RY(NMAX), RZ(NMAX)
      DOUBLE PRECISION        VX(NMAX), VY(NMAX), VZ(NMAX)
      DOUBLE PRECISION        FX(NMAX), FY(NMAX), FZ(NMAX)

      INTEGER     I, NATM

      DT2   = 0.5D0*DT
      DTSQ2 = DT*DT2

      DO 100 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)
100   CONTINUE

      END




      SUBROUTINE MOVEB

      IMPLICIT CHARACTER ( A-Z )

      COMMON/VELOC/VX,VY,VZ
      COMMON/STEPS/DT
      COMMON/ATOMS/NATM
      COMMON/ENERGY/V,K,W
      COMMON/FORCES/FX,FY,FZ

      INTEGER     NMAX, MMX, MMY, MMZ, NCELL, MAPSIZ

      INCLUDE 'md1.inc'

      DOUBLE PRECISION        DT, DT2
      DOUBLE PRECISION        V,K,W
      DOUBLE PRECISION        VX(NMAX), VY(NMAX), VZ(NMAX)
      DOUBLE PRECISION        FX(NMAX), FY(NMAX), FZ(NMAX)

      INTEGER     I, NATM

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

      DT2   = 0.5D0*DT
      K   = 0.0D0
    
      DO 100 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)

100   CONTINUE

      K = 0.5D0*K

      END



      SUBROUTINE FORCE

      IMPLICIT CHARACTER ( A-Z )

      COMMON/COORS/RX,RY,RZ
      COMMON/LSTDAT/LIST,HEAD,MAP
      COMMON/LSTDT2/SFX,SFY,SFZ,CELLIX,CELLIY,CELLIZ,MX,MY,MZ
      COMMON/ENERGY/V,K,W
      COMMON/ATOMS/NATM
      COMMON/LJPAR/SIGMA
      COMMON/ESHIFT/VRCUT,DVRCUT,DVRC12
      COMMON/CUTOFF/RCUT
      COMMON/FORCES/FX,FY,FZ

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

      INTEGER     NMAX, MMX, MMY, MMZ, NCELL, MAPSIZ

      INCLUDE 'md1.inc'

      DOUBLE PRECISION  RX(NMAX), RY(NMAX), RZ(NMAX)
      DOUBLE PRECISION  FX(NMAX), FY(NMAX), FZ(NMAX)
      DOUBLE PRECISION  RCUT, SIGMA
      DOUBLE PRECISION  VRCUT,DVRCUT,DVRC12
      DOUBLE PRECISION  V,K,W
      DOUBLE PRECISION  RXI, RYI, RZI, FXIJ, FYIJ, FZIJ, FIJ, RCUTSQ
      DOUBLE PRECISION  SIGSQ, FXI, FYI, FZI, SR2, SR6, VIJ, WIJ
      DOUBLE PRECISION  RIJ, RIJSQ, RXIJ, RYIJ, RZIJ
      DOUBLE PRECISION  SFX,SFY,SFZ,CELLIX,CELLIY,CELLIZ

      LOGICAL     EDGEI, EDGEJ

      INTEGER     HEAD(NCELL), LIST(NMAX), MAP(MAPSIZ), MX, MY,MZ
      INTEGER     ICELL, JCELL0, JCELL, I, J, NABOR, NATM

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

C    ** ZERO FORCS AND POTENTIAL **
      DO 10 I = 1, NATM
        FX(I) = 0.0D0
        FY(I) = 0.0D0
        FZ(I) = 0.0D0
10      CONTINUE

      V = 0.0D0
      W = 0.0D0

C    ** LOOP OVER ALL CELLS **

      DO 5000 ICELL = 1 , MX*MY*MZ

        I = HEAD(ICELL)

C**     TEST TO SEE IF CURRENT CELL IS AN EDGE CELL

        EDGEI=(I.GT.NATM)

C       ** LOOP OVER ALL MOLECULES IN THE CELL **

1000    IF (I .GT. 0) THEN
          RXI = RX(I)
          RYI = RY(I)
          RZI = RZ(I)

C**       IF AN EDGE CELL OMIT FORCE CALCULATIONS BETWEEN ATOMS WITHIN
C**       CURRENT CELL

          IF (.NOT.EDGEI) THEN
          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)
            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.0D0)-VRCUT-DVRC12*(RIJ-RCUT)
              WIJ   = SR6*(SR6-0.5D0)+DVRCUT*RIJ
              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
          ENDIF

C             ** LOOP OVER ALL MOLECULES IN NEIGHBOURING CELLS **

        JCELL0 = 13*(ICELL-1)

        DO 4000 NABOR = 1, 13

          JCELL = MAP(JCELL0+NABOR)
          J = HEAD(JCELL)

C**     TEST TO SEE IF NEIGHBOURING CELL IS AN EDGE CELL

          EDGEJ=(J.GT.NATM)

C**       IF BOTH EDGE CELLS OMIT FORCE CALCULATIONS BETWEEN PAIRS
C**       OF EDGE CELLS

          IF (EDGEI.AND.EDGEJ) GOTO 4000
3000      IF (J.NE.0) THEN
            RXIJ  = RXI-RX(J)
            RYIJ  = RYI-RY(J)
            RZIJ  = RZI-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.0D0)-VRCUT-DVRC12*(RIJ-RCUT))*0.5D0
              WIJ   = SR6*(SR6-0.5D0)+DVRCUT*RIJ
              FIJ   = WIJ/RIJSQ
              WIJ   = 0.5D0*WIJ
              FXIJ  = FIJ*RXIJ
              FYIJ  = FIJ*RYIJ
              FZIJ  = FIJ*RZIJ

C**   IF CURRENT CELL IS NOT AN EDGE CELL SUM FORCES AND POTENTIAL

              IF (.NOT.EDGEI) THEN
                V     = V+VIJ
                W     = W+WIJ
                FXI   = FXI+FXIJ
                FYI   = FYI+FYIJ
                FZI   = FZI+FZIJ
                ENDIF

C**   IF NEIGBOURING CELL IS NOT AN EDGE CELL SUM FORCES AND POTENTIAL

              IF (.NOT.EDGEJ) THEN
                V     = V+VIJ
                W     = W+WIJ
                FX(J) = FX(J)-FXIJ
                FY(J) = FY(J)-FYIJ
                FZ(J) = FZ(J)-FZIJ
                ENDIF
              ENDIF
            J = LIST(J)
            GO TO 3000
            ENDIF
4000      CONTINUE

C**   IF CURRENT CELL IS NOT AN EDGE CELL RETURN SUMMED FORCES TO ARRAYS

          IF(.NOT.EDGEI)THEN
            FX(I) = FXI
            FY(I) = FYI
            FZ(I) = FZI
          ENDIF


          I = LIST(I)
          GO TO 1000
          ENDIF
5000    CONTINUE

C    ** INCORPORATE ENERGY FACTORS **

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

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

      END




      SUBROUTINE SCALET

      IMPLICIT CHARACTER ( A-Z )

      COMMON/VELOC/VX,VY,VZ
      COMMON/TEMPER/EQTEMP
      COMMON/ENERGY/V,K,W
      COMMON/ATOMS/NATM

      INTEGER     NMAX, MMX, MMY, MMZ, NCELL, MAPSIZ

      INCLUDE 'md1.inc'

      DOUBLE PRECISION         VX(NMAX), VY(NMAX), VZ(NMAX)
      DOUBLE PRECISION         EQTEMP,V,K,W
      DOUBLE PRECISION         SCALEF

      INTEGER      NATM, I

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

C**   SCALE VELOCITIES

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

      END
