!     THIS SUBROUTINE SETS THE INITIAL PARAMETERS
!     MSE 6270, Leonid Zhigilei

      SUBROUTINE SetInit(DENSM,DENSN)
      INCLUDE 'common.h'

!     Counting particles of different types
      KTT(1:NTYPE)=0
      DO I=1,NAN
        DO J=1,NTYPE
          IF(KTYPE(I).EQ.J) KTT(J)=KTT(J)+1
        ENDDO
      ENDDO

      KTOT=0
      DO J=1,NTYPE
          KTOT=KTOT+KTT(J)
      ENDDO
      IF(KTOT.NE.NAN) Then
	  Print *,'KTYPE(I), NTYPE, and NAN do not match, Program stops'
        STOP
      ENDIF

      IF (NDIM.EQ.3) Then
        VOLUME=(XL*1.0D-08)*(YL*1.0D-08)*(ZL*1.0D-08) ! Volume in cm3
        TMASS=0.0d0
        DO J=1,NTYPE
          TMASS=TMASS+KTT(J)*XMASS(J)       ! Total mass in amu
        ENDDO
        DENSM=TMASS*AMUTOKG*1.0D+03/VOLUME  ! gr/cm3
        DENSN=NAN/VOLUME                    ! molecules/cm3
      Else
        DENSM=-1.
        DENSN=-1.
      EndIf

!     Let's define a few variables
      NO1=NAN-1
      NAN3=NAN+NAN+NAN
      XLHALF=XL*0.5d0
      YLHALF=YL*0.5d0
      ZLHALF=ZL*0.5d0

!     We will use DXR(I) to pick a value of force or energy that
!     corresponds to a given interparticle distance from the table
      DO I=1,NPOTS
       DXR(I)=NT/RM(I)
       RLIST(I) = RM(I)+RSKIN
       RLIST2(I) = RLIST(I)*RLIST(I)
      ENDDO

      FullList = .FALSE.                ! for pair potential
      IF(KEYBS.EQ.1) FullList = .TRUE.  ! for MC or S-W potential

!     Let's assume that we have up to 3 types of atoms
!     interacting via up to 6 pair potentials
      IJINDEX(1,1)=1
      IJINDEX(1,2)=3
      IJINDEX(1,3)=6
      IJINDEX(2,1)=3
      IJINDEX(2,2)=2
      IJINDEX(2,3)=5
      IJINDEX(3,1)=6
      IJINDEX(3,2)=5
      IJINDEX(3,3)=4

!     Defining the boundary layer
      NRIGID=0
      DO I=1,NAN
        VW(I)=0.0d0
        KRIGID(I)=0
        IF(KHIST(I).EQ.3) KRIGID(I)=KBOUND
        IF(KHIST(I).EQ.3) NRIGID=NRIGID+1
      ENDDO

!     It is convenient to multiply velocities by timestep DELTA
!     when a predictor-corrector algorithm is used.
      DO I=1,NAN3
        Q1(I)=Q1(I)*DELTA
        Q2(I)=0.0d0
        Q3(I)=0.0d0
        Q4(I)=0.0d0
        Q5(I)=0.0d0
      ENDDO
      IF (NDIM.EQ.2) Q1D(3,1:NAN)=0.0d0

!     G1 is used in Nordsieck Integration method
      DO J=1,NTYPE
        G1(J)=0.5d0*DELTA*DELTA/XMASS(J)
      ENDDO

!     Nordsieck predictor-corrector coefficients
      C1=3.d0/20.d0
      C2=251.d0/360.d0
      C3=11.d0/18.d0
      C4=1.d0/6.d0
      C5=1.d0/60.d0

      RETURN
      END