!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Select Scattering Mechanism & Perform the Scattering Part Modifying Particle Atributes (kx, ky, kz, iv, energy) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine scatter_carrier(Temp_wvk,Temp_ene,Temp_valley_i) use simulation use scattering use particles real ( kind = 8 ), dimension(DIM) :: Temp_wvk ! carrier position and wave vector k integer :: Temp_valley_i, Temp_valley_f real ( kind = 8 ) :: Temp_ene !, Temp_tau integer loc, i_fix, i_top, i, select_mech real ( kind = 8 ):: rr, bound_lower, bound_upper, ee, err1 ! Calculate index to the scattering table ee = Temp_ene loc = int(ee/StepE_eV) ! if (Temp_ene.ge.1) then ! err1 = 1.0 ! endif if(loc.eq.0)loc=1 ! if zero energy .. it has minimum energy if(loc.gt.num_EnergyLevel) loc = num_EnergyLevel ! maximum ! Select scattering mechanism i_top = max_scatt_mech(Temp_valley_i) ! i_count ... total # of scattering algorithms.. call random_number(harvest = RandomNum) rr = RandomNum if(rr.ge.ren_scatt_table(loc,i_top,Temp_valley_i))then ! if larger than total scattering rate = sum over all scattering algorithms ScatMechFreq(i_top+1,Temp_valley_i)=ScatMechFreq(i_top+1,Temp_valley_i)+1 ! counting for self scattering goto 222 ! self-scattering endif if(rr.lt.ren_scatt_table(loc,1,Temp_valley_i))then ! the first mechanism... i_fix = 1 ScatMechFreq(i_fix,Temp_valley_i) = ScatMechFreq(i_fix,Temp_valley_i)+1 ! counting goto 111 endif if(i_top.gt.1)then do i=1,i_top-1 bound_lower = ren_scatt_table(loc,i,Temp_valley_i) bound_upper = ren_scatt_table(loc,i+1,Temp_valley_i) if(rr.ge.bound_lower.and.rr.lt.bound_upper)then i_fix = i + 1 ! index ScatMechFreq(i_fix,Temp_valley_i)=ScatMechFreq(i_fix,Temp_valley_i)+1 ! count goto 111 endif enddo endif 111 continue ! Perform scattering (change energy and randomize momentum) select_mech = flag_mech(i_fix,Temp_valley_i) Temp_valley_f = i_valley(i_fix,Temp_valley_i) ! valley destination.. if (Temp_ene.ge.1) then err1 = 1.0 endif if(select_mech.eq.1)then call isotropic(i_fix,Temp_wvk,Temp_ene,Temp_valley_i,Temp_valley_f) ! i_fix... mechanism, elseif(select_mech.eq.2)then call polar_optical_angle(i_fix,Temp_wvk,Temp_ene,Temp_valley_i) ! elseif(select_mech.eq.3)then ! call Coulomb_angle_BH() endif Temp_valley_i = Temp_valley_f ! 222 continue ! self scattering return end !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! SCATTERING SUBROUTINES (CHANGE ENERGY AND WAVEVECTORS OF PARTICLES)In the definition of the scattering rates, a variable called 'flag_mech' has been defined. ! The values assigned to this variable correspond to: 1: Isotropic scattering (acoustic, intervalley) ! 2: Polar_optical ==> anisotropic scattering (small angle) 3: Coulomb scattering ==> small-angle scattering !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ISOTROPIC SCATTERING PROCESS - uniform probability density for scattering in all directions, Acoustic, intervalley !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine isotropic(i_fix,Temp_wvk,Temp_ene,Temp_valley_i,Temp_valley_f) use particles use scattering use simulation real ( kind = 8 ), dimension(DIM) :: Temp_wvk ! carrier position and wave vector k integer :: Temp_valley_i, Temp_valley_f real ( kind = 8 ) :: Temp_ene !, Temp_tau integer i_fix real ( kind = 8 ) :: fi, ct, st, rknew ! Test real ( kind = 8 ) :: Test_ene, err1 !, Temp_tau real ( kind = 8 ):: sq_wvk_sum, abs_wvk, factor integer j ! Update carrier energy Temp_ene = Temp_ene + E_change(i_fix,Temp_valley_i) ! Update carrier wavevector rknew = semh(Temp_valley_f)*sqrt(Temp_ene*(1.+nonParabolPara(Temp_valley_f)*Temp_ene)) call random_number(harvest = RandomNum) fi = 2.*pi*RandomNum call random_number(harvest = RandomNum) ct = 1. - 2.*RandomNum st = sqrt(1.-ct*ct) ! isotropic random... Temp_wvk(1) = rknew*st*cos(fi) Temp_wvk(2) = rknew*st*sin(fi) Temp_wvk(3) = rknew*ct !Test sq_wvk_sum = 0. do j = 1,DIM !col_len = RCutoff_SiIn sq_wvk_sum = sq_wvk_sum + Temp_wvk(j)*Temp_wvk(j) enddo abs_wvk = sqrt(sq_wvk_sum) factor = hhem(Temp_valley_f)*sq_wvk_sum Test_ene = 2*factor/(1.+sqrt(1.+nonParabolPara4(Temp_valley_f)*factor)) if (Test_ene.ne.Temp_ene) then err1 = 1.0 endif if (Test_ene.ge.1) then err1 = 1.0 endif return end !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! POLAR OPTICAL PHONONS SCATTERING ANGLE 0- Randomize the polar angle according to the notes !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine polar_optical_angle(i_fix,Temp_wvk,Temp_ene,Temp_valley_i) use particles use scattering use simulation integer i_fix real ( kind = 8 ), dimension(DIM) :: Temp_wvk ! carrier position and wave vector k integer :: Temp_valley_i real ( kind = 8 ) :: Temp_ene !, Temp_tau real ( kind = 8 ) :: k,kxy,kxp,kyp,kzp,kp real ( kind = 8 ) :: enew, cth0, sth0, cfi0, sfi0, ge, gnew, zeta, cth, sth, cfi, sfi ! Test real ( kind = 8 ) :: Test_ene, err1 !, Temp_tau real ( kind = 8 ):: sq_wvk_sum, abs_wvk, factor integer j ! Update carrier energy enew = Temp_ene + E_change(i_fix,Temp_valley_i) ! Calculate the rotation angles kxy = sqrt(Temp_wvk(1)*Temp_wvk(1)+Temp_wvk(2)*Temp_wvk(2)) k = sqrt(Temp_wvk(1)*Temp_wvk(1)+Temp_wvk(2)*Temp_wvk(2)+Temp_wvk(3)*Temp_wvk(3)) ! total k cth0 = Temp_wvk(3)/k ! cos theta = kz/k sth0 = kxy/k ! sin theta = kxy/k cfi0 = Temp_wvk(1)/kxy ! cos phi.. = kx/kxy sfi0 = Temp_wvk(2)/kxy ! sin phi.. = ky/kxy ! Randomize momentum in the rotated coordinate system kp = semh(Temp_valley_i)*sqrt(enew*(1.+nonParabolPara(Temp_valley_i)*enew)) ge = Temp_ene*(1.+nonParabolPara(Temp_valley_i)*Temp_ene) gnew = enew*(1.+nonParabolPara(Temp_valley_i)*enew) zeta = 2.*sqrt(ge*gnew)/(ge+gnew-2.*sqrt(ge*gnew)) call random_number(harvest = RandomNum) cth = ((zeta+1.)-(2.*zeta+1)**RandomNum)/zeta sth = sqrt(1.-cth*cth) call random_number(harvest = RandomNum) fi = 2.*pi*RandomNum cfi = cos(fi) sfi = sin(fi) kxp = kp*sth*cfi kyp = kp*sth*sfi kzp = kp*cth ! Return back to the original coordinate system Temp_wvk(1) = kxp*cfi0*cth0-kyp*sfi0+kzp*cfi0*sth0 Temp_wvk(2) = kxp*sfi0*cth0+kyp*cfi0+kzp*sfi0*sth0 Temp_wvk(3) = -kxp*sth0+kzp*cth0 Temp_ene = enew !Test sq_wvk_sum = 0. do j = 1,DIM !col_len = RCutoff_SiIn sq_wvk_sum = sq_wvk_sum + Temp_wvk(j)*Temp_wvk(j) enddo abs_wvk = sqrt(sq_wvk_sum) factor = hhem(Temp_valley_i)*sq_wvk_sum Test_ene = 2*factor/(1.+sqrt(1.+nonParabolPara4(Temp_valley_i)*factor)) if (Test_ene.ne.Temp_ene) then err1 = 1.0 endif if (Test_ene.ge.1) then err1 = 1.0 endif return end !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! COULOMB SCATTERING ANGLE - Randomize the polar angle according to the expressions derived in the notes ! The assumption is that the scattering process is elastic !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !subroutine Coulomb_angle_BH() ! integer iv,i_fix ! real kx,ky,kz,e ! real k,kxy ! real kxp,kyp,kzp !! Update carrier energy !! enew = e + w(i_fix,iv) ! Calculate the rotation angles ! kxy = sqrt(kx*kx+ky*ky) ! k = sqrt(kxy*kxy+kz*kz) ! cth0 = kz/k ! sth0 = kxy/k ! cfi0 = kx/kxy ! sfi0 = ky/kxy !! Randomize momentum in the rotated coordinate system ! ge = e*(1.+af(iv)*e) ! rr = ran(iso) ! cth = 1. - 2*rr/(1.+4.*(1-rr)*ge/Energy_debye) ! sth = sqrt(1.-cth*cth) ! fi = two_pi*ran(iso) ! cfi = cos(fi) ! sfi = sin(fi) ! kxp = k*sth*cfi ! kyp = k*sth*sfi ! kzp = k*cth !! Return back to the original coordinate system ! kx = kxp*cfi0*cth0-kyp*sfi0+kzp*cfi0*sth0 ! ky = kxp*sfi0*cth0+kyp*cfi0+kzp*sfi0*sth0 ! kz = -kxp*sth0+kzp*cth0 ! e = enew ! return !end