C *******************************************************************
C    ** WORKER PROGRAM                                                **
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, WHEN USED IN CONJUNCTION ITS CONTROLLER         **
C    ** 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 NODE

      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

      PARAMETER ( NMAX = 60000 )
      PARAMETER ( MMX=50, MMY=50, MMZ=50, NCELL=MMX*MMY*MMZ )
      PARAMETER ( MAPSIZ = 13*NCELL )
      PARAMETER(NW1 = 11)

      REAL        RX(NMAX), RY(NMAX), RZ(NMAX)
      REAL        VX(NMAX), VY(NMAX), VZ(NMAX)
      REAL        FX(NMAX), FY(NMAX), FZ(NMAX)
      REAL        RCUT, DT, SIGMA
      REAL        VRCUT,DVRCUT,DVRC12
      REAL        SFX,SFY,SFZ,CELLIX,CELLIY,CELLIZ
      REAL        DENS, DENLJ, EQTEMP
      REAL        MASS, FREE, PI, TEMP
      REAL        E, V, K, W
      REAL        EN, VN, KN, PRES
      REAL        ACE, ACV, ACK, ACP
      REAL        ACESQ, ACVSQ, ACKSQ, ACPSQ
      REAL        AVE, AVV, AVK, AVP, AVT
      REAL        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.)

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.8
C  ENTER THE DENSITY
      DENS = 0.75
C  ENTER THE POTENTIAL CUTOFF DISTANCE
      RCUT = 2.0**(1.0/6.0)
C  ENTER THE TIMESTEP
      DT = 0.01
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/REAL(NATM))**(1.0/3.0)
      RCUT   = RCUT*SIGMA
      M      = INT(1./RCUT)
      MASS   = 1.0
      DENLJ  = DENS
      DENS   = DENS/(SIGMA**3)
      DT     = DT*SIGMA
      FREE   = REAL(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.0
      ACV    = 0.0
      ACK    = 0.0
      ACP    = 0.0

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

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

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**   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)
        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

        E    = K+V
        EN   = E/REAL(NATM)
        VN   = V/REAL(NATM)
        KN   = K/REAL(NATM)
        TEMP = 2.0*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/REAL(NATM)
        VN   = V/REAL(NATM)
        KN   = K/REAL(NATM)
        TEMP = 2.0*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 = REAL(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.0) FLE = SQRT(ACESQ)
      IF (ACKSQ.GT.0.0) FLK = SQRT(ACKSQ)
      IF (ACVSQ.GT.0.0) FLV = SQRT(ACVSQ)
      IF (ACPSQ.GT.0.0) FLP = SQRT(ACPSQ)

      AVT = AVK/1.5
      FLT = FLK/1.5

      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
      PARAMETER ( NMAX = 60000 )
      PARAMETER ( MMX=50, MMY=50, MMZ=50, NCELL=MMX*MMY*MMZ )
      PARAMETER ( MAPSIZ = 13*NCELL )

      REAL         DT, RCUT, SIGMA, EQTEMP
      REAL         SR6,VRCUT,DVRCUT,DVRC12
      REAL         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,INT(1.0/RCUT)+2)
      MY    = MIN(MMY,INT(1.0/RCUT)+2)
      MZ    = MIN(MMZ,INT(1.0/RCUT)+2)

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

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

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

      CELLIX = REAL(MX)
      CELLIY = REAL(MY)
      CELLIZ = REAL(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.0)
      DVRCUT = -1.0*SR6*(SR6-0.5)/RCUT
      DVRC12 = 12.0*DVRCUT

      RETURN
      END



      SUBROUTINE FCC(NC)

C*********************************************************************
C**   ROUTINE TO CREATE AN FCC CONFIGUATION 
C*********************************************************************

      IMPLICIT CHARACTER ( A - Z )

      INTEGER     NMAX

      PARAMETER ( NMAX = 60000 )

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

      REAL        RX(NMAX), RY(NMAX), RZ(NMAX)
      REAL        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) = REAL(LX)/REAL(NC)+1.0E-6-0.5
            RY(NATM) = REAL(LY)/REAL(NC)+1.0E-6-0.5
            RZ(NATM) = REAL(LZ)/REAL(NC)+1.0E-6-0.5

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

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

            NATM = NATM+1
            RX(NATM) = (REAL(LX)+0.5)/REAL(NC)-1.0E-6-0.5
            RY(NATM) = (REAL(LY)+0.5)/REAL(NC)-1.0E-6-0.5
            RZ(NATM) = REAL(LZ)/REAL(NC)+1.0E-6-0.5
            CALL SRAND(IS,VX(NATM))
            CALL SRAND(IS,VY(NATM))
            CALL SRAND(IS,VZ(NATM))

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

 60         CONTINUE
 50       CONTINUE
 40     CONTINUE

      RETURN
      END



      SUBROUTINE SRAND(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*********************************************************************

      REAL*8 SCALE

      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.0*(SCALE*FLOAT(IS)-0.5)*1.0e-6

      RETURN
      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

      PARAMETER ( NMAX = 60000 )
      PARAMETER ( MMX=50, MMY=50, MMZ=50, NCELL=MMX*MMY*MMZ )
      PARAMETER ( MAPSIZ = 13*NCELL )

      REAL        RX(NMAX), RY(NMAX), RZ(NMAX)
      REAL        VX(NMAX), VY(NMAX), VZ(NMAX)
      REAL        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+INT((RX(I)*SFX+0.5)*CELLIX)
     &           +INT((RY(I)*SFY+0.5)*CELLIY)*MX
     &           +INT((RZ(I)*SFZ+0.5)*CELLIZ)*MX*MY
        LIST(I) = HEAD(ICELL)
        HEAD(ICELL) = I
 10     CONTINUE

      NEDGE = NMAX+1

C**   CREATE WEST FACE OF BOX

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

          J = LIST(J)
          GOTO 21
 20       CONTINUE

C**   CREATE EAST FACE OF BOX

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

          J = LIST(J)
          GOTO 31
 30       CONTINUE

C**   CREATE SOUTH FACE OF BOX

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

          J = LIST(J)
          GOTO 101
 100      CONTINUE

C**   CREATE NORTH FACE OF BOX

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

          J = LIST(J)
          GOTO 131
 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.EQ.0) GOTO 220
        NEDGE = NEDGE-1
        RX(NEDGE) = RX(J)
        RY(NEDGE) = RY(J)
        RZ(NEDGE) = RZ(J)-1.0
        ICELL = 1+INT((RX(NEDGE)*SFX+0.5)*CELLIX)
     &           +INT((RY(NEDGE)*SFY+0.5)*CELLIY)*MX
        LIST(NEDGE) = HEAD(ICELL)
        HEAD(ICELL) = NEDGE

        J = LIST(J)
        GOTO 221
 220    CONTINUE


C**   CREATE UP FACE OF BOX

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

        J = LIST(J)
        GOTO 271
 270    CONTINUE

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

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

      RETURN
      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
      PARAMETER ( NMAX = 60000 )
      PARAMETER ( MMX=50, MMY=50, MMZ=50, NCELL=MMX*MMY*MMZ )
      PARAMETER ( MAPSIZ = 13*NCELL )

      REAL        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 40 IY = 1, MY
          DO 40 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)
40          CONTINUE

      RETURN
      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
      PARAMETER ( NMAX = 60000 )

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

      INTEGER     I, NATM

      DT2   = DT/2.0
      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)-AINT(RX(I)+RX(I))
         RY(I) = RY(I)-AINT(RY(I)+RY(I))
         RZ(I) = RZ(I)-AINT(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

      RETURN
      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
      PARAMETER ( NMAX = 60000 )

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

      INTEGER     I, NATM

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

      DT2   = DT/2.0
      K   = 0.0

      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.5*K

      RETURN
      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

      PARAMETER ( NMAX = 60000 )
      PARAMETER ( MMX=50, MMY=50, MMZ=50, NCELL=MMX*MMY*MMZ )
      PARAMETER ( MAPSIZ = 13*NCELL)

      REAL        RX(NMAX), RY(NMAX), RZ(NMAX)
      REAL        FX(NMAX), FY(NMAX), FZ(NMAX)
      REAL        RCUT, SIGMA
      REAL        VRCUT,DVRCUT,DVRC12
      REAL        V,K,W
      REAL        RXI, RYI, RZI, FXIJ, FYIJ, FZIJ, FIJ, RCUTSQ
      REAL        SIGSQ, FXI, FYI, FZI, SR2, SR6, VIJ, WIJ
      REAL        RIJ, RIJSQ, RXIJ, RYIJ, RZIJ
      REAL        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.0
        FY(I) = 0.0
        FZ(I) = 0.0
10      CONTINUE

      V = 0.0
      W = 0.0

C    ** LOOP OVER ALL CELLS **

      DO 5000 ICELL = 1 , 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 (EDGEI) GOTO 2200
          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   = SQRT(RIJSQ)
              SR2   = SIGSQ/RIJSQ
              SR6   = SR2*SR2*SR2
              VIJ   = SR6*(SR6-1.0)-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
              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
 2200     CONTINUE

C          ** LOOP OVER NEIGHBOURING CELLS **

        JCELL0 = 13*(ICELL-1)
        DO 4000 NABOR = 1, 13

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

          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   = SQRT(RIJSQ)
              SR2   = SIGSQ/RIJSQ
              SR6   = SR2*SR2*SR2
              VIJ   = (SR6*(SR6-1.0)-VRCUT-DVRC12*(RIJ-RCUT))*0.5
              WIJ   = SR6*(SR6-0.5)+DVRCUT*RIJ
              FIJ   = WIJ/RIJSQ
              WIJ   = 0.5*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.0 *V
      W = 48.0*W/3.0

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

      RETURN
      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
      PARAMETER ( NMAX = 60000 )

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

      INTEGER      NATM, I

      SCALEF = SQRT(EQTEMP/(2.0*K/(3.0*REAL(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

      RETURN
      END
