!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! INITIALIZE CARRIER ENERGY AND WAVEVECTOR ACCORDING TO THE MAXWELL-BOLTZMANN STATISTICS ! Assumption: All carriers are initially in valley 1 for GaAs !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine initial_distribution() use particles use scattering use simulation real ( kind = 8 ):: k,kx,ky,kz,e, sumEeV, aveEeV, rr, tc, binFreq, fai, ct, st, err1 !! local real ( kind = 8 ):: tempPosX,LowAlBound,UpAlBound real ( kind = 8 ), dimension(DIM) :: Temp_wvk real ( kind = 8 ), dimension(DIM) :: Temp_vel ! carrier position and wave vector k real ( kind = 8 ):: sq_wvk_sum, abs_wvk, FactorEne, FactorVel integer valley, i ,j, ebin, bin, min, minbin, binTemp, Temp_nAlBin sumEeV = 0. bin_ParNumDen(:) = 0.0d0 bin_ParNumber(:) = 0 if (Use_prerun.eq.0) then call initial_energy() endif do i = 1, num_Carriers call random_number(harvest = RandomNum) ! init position setting - each bininitially tempPosX = box_len(1)*RandomNum + box_coor(1,1) ! or equal bin = 1 do j = 1,num_Bins ! if ((tempPosX.ge.(loc_bin(j)-0.5*BinSize)).and.(tempPosX.lt.(loc_bin(j)+0.5*BinSize))) then bin = j endif enddo if (InitialPosDist.eq.1) then Par_pos(i,1) = tempPosX bin_ParNumber(bin) = bin_ParNumber(bin) + 1 else if (bin_ParNumber(bin).ge.num_InitBin) then min = num_InitBin do j = 1,num_Bins ! binTemp = bin_ParNumber(j) if (PropertyNode.eq.1) then if ((j.eq.1).or.(j.eq.num_Bins)) then binTemp = binTemp*2 endif endif if (binTemp.lt.min) then min = binTemp minbin = j endif enddo if (PropertyNode.eq.1) then if (minbin.eq.1)then Par_pos(i,1) = loc_bin(minbin) + 0.5*BinSize*RandomNum elseif (minbin.eq.num_Bins) then Par_pos(i,1) = loc_bin(minbin)-0.5*BinSize + 0.5*BinSize*RandomNum else Par_pos(i,1) = loc_bin(minbin)-0.5*BinSize + BinSize*RandomNum endif else Par_pos(i,1) = loc_bin(minbin)-0.5*BinSize + BinSize*RandomNum endif bin_ParNumber(minbin) = bin_ParNumber(minbin) + 1 else Par_pos(i,1) = tempPosX bin_ParNumber(bin) = bin_ParNumber(bin) + 1 endif endif call random_number(harvest = RandomNum) Par_pos(i,2) = box_len(2)*RandomNum + box_coor(1,2) call random_number(harvest = RandomNum) Par_pos(i,3) = box_len(3)*RandomNum + box_coor(1,3) call DecideAlBin(Par_pos(i,1),Temp_nAlBin,LowAlBound,UpAlBound) call random_number(harvest = RandomNum) do while(RandomNum.eq.0) call random_number(harvest = RandomNum) enddo ebin = 1 do j = 2, num_EnergyLevel if ((RandomNum.gt.RelCumFreqDist(Temp_nAlBin,j-1)).and.(RandomNum.le.RelCumFreqDist(Temp_nAlBin,j))) then ebin = j endif enddo if (ebin.eq.1) then binFreq = RelCumFreqDist(Temp_nAlBin,ebin) e = StepE_eV*(RandomNum/RelCumFreqDist(Temp_nAlBin,ebin)) !(log(RelCumFreqDist(ebin))/log(RandomNum)) else binFreq = RelCumFreqDist(Temp_nAlBin,ebin) - RelCumFreqDist(Temp_nAlBin,ebin-1) e = StepE_eV*real((ebin-1),kind = 8) + StepE_eV *(1.-(log(RelCumFreqDist(Temp_nAlBin,ebin))-log(RandomNum))/(log(RelCumFreqDist(Temp_nAlBin,ebin))-log(RelCumFreqDist(Temp_nAlBin,ebin-1)))) endif call random_number(harvest = RandomNum) if (Use_prerun.ne.0) then valley = 1 if (RandomNum.le.PRValley(1)) then valley = 1 elseif (RandomNum.le.PRValley(2)) then valley = 2 elseif (RandomNum.le.PRValley(3)) then valley = 3 endif else valley = 1 ! this is for GaAs materials endif ! Initial wavevector k=semh(Temp_nAlBin,valley)*sqrt(e*(1.+ nonParabolPara(Temp_nAlBin,valley)*e)) call random_number(harvest = RandomNum) fai=2.*pi*RandomNum call random_number(harvest = RandomNum) ct=1.-2.*RandomNum st=sqrt(1.-ct*ct) kx=k*st*cos(fai) ky=k*st*sin(fai) kz=k*ct ! Initial free-flight 103 call random_number(harvest = RandomNum) if(RandomNum.le.1.e-5)go to 103 tc = -(log(RandomNum))*tau_max(Temp_nAlBin,valley) ! random free flight time.. ! vel_InitBias velocity bias... for initial current... 1. v 2.add 3. new k vector, new e. .. FactorEne, FactorVel Temp_wvk(1) = kx Temp_wvk(2) = ky Temp_wvk(3) = kz FactorVel = 1./(1.+nonParabolPara2(Temp_nAlBin,valley)*e) do j = 1, DIM Temp_vel(j) = hbar*Temp_wvk(j)*FactorVel/eff_mass(Temp_nAlBin,valley) enddo Temp_vel(1) = Temp_vel(1)+ vel_InitBias ! new wvk Temp_wvk(1) = eff_mass(Temp_nAlBin,valley)*Temp_vel(1)/hbar/FactorVel sq_wvk_sum = Temp_wvk(1)*Temp_wvk(1) + Temp_wvk(2)*Temp_wvk(2) + Temp_wvk(3)*Temp_wvk(3) abs_wvk = sqrt(sq_wvk_sum) FactorEne = hhem(Temp_nAlBin,valley)*sq_wvk_sum ! new energy e = 2*FactorEne/(1.+sqrt(1.+nonParabolPara4(Temp_nAlBin,valley)*FactorEne)) ! Map particle atributes Par_Al(i) = Temp_nAlBin Par_wvk(i,1) = Temp_wvk(1) Par_wvk(i,2) = Temp_wvk(2) Par_wvk(i,3) = Temp_wvk(3) Par_tau(i) = tc Par_valley(i) = valley Par_ene(i) = e Par_use(i) = 1 enddo ! if ((Par_ene(i).le.MaxE_eV).and.(Par_ene(i).gt.0)) then else err1 = 1 ! for debugging endif return end subroutine initial_energy() use scattering use simulation integer i,valley,nAlBin real ( kind = 8 ):: NonparaCoeffG, D_eeV, fe, Temp_EeV real ( kind = 8 ), dimension(:), allocatable :: totalFreq allocate ( totalFreq(num_BinsForAl+1) ) totalFreq(:)= 0. valley = 1 ! initially located in Gamma CumFreqDist(:,:) = 0. do nAlBin = 1, num_BinsForAl+1 do i = 1, num_EnergyLevel Temp_EeV = StepE_eV*real(i,kind = 8) NonparaCoeffG = (1+nonParabolPara2(nAlBin,valley)*Temp_EeV)*(sqrt(1+nonParabolPara(nAlBin,valley)*Temp_EeV)) D_eeV = NonparaCoeffG*e_c*Al_PrimV(nAlBin)*sqrt(2*Temp_EeV*e_c)*eff_mass(nAlBin,valley)*sqrt(eff_mass(nAlBin,valley))/(pi*pi*hbar*hbar*hbar) exp_e = exp(((Temp_EeV-E_F)*e_c)/(kB*Te)) fe = 1/(exp_e + 1) FreqDist(nAlBin,i) = fe*D_eeV*StepE_eV if (i.eq.1) then CumFreqDist(nAlBin,i) = FreqDist(nAlBin,i) else CumFreqDist(nAlBin,i) = FreqDist(nAlBin,i)+CumFreqDist(nAlBin,i-1) endif totalFreq(nAlBin) = totalFreq(nAlBin) + FreqDist(nAlBin,i) enddo do i = 1, num_EnergyLevel RelCumFreqDist(nAlBin,i) = CumFreqDist(nAlBin,i)/totalFreq(nAlBin) enddo enddo return end