!****** B O L T Z M A N N  D I S T R I B !*********************************************************
!
!  F. Nattino, 04/2012
!
!  Select vibrational state according to Boltzmann distribution
!
!  Units: T = K, emax = eV  
!
!  USE: Constants (wn2J, eV2wn and kb) and Random ( function RAN0)
!
!  NOTE: we treat the doubly degenerate modes as two single degenerate modes
!        
  MODULE  BoltzmannDistribution
   USE Constants
   USE Random

! Frequencies from O.N. Ulenikov, E.S. Bekhtereva, S. Albert, S. Bauerecker, H. Hollenstein and M. Quack, mol. phys. 108, 1209 (2010)

  IMPLICIT NONE

  INTEGER, PARAMETER                    :: Nfreq = 9  ! Changing of freq implies changing also the nested do in FractionalPopulations!!
  REAL, DIMENSION(nfreq), PARAMETER     :: omega = (/ 2992.786D0, & 
    2142.583D0, 1004.548D0, 2250.828D0, 2250.828D0, 1292.500D0, 1292.500D0, 1035.920D0, 1035.920D0 /)
  INTEGER                               :: Ndeger
  INTEGER, DIMENSION(nfreq), PARAMETER  ::  g0   = (/ (1, Ndeger =1,Nfreq)   /) ! g0( i ) contains degeneracy of i-th mode 
  REAL,    ALLOCATABLE, SAVE            :: FractionalPopulation(:), VibrationalEnergy(:), BoltzmannFactor(:)
  INTEGER, ALLOCATABLE, SAVE            :: VibrationalStates(:,:)
  INTEGER, SAVE                         :: Numstates
  CHARACTER(*), PARAMETER               :: BoltzmannDistributionFile = 'BoltzmannDistribution.dat'

  CONTAINS

! -------------------------------------------------------------------------------------------------

  SUBROUTINE ComputeBoltzmannFactors( T , emax ) 

  IMPLICIT NONE
  REAL, INTENT(IN)          :: T , emax
  REAL,    ALLOCATABLE      :: FvT(:),Estate(:)
  INTEGER, ALLOCATABLE      :: state(:,:), order(:)
  LOGICAL, ALLOCATABLE      :: Mask(:)
  INTEGER, DIMENSION(nfreq) :: v,nvmax
  INTEGER                   :: i,j,k,l,m,n,o,p,q
  INTEGER                   :: NC, Nmode
  REAL                      :: glevel, g, Elevel, Den

! Find vmax per each vibrational eigenstate       
  DO i = 1, nfreq
    nvmax(i) = FLOOR( emax*eV2wn / omega(i) )
  ENDDO

  numstates=0

! Find how many levels we do have to consider
  DO i=0,nvmax(1)
   DO j=0,nvmax(2)
    DO k=0,nvmax(3)  
     DO l=0,nvmax(4)
      DO m=0,nvmax(5)
       DO n=0,nvmax(6)
        DO o=0,nvmax(7)
         DO p=0,nvmax(8)
          DO q=0,nvmax(9)
           v = (/ i,j,k,l,m,n,o,p,q /)
           Elevel = dot_product( dble(v) , omega)
           IF ( Elevel<= Emax * eV2wn) THEN
            Numstates = Numstates + 1
           ENDIF
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
     ENDDO
    ENDDO
   ENDDO
  ENDDO

! Allocate matrixes for storing data of the vibrational states
  ALLOCATE(FvT(numstates),order(numstates),state(numstates,nfreq),Estate(numstates),Mask(numstates))
  ALLOCATE(VibrationalEnergy(numstates),FractionalPopulation(numstates),VibrationalStates(numstates,nfreq), &
           BoltzmannFactor(numstates))

  NC = 0

! Compute boltzmann factors per each level
  DO i=0,nvmax(1)
   DO j=0,nvmax(2)
    DO k=0,nvmax(3)
     DO l=0,nvmax(4)
      DO m=0,nvmax(5)
       DO n=0,nvmax(6)
        DO o=0,nvmax(7)
         DO p=0,nvmax(8)
          DO q=0,nvmax(9)
           v = (/ i,j,k,l,m,n,o,p,q /)
           Elevel = dot_product( real(v) , omega)
           IF ( Elevel <= Emax* eV2wn) THEN
            NC = NC + 1
            State ( NC , : ) = v
            EState( NC )  = Elevel 
            glevel = 1.
            DO NMode=1,nfreq
             IF ( State( NC, NMode) >= 1 ) THEN
              IF(g0(nmode)==1) g=dble(1)
              IF(g0(nmode)==2) g=dble(state(nc,nmode))+1.D0
              IF(g0(nmode)==3) g=(dble(state(nc,nmode))+1.D0)*(dble(state(nc,nmode))+2.D0)/2.D0
              glevel=glevel*g
             ENDIF
            ENDDO
            FvT(nc)=glevel*exp(- Estate( NC ) * wn2J / ( kb * T ))
           ENDIF
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
     ENDDO
    ENDDO
   ENDDO
  ENDDO

! Order levels 
  Mask = .TRUE.

  DO i = 1, numstates  
    order(i) = minloc( Estate(:), 1, Mask(:))
    Mask( order(i) ) = .FALSE.
  ENDDO

  DO i = 1, numstates
    VibrationalEnergy(i)   = Estate( order(i) )
    BoltzmannFactor(i)     = FvT(    order(i) )
    VibrationalStates(i,:) = state(  order(i) , : ) 
  ENDDO


! Compute fractional population and write them in file
  OPEN(1,FILE=BoltzmannDistributionFile)
  WRITE(1,*)'   v1   v2   v3 v4_1 v4_2 v5_1 v5_2 v6_1 v6_2  Frac.Pop. En. rel(eV)'
  Den = SUM ( BoltzmannFactor )
  DO i = 1, numstates
    FractionalPopulation( i ) = BoltzmannFactor( i ) / Den
    WRITE(1,'(9(4x,i1),2(4x,f12.10))') ( VibrationalStates(i,j), j=1,9),FractionalPopulation(i),VibrationalEnergy(i)/eV2wn
  ENDDO
  CLOSE(1)

  PRINT '(1/,27("*"),A,27("*")1/)','  B O L T Z M A N N  D I S T R I B U T I O N  '
  PRINT '(A,f8.3)',' Frac. pop. of the first 10 vib. levels at T (K) =', T
  PRINT *,'   v1   v2   v3 v4_1 v4_2 v5_1 v5_2 v6_1 v6_2  Frac. Pop'
  DO i = 1,10
   PRINT '(9(4x,i1),4x,f7.5)', ( VibrationalStates(i,j), j=1,9),FractionalPopulation(i)
  ENDDO
  PRINT '(1/,A,f7.4,2A,1/)',' All other states up to ',emax,' eV (relative Vib. Energy) printed in ',BoltzmannDistributionFile
  
  IF (FractionalPopulation( Numstates ) >= 1.E-3 ) PRINT *,'WARNING! last state with population >= 1.E-3!'
  
  RETURN
  END SUBROUTINE  

! -------------------------------------------------------------------------------------------------

  SUBROUTINE SelectVibrationalState( vibstate )
  
  IMPLICIT NONE

  INTEGER, DIMENSION(nfreq), INTENT(OUT) :: vibstate
  INTEGER                                :: i
  REAL                                   :: low, high, Dice

  low=0.
  high=0.

  Dice = RAN0()
  DO i = 1, Numstates
    low  = high
    high = low + FractionalPopulation( i )
    IF ( dice >= low .and. dice<high ) THEN
      vibstate = VibrationalStates(i,:)
      GOTO 2001
    ENDIF
  ENDDO

2001   CONTINUE

  RETURN
  END SUBROUTINE

  END MODULE