Photon distribution (nonmonochromatic FEL) + bugfix

Source_files/Little_subroutines.f90

@@ -585,6 +585,29 @@ pure subroutine Gaussian(mu, sigma, x, Gaus, normalized_max) ! at the time x acc
 end subroutine Gaussian
 
 
+function sample_gaussian(x0, sigma, if_FWHM) result(x_out)
+   real(8) x_out
+   real(8), intent(in) :: x0, sigma ! gaussian mean and sigma (by default, it's not FWHM)
+   logical, intent(in), optional :: if_FWHM   ! if given sigma is FWHM (instead of standard deviation)
+   !---------------------
+   real(8), parameter :: Pi = 3.1415926535897932384626433832795d0
+   real(8) :: rvt1, rvt2, t0, s0
+
+   s0 = sigma
+   if (present(if_FWHM)) then ! check if it is FWHM instead of sigma
+      if (if_FWHM) s0 = sigma/( 2.0d0*sqrt( 2.0d0*log(2.0d0) ) )
+   endif
+
+   t0 = -10.0d0 * s0 ! to start with
+   do while (abs(t0) > 3.0e0*s0) ! accept only within 3*sigma interval
+      call random_number(rvt1)   ! Random value 1 for generated time t0
+      call random_number(rvt2)   ! Random value 2 for generated time t0
+      t0 = s0 * COS(2.0d0*Pi*rvt2)*SQRT(-2.0d0*LOG(rvt1))    ! Gaussian distribution for t0 around 0
+      x_out = x0 + t0   ! around the given mean value x0
+   enddo
+end function sample_gaussian
+
+
 
 subroutine interpolate_data_on_grid(given_array_x, given_array_y, given_grid, array_to_fill)
    real(8), dimension(:), intent(in) :: given_array_x, given_array_y, given_grid

Source_files/Monte_Carlo.f90

@@ -647,8 +647,8 @@ subroutine MC_for_photon(tim, MC, Nph, numpar, matter, laser, Scell, Eetot_cur,
          call random_number(RN)
          MC%photons(j)%ti = tim + RN*numpar%dt ! [fs] time of photoabsorbtion within this timestep
          t_cur = MC%photons(j)%ti	! [fs] time of photoabsorbtion
-         ! find number of pulse to which this photon belongs:
-         call Find_in_1D_array(laser(:)%hw, MC%photons(j)%E, PN) ! module 'Little_subroutines'
+         !! find number of pulse to which this photon belongs:
+         !call Find_in_1D_array(laser(:)%hw, MC%photons(j)%E, PN) ! module 'Little_subroutines'
          !call Photon_which_shell(PN, matter, shl) ! finds shell from which electron is ionized
          ! Find by which shell this photon is absorbed:
          call which_shell(MC%photons(j)%E, matter%Atoms, matter, 0, KOA, SHL) ! module 'MC_cross_sections'
@@ -747,7 +747,14 @@ subroutine choose_photon_energy(numpar, laser, MC, tim) ! see below
          coun = coun + 1 ! next pulse
          sum_phot = sum_phot + Nphot_pulses(coun) ! total number of photons
       enddo
-      MC%photons(i)%E = laser(coun)%hw	! [eV] photon energy in this pulse #coun
+
+      ! Distribution of the photon energy:
+      if (laser(coun)%FWHM_hw > 0.0d0) then  ! if it is a distribution
+         MC%photons(i)%E = sample_gaussian(laser(coun)%hw, laser(coun)%FWHM_hw, .true.) ! module "Little_subroutine"
+      else  ! single photon energy
+         MC%photons(i)%E = laser(coun)%hw ! [eV] photon energy in this pulse #coun
+      endif
+!       print*, MC%photons(i)%E
    enddo
 end subroutine choose_photon_energy
 

Source_files/Objects.f90

@@ -625,6 +625,7 @@ MODULE Objects
 type Pulse
    integer :: KOP	! kind of pulse: 0 = flat-top, 1 = Gaussian, 2 = SASE
    real(8) :: hw	! [eV] photon energy
+   real(8) :: FWHM_hw   ! [eV] distribution of photon energy spectrum (assumed gaussian)
    real(8) :: t		! [fs] pulse duration
    real(8) :: t0	! [fs] pulse maximum position
    real(8) :: F		! [eV/atom] absorbed fluence

Source_files/Read_input_data.f90

@@ -204,6 +204,7 @@ subroutine initialize_default_single_laser(laser, i) ! Must be already allocated
    integer, intent(in) :: i	! number of the pulse
    laser(i)%F = 0.0d0  ! ABSORBED DOSE IN [eV/atom]
    laser(i)%hw = 100.0d0  ! PHOTON ENERGY IN [eV]
+   laser(i)%FWHM_hw = 0.0d0  ! distribution of photon energies [eV]
    laser(i)%t = 10.0d0	  ! PULSE FWHM-DURATION IN
    laser(i)%KOP = 1  	  ! type of pulse: 0=rectangular, 1=Gaussian, 2=SASE
    !laser(i)%t = laser(i)%t/2.35482	! make a gaussian parameter out of it
@@ -4807,8 +4808,13 @@ subroutine read_input_material(File_name, Scell, matter, numpar, laser, Err)
          endif
       endif
 
-      read(FN,*,IOSTAT=Reason) laser(i)%hw  ! PHOTON ENERGY IN [eV]
-      call read_file(Reason, count_lines, read_well)
+      read(FN,*,IOSTAT=Reason) laser(i)%hw, laser(i)%FWHM_hw  ! photon energy in [eV], and FWHM width of distribution [eV]
+      if (Reason /= 0) then ! probably only energy provided, not the distribution width, so assume FWHM_hw=0
+         laser(i)%FWHM_hw = 0.0d0
+         backspace(FN)  ! get back and try to read the line again:
+         read(FN,*,IOSTAT=Reason) laser(i)%hw  ! photon energy in [eV] only, no distribution
+         call read_file(Reason, count_lines, read_well)
+      endif
       if (.not. read_well) then
          write(Error_descript,'(a,i5,a,$)') 'Could not read line ', count_lines, ' in file '//trim(adjustl(File_name))
          call Save_error_details(Err, 3, Error_descript)

Source_files/TB.f90

@@ -1303,7 +1303,7 @@ subroutine get_DOS_masks(Scell, matter, numpar, only_coupling, do_cartesian)
    nat = size(Scell(NSC)%MDatoms) ! number of atoms
    Nsiz = size(Scell(NSC)%Ha,1) ! total number of orbitals
 
-   BS:if (do_cart) then ! Cartesian basis set:
+   BS:if (do_cart) then ! Cartesian basis set (UNUSED FOR NOW):
       ! Find number of different orbital types:
       norb = identify_xTB_orbitals_per_atom(numpar%N_basis_size) ! module "TB_xTB"
       n_types = number_of_types_of_orbitals(norb, cartesian=.true.)  ! module "Little_subroutines"