subroutine CalculateProperties()
  use simulation
  use scattering
  use particles
  use nodes
  implicit none
  real ( kind = 8 ) :: TempEnergy, BoundaryEeV, prevBoundaryEeV, factor
  real ( kind = 8 ), dimension(DIM) :: TempTotVel
  real ( kind = 8 ), dimension(num_Valley,DIM) :: TempTotValleyVel
  integer:: i, j, TempEneLevel
! Total Energy, average velocity in x, y, z direction , average velocity in x, y, z, Valley freq
  TotParEne = 0.
  TotEneValley(:) = 0.
  TempTotVel(:) = 0.
  TempTotValleyVel(:,:) = 0.
  FreqValley(:) = 0.
  EneLevelHist(:) = 0.
  EneLevelValleyHist(:,:) = 0.
  do i = 1,num_Carriers 	
     FreqValley(Par_valley(i)) = FreqValley(Par_valley(i)) + 1
     if (Par_valley(i).eq.1) then ! Gamma Valley 
        TempEnergy = Par_ene(i)
     elseif (Par_valley(i).eq.2) then ! L Valley 
        TempEnergy = Par_ene(i) + split_L_G
     elseif (Par_valley(i).eq.3) then ! X Valley 
        TempEnergy = Par_ene(i) + split_X_G
     else 
     endif
     TempEneLevel = int(TempEnergy/StepE_eV)  ! for level frequency 
     if (TempEneLevel.lt.3) then  
         prevBoundaryEeV = 0.
         do j = 1,TempEneLevel+2
            BoundaryEeV = StepE_eV*real(j,kind = 8)
            if ((TempEnergy.ge.prevBoundaryEeV).and.(TempEnergy.lt.BoundaryEeV)) then  
                EneLevelHist(j) = EneLevelHist(j)+1
                EneLevelValleyHist(Par_valley(i),j) = EneLevelValleyHist(Par_valley(i),j) +1
            endif 
            prevBoundaryEeV = BoundaryEeV
         enddo
     elseif (TempEneLevel.gt.(num_EnergyLevel-2))then
        prevBoundaryEeV = StepE_eV*real((TempEneLevel-3),kind = 8)
        do j = TempEneLevel-2,num_EnergyLevel
            BoundaryEeV = StepE_eV*real(j,kind = 8)
            if ((TempEnergy.ge.prevBoundaryEeV).and.(TempEnergy.lt.BoundaryEeV)) then  
                EneLevelHist(j) = EneLevelHist(j)+1
                EneLevelValleyHist(Par_valley(i),j) = EneLevelValleyHist(Par_valley(i),j) +1
            endif
            if (TempEnergy.ge.MaxE_eV) then  
                EneLevelHist(num_EnergyLevel) = EneLevelHist(num_EnergyLevel)+1
                EneLevelValleyHist(Par_valley(i),num_EnergyLevel) = EneLevelValleyHist(Par_valley(i),num_EnergyLevel) +1
            endif
            prevBoundaryEeV = BoundaryEeV
        enddo
     else 
        prevBoundaryEeV = StepE_eV*real((TempEneLevel-3),kind = 8)	
        do j = TempEneLevel-2,TempEneLevel+2
            BoundaryEeV = StepE_eV*real(j,kind = 8)
            if ((TempEnergy.ge.prevBoundaryEeV).and.(TempEnergy.lt.BoundaryEeV)) then  
                EneLevelHist(j) = EneLevelHist(j)+1
                EneLevelValleyHist(Par_valley(i),j) = EneLevelValleyHist(Par_valley(i),j) +1
            endif
            prevBoundaryEeV = BoundaryEeV
        enddo
     endif
     TotParEne = TotParEne + TempEnergy
     TotEneValley(Par_valley(i)) = TotEneValley(Par_valley(i)) + Par_ene(i)
     factor = 1./(1.+nonParabolPara2(Par_valley(i))*Par_ene(i))
     do j = 1, DIM
        Par_vel(i,j) = hbar*Par_wvk(i,j)*factor/eff_mass(Par_valley(i))
        TempTotVel(j) = TempTotVel(j) + Par_vel(i,j)
        TempTotValleyVel(Par_valley(i),j) = TempTotValleyVel(Par_valley(i),j) + Par_vel(i,j)
     enddo
!    Par_tau 
  enddo  
  ! average velocity in x, y, z direction 
  AveParEne =  TotParEne/real(num_Carriers,kind = 8)  
  do i = 1, num_Valley
    AveEneValley(i) = TotEneValley(i)/real(FreqValley(i),kind = 8)
    RelFreqValley(i) =  real(FreqValley(i),kind = 8)/real(num_Carriers,kind = 8) 
    AveParVel(i) = TempTotVel(i)/real(num_Carriers,kind = 8) 
    do j = 1, DIM
        AveParValleyVel(i,j) = TempTotValleyVel(i,j)/real(FreqValley(i),kind = 8)
    enddo
  enddo 
! Histogram
  do j = 1,num_EnergyLevel 
    RelEneLevelHist(j) = real(EneLevelHist(j),kind = 8)/real(num_Carriers,kind = 8) 
    do i = 1, num_Valley
        RelEneLevelValleyHist(i,j) = real(EneLevelValleyHist(i,j),kind = 8)/real(FreqValley(i),kind = 8)
    enddo 
  enddo
! Carrier density profile Node 
return
end