!     This subroutine distributes velocities to NAN particles
!     in order to make even distribution of kinetic energies
!     from 0 to 6kT. Half of this energy will be transferred
!     to the potential energy. (I Assume that the initial
!     configuration is quenched [POT.EN.=0])
!     More about initialization of velocities in Frenkel & Smit p 56-57
!     MSE 6270, Leonid Zhigilei

      SUBROUTINE Vel()
      INCLUDE 'common.h'
      DIMENSION V(NAN)

      xyz_loop: DO IXYZ=1,3
        SK=0.0d0
        QIN1=0.0d0
        DO I=1,NAN
          call random_number(rww)
          V(I)=2.*DELTA*SQRT(QTEM*BK*RWW/ENUNIT/XMASS(KTYPE(I)))
          call random_number(rww1)
          AWW=RWW1-0.5
          V(I)=SIGN(V(I),AWW)
          SK=SK+RWW
        END DO

!       FINITE SIZE CORRECTION
        P=NAN/2.0d0/SK
        VM=0.0d0
        DO I=1,NAN
          V(I)=V(I)*P
          VM=VM+V(I)*XMASS(KTYPE(I))
        END DO

!       TOTAL MOMENTUM CORRECTION
        VM=VM/NAN
        DO I=1,NAN
          V(I)=V(I)-VM/XMASS(KTYPE(I))
          QIN1=QIN1+V(I)*V(I)*XMASS(KTYPE(I))/2.0d0
        END DO

!       TEMPERATURE CORRECTION
        T1=QIN1*ENUNIT/BK/NAN/DELTA/DELTA
        WRITE(6,*) 'VEL------> T=',QTEM,', Tvel=',T1,' KELVIN.'
        PR=SQRT(QTEM/T1)
        QIN1=0.0d0
        DO I=1,NAN
          V(I)=V(I)*PR
          QIN1=QIN1+V(I)*V(I)*XMASS(KTYPE(I))/2.0d0
        END DO
        T2=QIN1*ENUNIT/BK/NAN/DELTA/DELTA
        WRITE(6,*) 'VEL------> T=',QTEM,',Corrected Tvel=',T2,' KELVIN.'

!       Now let's distribute the velocities to Q1D(1:3, 1:NAN)
        IF(IXYZ.LE.2) Then
          DO I=1,NAN
            Q1D(IXYZ,I)=V(I)
          ENDDO
        ELSEIF(NDIM.EQ.3) Then
          DO I=1,NAN
            Q1D(IXYZ,I)=V(I)
          ENDDO
        ENDIF

      END DO xyz_loop

      RETURN
      END