subroutine CalculateProperties()
  use simulation
  use scattering
  use particles
  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, k, bin, TempEneLevel,EneLevel,Temp_nAlBin
! 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.
  bin_ParNumber(:) = 0
  bin_ParAveEne(:) = 0.0d0
  bin_ParAveEneV(:) = 0.0d0
  bin_ParDrftVel(:) = 0.0d0
  bin_ParVal(:,:) = 0.0d0
  binhist_ParNumber(:,:) = 0
  num_Carriers = 0
  do i = 1,num_maxCarriers
    if (Par_use(i).eq.1) then ! used particle
     num_Carriers = num_Carriers + 1
     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) + Al_split_L_G(Par_Al(i))
     elseif (Par_valley(i).eq.3) then ! X Valley 
        TempEnergy = Par_ene(i) + Al_split_X_G(Par_Al(i))
     else 
     endif
     TempEneLevel = int(TempEnergy/StepE_eV)  ! for level frequency 
     if (TempEneLevel.ge.num_EnergyLevel) then
        TempEneLevel = num_EnergyLevel
     endif        
     EneLevel = TempEneLevel
     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
                EneLevel = j
            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 
                EneLevel = j
            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
                EneLevel = num_EnergyLevel
            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
                EneLevel = j
            endif
            prevBoundaryEeV = BoundaryEeV
        enddo
     endif
     TotParEne = TotParEne + TempEnergy
     TotEneValley(Par_valley(i)) = TotEneValley(Par_valley(i)) + Par_ene(i)
     factor = 1./(1.+nonParabolPara2(Par_Al(i),Par_valley(i))*Par_ene(i))
     do j = 1, DIM
        Par_vel(i,j) = hbar*Par_wvk(i,j)*factor/eff_mass(Par_Al(i),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
! Carrier density profile bin profile
     bin = num_Bins
     do j = 1,num_Bins ! 
         if ((Par_pos(i,1).ge.(loc_bin(j)-0.5*BinSize)).and.(Par_pos(i,1).lt.(loc_bin(j)+0.5*BinSize))) then
            bin = j
         endif           
     enddo
     bin_ParNumber(bin) = bin_ParNumber(bin) + 1 
     bin_ParAveEne(bin) = bin_ParAveEne(bin) + TempEnergy
     bin_ParAveEneV(bin) = bin_ParAveEneV(bin) + Par_ene(i)
     bin_ParDrftVel(bin) = bin_ParDrftVel(bin) +  Par_vel(i,1)
     bin_ParVal(Par_valley(i),bin) = bin_ParVal(Par_valley(i),bin) + 1.0
     binhist_ParNumber(bin,EneLevel) = binhist_ParNumber(bin,EneLevel) + 1 
    endif
  enddo  ! num_Carriers
  do j = 1,num_Bins ! 
     bin_ParNumDen(j) =  real(bin_ParNumber(j),kind = 8)/real(num_Carriers,kind = 8)  
     if (bin_ParNumber(j).gt.0) then
        bin_ParAveEne(j) = bin_ParAveEne(j)/real(bin_ParNumber(j),kind = 8)    
        bin_ParAveEneV(j) = bin_ParAveEneV(j)/real(bin_ParNumber(j),kind = 8)        
        bin_ParDrftVel(j) = bin_ParDrftVel(j)/real(bin_ParNumber(j),kind = 8)        
     else 
        bin_ParAveEne(j) = 0.0d0
        bin_ParAveEneV(j) = 0.0d0
        bin_ParDrftVel(j) = 0.0d0
     endif
     do i = 1, num_EnergyLevel
        if (bin_ParNumber(j).gt.0) then
            binhist_ParNumDen(j,i) = real(binhist_ParNumber(j,i),kind = 8)/real(bin_ParNumber(j),kind = 8)
        else
            binhist_ParNumDen(j,i) = 0.0d0
        endif
     enddo
     do i = 1, num_Valley
        if (bin_ParNumber(j).gt.0) then
            bin_ParVal(i,j) = bin_ParVal(i,j)/real(bin_ParNumber(j),kind = 8)
        else
            bin_ParVal(i,j) = 0.0d0
        endif
     enddo
  enddo
    ! average velocity in x, y, z direction 
  AveParEne =  TotParEne/real(num_Carriers,kind = 8)  
  do i = 1, num_Valley
    if (FreqValley(i).gt.0) then
        AveEneValley(i) = TotEneValley(i)/real(FreqValley(i),kind = 8)
    else
        AveEneValley(i) = 0.0d0
    endif
    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
        if (FreqValley(i).gt.0) then
            AveParValleyVel(i,j) = TempTotValleyVel(i,j)/real(FreqValley(i),kind = 8)
        else
            AveParValleyVel(i,j) = 0.0d0
        endif        
    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
        if (FreqValley(i).gt.0) then
            RelEneLevelValleyHist(i,j) = real(EneLevelValleyHist(i,j),kind = 8)/real(FreqValley(i),kind = 8)
        else
            RelEneLevelValleyHist(i,j) = 0.0d0
        endif
    enddo 
  enddo
! if properties for node center
  if (PropertyNode.eq.1) then
    bin_ParNumDen(1) =  2.0*bin_ParNumDen(1)  
    bin_ParNumber(1) =  2*bin_ParNumber(1) 
    bin_ParNumDen(num_Bins) =  2.0*bin_ParNumDen(num_Bins)
    bin_ParNumber(num_Bins) =  2*bin_ParNumber(num_Bins) 
  endif
return
end

subroutine PhononCountingAndEnergy(Time_step_output)
  use simulation
  use scattering
  use particles !   real ( kind = 8 ), dimension(num_Valley,DIM) :: TempTotValleyVel
  integer:: i, j, k !num_ScatMech,num_Valley
  real ( kind = 8 )::Time_step_output
  TA_bin_PhononEne(:) = 0.0d0
  TA_bin_PhononNum(:) = 0
  TotalP_Num = 0
  TotalP_Ene = 0.0d0
  TotalP_EneWm = 0.0d0
  do j = 1,num_Bins 
    TA_bin_PhononNum(j)= TA_bin_PhononNum(j)+BinScatCount(j,2,1)+BinScatCount(j,4,1)+BinScatCount(j,6,1)-(BinScatCount(j,3,1)+BinScatCount(j,5,1)+BinScatCount(j,7,1))
    TA_bin_PhononNum(j)= TA_bin_PhononNum(j)+BinScatCount(j,2,2)+BinScatCount(j,4,2)+BinScatCount(j,6,2)+BinScatCount(j,8,2)-(BinScatCount(j,3,2)+BinScatCount(j,5,2)+BinScatCount(j,7,2)+BinScatCount(j,9,2))
    TA_bin_PhononNum(j)= TA_bin_PhononNum(j)+BinScatCount(j,2,3)+BinScatCount(j,4,3)+BinScatCount(j,6,3)+BinScatCount(j,8,3)-(BinScatCount(j,3,3)+BinScatCount(j,5,3)+BinScatCount(j,7,3)+BinScatCount(j,9,3))  
    TA_bin_PhononEne(j)= TA_bin_PhononEne(j)+real((BinScatCount(j,2,1)-BinScatCount(j,3,1)),kind = 8)*polar_Ep_G+real((BinScatCount(j,4,1)-BinScatCount(j,5,1)),kind = 8)*IntVal_Ep_G_L+real((BinScatCount(j,6,1)-BinScatCount(j,7,1)),kind = 8)*IntVal_Ep_G_X
    TA_bin_PhononEne(j)= TA_bin_PhononEne(j)+real((BinScatCount(j,2,2)-BinScatCount(j,3,2)),kind = 8)*polar_Ep_L+real((BinScatCount(j,4,2)-BinScatCount(j,5,2)),kind = 8)*IntVal_Ep_L_G+real((BinScatCount(j,6,2)-BinScatCount(j,7,2)),kind = 8)*IntVal_Ep_L_L+real((BinScatCount(j,8,2)-BinScatCount(j,9,2)),kind = 8)*IntVal_Ep_L_X
    TA_bin_PhononEne(j)= TA_bin_PhononEne(j)+real((BinScatCount(j,2,3)-BinScatCount(j,3,3)),kind = 8)*polar_Ep_X+real((BinScatCount(j,4,3)-BinScatCount(j,5,3)),kind = 8)*IntVal_Ep_X_G+real((BinScatCount(j,6,3)-BinScatCount(j,7,3)),kind = 8)*IntVal_Ep_X_L+real((BinScatCount(j,8,3)-BinScatCount(j,9,3)),kind = 8)*IntVal_Ep_X_X
    TotalP_Num = TotalP_Num + TA_bin_PhononNum(j)
    TotalP_Ene = TotalP_Ene + TA_bin_PhononEne(j)          
    TA_bin_PhononEneWm(j) = (TA_bin_PhononEne(j)*e_c)*n_d/real(num_InitNode,kind = 8)/Time_step_output ! J/s- BinSize*TA_bin_ParNumber(j)
    if ((j.gt.num_boundBin).and.(j.le.(num_Bins-num_boundBin))) then
      TotalP_EneWm = TotalP_EneWm + TA_bin_PhononEneWm(j)
    endif
  enddo
  ! 0714 
  TotalScatCount = 0
  TotalopScatCount = 0
  do i = 1,num_Valley 
    do j= 1,num_ScatMech
      TotalScatCount = TotalScatCount + AllScatCount(j,i)
    enddo 
  enddo
  TotalopScatCount = TotalScatCount - AllScatCount(1,1)- AllScatCount(1,2) - AllScatCount(1,3)
  
  TotalP_EneWm = TotalP_EneWm/real((num_Bins-num_boundBin-num_boundBin),kind = 8) ! except 4 bins 
 ! 1: acoustic phonons, 2-3: ab -> em  polar optical phonons  4-5: G-L iv  6-7: G-X iv 1. acoustic phonons 2-3 polar optical phonons, 4-5: L-G 6-7: L-L 8-9: L-X intervalley
 ! 1- acoustic phonons, 2-3: polar optical phonons, 4-5: X-G, 6-7: X-L, 8-9: X-X 
 ! polar_Ep_G, polar_Ep_L, polar_Ep_X, IntVal_Ep_G_L, IntVal_Ep_G_X, IntVal_Ep_L_G, IntVal_Ep_L_L, IntVal_Ep_L_X, IntVal_Ep_X_G, IntVal_Ep_X_L, IntVal_Ep_X_X
  polar_G_Num = AllScatCount(2,1)-AllScatCount(3,1)
  polar_G_Ene = real((AllScatCount(2,1)-AllScatCount(3,1)),kind = 8)*polar_Ep_G
  IntVal_G_L_Num = AllScatCount(4,1)-AllScatCount(5,1)
  IntVal_G_L_Ene = real((AllScatCount(4,1)-AllScatCount(5,1)),kind = 8)*IntVal_Ep_G_L
  IntVal_G_X_Num = AllScatCount(6,1)-AllScatCount(7,1)
  IntVal_G_X_Ene = real((AllScatCount(6,1)-AllScatCount(7,1)),kind = 8)*IntVal_Ep_G_X
  polar_L_Num = AllScatCount(2,2)-AllScatCount(3,2)
  polar_L_Ene = real((AllScatCount(2,2)-AllScatCount(3,2)),kind = 8)*polar_Ep_L
  IntVal_L_G_Num = AllScatCount(4,2)-AllScatCount(5,2)
  IntVal_L_G_Ene = real((AllScatCount(4,2)-AllScatCount(5,2)),kind = 8)*IntVal_Ep_L_G
  IntVal_L_L_Num = AllScatCount(6,2)-AllScatCount(7,2)
  IntVal_L_L_Ene = real((AllScatCount(6,2)-AllScatCount(7,2)),kind = 8)*IntVal_Ep_L_L
  IntVal_L_X_Num = AllScatCount(8,2)-AllScatCount(9,2)
  IntVal_L_X_Ene = real((AllScatCount(8,2)-AllScatCount(9,2)),kind = 8)*IntVal_Ep_L_X
  polar_X_Num = AllScatCount(2,3)-AllScatCount(3,3)
  polar_X_Ene = real((AllScatCount(2,3)-AllScatCount(3,3)),kind = 8)*polar_Ep_X
  IntVal_X_G_Num = AllScatCount(4,3)-AllScatCount(5,3)
  IntVal_X_G_Ene = real((AllScatCount(4,3)-AllScatCount(5,3)),kind = 8)*IntVal_Ep_X_G
  IntVal_X_L_Num = AllScatCount(6,3)-AllScatCount(7,3)
  IntVal_X_L_Ene = real((AllScatCount(6,3)-AllScatCount(7,3)),kind = 8)*IntVal_Ep_X_L
  IntVal_X_X_Num = AllScatCount(8,3)-AllScatCount(9,3)
  IntVal_X_X_Ene = real((AllScatCount(8,3)-AllScatCount(9,3)),kind = 8)*IntVal_Ep_X_X
  return
end

! 0720 
subroutine PropertiesForSummary()
  use simulation
  use scattering
  use particles !   real ( kind = 8 ), dimension(num_Valley,DIM) :: TempTotValleyVel
  integer:: i   
  time_fave_period = time_fave_period + 1.0 
  Tfave_PhonEWm = 0.0
  Tfave_eNum = 0.0
  Tfave_vel = 0.0
  Tfave_TotEne =  0.0
  Tfave_PhonEWm_eb = 0.0
  Tfave_eNum_eb  = 0.0
  Tfave_vel_eb  = 0.0
  Tfave_TotEne_eb  = 0.0
  do i = 1, num_Bins   
    if (PhononSim.ne.0) then  !  Phonon simulation 11/2012  
      Tfave_bin_TpA(i) = Tfave_bin_TpA(i) + TABin_TemperAP(i)
      Tfave_bin_TpO(i) = Tfave_bin_TpO(i) + TABin_TemperOP(i)
      Tfave_bin_qsupA(i) = Tfave_bin_qsupA(i) + TABinQsupA(i)
      Tfave_bin_qsupO(i) = Tfave_bin_qsupO(i) + TABinQsupO(i)
      Tfave_bin_qsupT(i) = Tfave_bin_qsupT(i) + TABinQsupTotal(i)
      Tfave_bin_epWm2(i) = Tfave_bin_epWm2(i) + TAepInterEneWm2(i)
      Tfave_bin_ppWm2(i) = Tfave_bin_ppWm2(i) + TAppInterEneWm2(i)
      fave_bin_TpA(i) = Tfave_bin_TpA(i)/time_fave_period
      fave_bin_TpO(i) = Tfave_bin_TpO(i)/time_fave_period
      fave_bin_qsupA(i) = Tfave_bin_qsupA(i)/time_fave_period
      fave_bin_qsupO(i) = Tfave_bin_qsupO(i)/time_fave_period
      fave_bin_qsupT(i) = Tfave_bin_qsupT(i)/time_fave_period
      fave_bin_epWm2(i) = Tfave_bin_epWm2(i)/time_fave_period
      fave_bin_ppWm2(i) = Tfave_bin_ppWm2(i)/time_fave_period
    endif
    Tfave_bin_eNum(i) = Tfave_bin_eNum(i) + TA_bin_ParNumber(i)
    Tfave_bin_vel(i) = Tfave_bin_vel(i) + TA_bin_ParDrftVel(i)
    Tfave_bin_PhonEWm(i) = Tfave_bin_PhonEWm(i) + TA_bin_PhononEneWm(i)
    Tfave_bin_Eneev(i) = Tfave_bin_Eneev(i) + TA_bin_ParAveEne(i)  
    fave_bin_eNum(i) = Tfave_bin_eNum(i)/time_fave_period
    fave_bin_vel(i) = Tfave_bin_vel(i)/time_fave_period
    fave_bin_PhonEWm(i) = Tfave_bin_PhonEWm(i)/time_fave_period
    fave_bin_Eneev(i) = Tfave_bin_Eneev(i)/time_fave_period
    BoundFactor = 1.0
    if (PoissonSetting.eq.1) then      
      Tfave_bin_PoiPot(i) = Tfave_bin_PoiPot(i) + TA_Node_PoiPot(i)
      Tfave_bin_AllPot(i) = Tfave_bin_AllPot(i) + TA_Node_Potential(i)
      fave_bin_PoiPot(i) = Tfave_bin_PoiPot(i)/time_fave_period
      fave_bin_AllPot(i) = Tfave_bin_AllPot(i)/time_fave_period
    endif
    if (PropertyNode.eq.1) then 
      if ((i.eq.1).or.(i.eq.num_Bins)) then
        BoundFactor = 0.5
      endif
    endif
! if (PoissonSetting.eq.1) then      Tfave_bin_PoiPot(i) = Tfave_bin_PoiPot(i) + TA_Node_PoiPot(i)      Tfave_bin_AllPot(i) = Tfave_bin_AllPot(i) + TA_Node_Potential(i)      fave_bin_PoiPot(i) = Tfave_bin_PoiPot(i)/time_fave_period      fave_bin_AllPot(i) = Tfave_bin_AllPot(i)/time_fave_period      if ((i.eq.1).or.(i.eq.num_Bins)) then        BoundFactor = 0.5      endif endif
    ! Phonon simulation 11/2012    
    Tfave_PhonEWm = Tfave_PhonEWm + fave_bin_eNum(i)*fave_bin_PhonEWm(i)*BoundFactor
    Tfave_eNum = Tfave_eNum + fave_bin_eNum(i)*BoundFactor
    Tfave_vel = Tfave_vel + fave_bin_eNum(i)*fave_bin_vel(i)*BoundFactor
    Tfave_TotEne =  Tfave_TotEne + fave_bin_eNum(i)*fave_bin_Eneev(i)*BoundFactor
    if ((i.gt.num_boundBin).and.(i.lt.(num_Bins-num_boundBin))) then 
      Tfave_PhonEWm_eb = Tfave_PhonEWm_eb  + fave_bin_eNum(i)*fave_bin_PhonEWm(i)
      Tfave_eNum_eb  = Tfave_eNum_eb + fave_bin_eNum(i)
      Tfave_vel_eb  = Tfave_vel_eb  + fave_bin_eNum(i)*fave_bin_vel(i)
      Tfave_TotEne_eb  =  Tfave_TotEne_eb + fave_bin_eNum(i)*fave_bin_Eneev(i)
    endif
  enddo
  fave_PhonEWm = Tfave_PhonEWm/Tfave_eNum
  fave_eNum = Tfave_eNum
  fave_vel = Tfave_vel/Tfave_eNum
  fave_TotEne =  Tfave_TotEne/Tfave_eNum    
  fave_PhonEWm_eb = Tfave_PhonEWm_eb/Tfave_eNum_eb
  fave_eNum_eb  = Tfave_eNum_eb
  fave_vel_eb  = Tfave_vel_eb/Tfave_eNum_eb
  fave_TotEne_eb  =  Tfave_TotEne_eb/Tfave_eNum_eb
  fave_je = n_d*fave_eNum*fave_vel*e_c*0.0001/real(num_Carriers,kind = 8)
  fave_jeA = fave_je*box_len(1)
  if (PropertyNode.eq.1) then 
    fave_je_eb =  n_d*fave_eNum_eb*fave_vel_eb*e_c*0.0001/real((num_Bins-(num_boundBin+num_boundBin-1))*num_InitBin,kind = 8)
    fave_jeA_eb = fave_je_eb*(box_len(1)-BinSize*real((num_boundBin+num_boundBin-1),kind = 8))
  else
    fave_je_eb =  n_d*fave_eNum_eb*fave_vel_eb*e_c*0.0001/real((num_Bins-(num_boundBin+num_boundBin))*num_InitBin,kind = 8)
    fave_jeA_eb = fave_je_eb*(box_len(1)-BinSize*real((num_boundBin+num_boundBin),kind = 8))
  endif
  Tfave_InOut_B = Tfave_InOut_B + EneBound/(real(step_output2,kind = 8)*time_step*DerivedCellV)
  Tfave_InOut_LB = Tfave_InOut_LB + EneBoundL/(real(step_output2,kind = 8)*time_step)
  Tfave_InOut_RB = Tfave_InOut_RB + EneBoundR/(real(step_output2,kind = 8)*time_step)
  Tfave_G_ratio = Tfave_G_ratio + RelFreqValley(1)
  fave_InOut_B = Tfave_InOut_B/time_fave_period
  fave_InOut_LB = Tfave_InOut_LB/time_fave_period
  fave_InOut_RB = Tfave_InOut_RB/time_fave_period
  fave_G_ratio = Tfave_G_ratio/time_fave_period 
  return
end

! Initialize vel_InitBias, Initial_je