      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

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

      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 **
c--      OPEN (NW1,FILE='result')
c--      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
      WRITE(*,*) 'NC = ',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 'SYSTEM 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)

c--      CLOSE(NW1)

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

      END
