!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! To find initial state !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine PreRunBulk(nAlBin) use particles use scattering use simulation implicit none integer PRstep, PRcount_step, i, j,nAlBin real ( kind = 8 ):: PRtotalEnergyFreq, PRtotalValley,TempEInjLB,TempEInjRB real ( kind = 8 ):: PRwvecAbs,PRFactorVel,BiasWvecX,theta,upperLim,lowerLim,Portion,wvecX,velX real ( kind = 8 ):: FactorVel, TempPRVelp, TempPRVeln, TempTimePRVelp, TempTimePRVeln, TempPRVelx, TempTimePRVelx, LowerEne, UpperEne real ( kind = 8 ), dimension(3) :: TempPRval, TempTimePRval real ( kind = 8 ), dimension(:), allocatable :: TempPREne, TempTimePREne PRtime = 0. ! in simulation module PRstep = 0 PRcount_step = 0 allocate ( TempPREne(num_EnergyLevel) ) allocate ( TempTimePREne(num_EnergyLevel) ) TempTimePRVelp = 0.0 TempTimePRVeln = 0.0 TempTimePRVelx = 0.0 TempTimePRval(1:3) = 0.0 TempTimePREne(1:num_EnergyLevel) = 0.0 PRtotalEnergyFreq = 0.0 PRtotalValley = 0.0 do while(PRtime.le.PRtime_max) PRstep = PRstep + 1 PRtime = time_step*real(PRstep,kind = 8) call PreRunfreeflightscatter(nAlBin) ! properties calculation and fileoutput if (PRtime.ge.PRtime_early) then ! for early time ! Energy distribution TempPRVelx = 0.0 TempPRVelp = 0.0 TempPRVeln = 0.0 TempPRval(1:3) = 0.0 TempPREne(1:num_EnergyLevel) = 0.0 do i = 1, PRParNum FactorVel = 1./(1.+nonParabolPara2(nAlBin,PreRun_Par_valley(i))*PreRun_Par_ene(i)) do j = 1, DIM PreRun_Par_vel(i,j) = hbar*PreRun_Par_wvk(i,j)*FactorVel/eff_mass(nAlBin,PreRun_Par_valley(i)) enddo TempPRVelx = TempPRVelx + PreRun_Par_vel(i,1) if (PreRun_Par_vel(i,1).ge.0) then TempPRVelp = TempPRVelp + 1.0 else TempPRVeln = TempPRVeln + 1.0 endif TempPRval(PreRun_Par_valley(i)) = TempPRval(PreRun_Par_valley(i)) + 1.0 if (PreRun_Par_ene(i).lt.0) then PreRun_Par_ene(i) = 0.0 endif if (PreRun_Par_ene(i).ge.1.0) then PreRun_Par_ene(i) = 0.999999999 endif do j = 1, num_EnergyLevel LowerEne = StepE_eV*real((j-1),kind = 8) UpperEne = StepE_eV*real(j,kind = 8) if ((PreRun_Par_ene(i).ge.LowerEne).and.(PreRun_Par_ene(i).lt.UpperEne)) then TempPREne(j) = TempPREne(j) + 1.0 endif enddo enddo ! properties TempTimePRVelx = TempTimePRVelx + TempPRVelx/real(PRParNum,kind = 8) TempTimePRVelp = TempTimePRVelp + TempPRVelp/real(PRParNum,kind = 8) TempTimePRVeln = TempTimePRVeln + TempPRVeln/real(PRParNum,kind = 8) do j = 1, 3 TempTimePRval(j) = TempTimePRval(j) + TempPRval(j)/real(PRParNum,kind = 8) enddo do j = 1, num_EnergyLevel TempTimePREne(j) = TempTimePREne(j) + TempPREne(j)/real(PRParNum,kind = 8) enddo PRcount_step = PRcount_step + 1 endif enddo vel_InitBias = TempTimePRVelx/real(PRcount_step,kind = 8) TempTimePRVelp = TempTimePRVelp/real(PRcount_step,kind = 8) ! Positive Speed TempTimePRVeln = TempTimePRVeln/real(PRcount_step,kind = 8) ! negative Speed ! num_InjLB = (num_InjLB+num_InjRB)*TempTimePRVelp/(TempTimePRVelp+TempTimePRVeln) ! num_InjRB = (num_InjLB+num_InjRB)*TempTimePRVeln/(TempTimePRVelp+TempTimePRVeln) InjParPerStep = vel_InitBias*time_step*real(num_Carriers,kind = 8)/box_len(1) ! BinSize = box_len(1) RelCumFreqDist(:,:) = 0.0 do j = 1,num_EnergyLevel PRtotalEnergyFreq = PRtotalEnergyFreq + TempTimePREne(j) enddo PRtotalValley = TempTimePRval(1) + TempTimePRval(2) + TempTimePRval(3) PRValley(1) = TempTimePRval(1)/PRtotalValley do j = 2, 3 PRValley(j) = PRValley(j-1) + TempTimePRval(j)/PRtotalValley enddo RelCumFreqDist(nAlBin,1) = TempTimePREne(1)/PRtotalEnergyFreq do j = 2, num_EnergyLevel RelCumFreqDist(nAlBin,j) = RelCumFreqDist(nAlBin,j-1) + TempTimePREne(j)/PRtotalEnergyFreq enddo do j = 1,num_EnergyLevel TempTimePREne(j) = TempTimePREne(j)/PRtotalEnergyFreq enddo if (nAlBin.eq.(num_BinsForAl+1)) then ! Calculation for Injection num_InjLB = 0.0 endif num_InjRB = 0.0 do i = 1,num_EnergyLevel UpperEne = StepE_eV*real(i,kind = 8) PRwvecAbs = semh(nAlBin,1)*sqrt(UpperEne*(1.+ nonParabolPara(nAlBin,1)*UpperEne)) ! semh(valley)*sqrt(e*(1.+ nonParabolPara(valley)*e)) PRFactorVel = 1./(1.+2*nonParabolPara(nAlBin,1)*UpperEne) BiasWvecX = vel_InitBias*eff_mass(nAlBin,1)/hbar/PRFactorVel TempEInjLB = 0.0 TempEInjRB = 0.0 do j = 1,num_EnergyLevel theta = 0.5*real((2*j-1),kind = 8)*pi/real(num_EnergyLevel,kind = 8) upperLim = theta+0.5*pi/real(num_EnergyLevel,kind = 8) lowerLim = theta-0.5*pi/real(num_EnergyLevel,kind = 8) Portion = 0.5*(cos(lowerLim)-cos(upperLim)) wvecX = PRwvecAbs*cos(theta)+BiasWvecX velX = hbar*wvecX*PRFactorVel/eff_mass(nAlBin,1) if (velX.gt.0) then TempEInjLB = TempEInjLB+Portion*velX*time_step*(real(num_InitBin,kind = 8)/BinSize) endif if (velX.lt.0) then TempEInjRB = TempEInjRB-Portion*velX*time_step*(real(num_InitBin,kind = 8)/BinSize) endif enddo if (nAlBin.eq.(num_BinsForAl+1)) then num_InjLB = num_InjLB + TempTimePREne(i)*TempEInjLB endif num_InjRB = num_InjRB + TempTimePREne(i)*TempEInjRB enddo if (nAlBin.eq.num_BinsForAl) then open ( unit = 46, file = 'PreRunResultR.txt', status = 'replace', action = 'write' ) ! 0720 else open ( unit = 46, file = 'PreRunResultL.txt', status = 'replace', action = 'write' ) ! 0720 endif write(46,'(a,e16.8,a,e16.8,a,e16.8)'),'PRAveVelx = ',vel_InitBias,' num_InjLB =',num_InjLB, ' num_InjRB =',num_InjRB write(46,'(a,e16.8,a,e16.8,a,e16.8)'),'PRValley_G = ',PRValley(1),' PRValley_G+L =',PRValley(2), ' PRValley_G+L+X =',PRValley(3) write(46,'(a,a)'),'Ene',' Rel Cumulative ' do j = 1, num_EnergyLevel UpperEne = StepE_eV*real(j,kind = 8) write(46,'(f10.5,2e20.9)'),UpperEne,TempTimePREne(j),RelCumFreqDist(nAlBin,j) enddo close(46) ! 0720 return end subroutine PreRunfreeflightscatter(nAlBin) use particles use scattering use simulation integer valley,i,j real ( kind = 8 ), dimension(DIM) :: Temp_wvk ! carrier position and wave vector k integer :: Temp_valley, Temp_ScatMech,nAlBin real ( kind = 8 ) :: Temp_ene, Temp_tau, timeFreeFlight1,timeFreeFlight2,timeDrift,timeNewFlight, dt_remain, err1 do i = 1,num_ScatMech ! Reset counter for scattering frequency do j = 1,num_Valley ScatMechFreq(i,j) = 0. enddo enddo do i = 1, PRParNum ! num_Carriers ! loop for all carriers do j = 1,DIM !col_len = RCutoff_SiIn Temp_wvk(j) = PreRun_Par_wvk(i,j) enddo Temp_tau = PreRun_Par_tau(i) Temp_valley = PreRun_Par_valley(i) Temp_ene = PreRun_Par_ene(i) timeFreeFlight1 = Temp_tau ! Initial free-flight of the carriers if(timeFreeFlight1.ge.time_step)then timeDrift = time_step else timeDrift = timeFreeFlight1 endif call PreRunDrift(timeDrift,Temp_wvk,Temp_ene,Temp_valley,nAlBin) if(timeFreeFlight1.gt.time_step) goto 451 452 timeFreeFlight2=timeFreeFlight1 ! Free-flight and scatter part call scatter_carrier(nAlBin,Temp_wvk,Temp_ene,Temp_valley,Temp_ScatMech) 269 call random_number(harvest = RandomNum) if(RandomNum.le.1e-6) goto 269 timeNewFlight=-(log(RandomNum))*tau_max(nAlBin,Temp_valley) dt_remain = time_step - timeFreeFlight2 ! remaining time to scatter in dt-interval if(timeNewFlight.le.dt_remain)then timeDrift = timeNewFlight else timeDrift = dt_remain endif call PreRunDrift(timeDrift,Temp_wvk,Temp_ene,Temp_valley,nAlBin) timeFreeFlight2 = timeFreeFlight2 + timeNewFlight timeFreeFlight1 = timeFreeFlight2 if(timeFreeFlight1.lt.time_step) goto 452 451 timeFreeFlight1 = timeFreeFlight1 - time_step Temp_tau = timeFreeFlight1 do j = 1,DIM ! Map particle atributes PreRun_Par_wvk(i,j) = Temp_wvk(j) enddo PreRun_Par_tau(i) = Temp_tau PreRun_Par_valley(i) = Temp_valley PreRun_Par_ene(i) = Temp_ene enddo return end subroutine PreRunDrift(tau,Temp_wvk,Temp_ene,Temp_valley,nAlBin) use particles use scattering use simulation real ( kind = 8 ), dimension(DIM) :: Temp_wvk ! carrier position and wave vector k integer :: Temp_valley, i,nAlBin real ( kind = 8 ) :: Temp_ene, tau, sq_wvk_sum, abs_wvk, factor do i = 1, DIM Temp_wvk(i) = Temp_wvk(i) - echbar*tau*EFieldBg(i) enddo 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) factor = hhem(nAlBin,Temp_valley)*sq_wvk_sum Temp_ene = 2*factor/(1.+sqrt(1.+nonParabolPara4(nAlBin,Temp_valley)*factor)) return end subroutine PRInitialDist(nAlBin) use particles use scattering use simulation real ( kind = 8 ):: k,kx,ky,kz,e,tc,binFreq,fai,ct,st,err1 !! local integer nAlBin,valley,i,j,ebin call initial_energy() do i = 1, PRParNum 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(nAlBin,j-1)).and.(RandomNum.le.RelCumFreqDist(nAlBin,j))) then ebin = j endif enddo if (ebin.eq.1) then binFreq = RelCumFreqDist(nAlBin,ebin) e = StepE_eV*(RandomNum/RelCumFreqDist(nAlBin,ebin)) !(log(RelCumFreqDist(ebin))/log(RandomNum)) else binFreq = RelCumFreqDist(nAlBin,ebin) - RelCumFreqDist(nAlBin,ebin-1) e = StepE_eV*real((ebin-1),kind = 8) + StepE_eV *(1.-(log(RelCumFreqDist(nAlBin,ebin))-log(RandomNum))/(log(RelCumFreqDist(nAlBin,ebin))-log(RelCumFreqDist(nAlBin,ebin-1)))) endif call random_number(harvest = RandomNum) valley = 1 ! this is for GaAs materials k=semh(nAlBin,valley)*sqrt(e*(1.+ nonParabolPara(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 123 call random_number(harvest = RandomNum) ! Initial free-flight if(RandomNum.le.1.e-5) goto 123 tc = -(log(RandomNum))*tau_max(nAlBin,valley) ! random free flight time.. PreRun_Par_wvk(i,1) = kx ! Map particle atributes PreRun_Par_wvk(i,2) = ky PreRun_Par_wvk(i,3) = kz PreRun_Par_tau(i) = tc PreRun_Par_valley(i) = valley PreRun_Par_ene(i) = e ! if (Par_ene(i).ge.1) then err1 = 1.0 endif enddo return end