1subroutine build_spectrum(ampl,en,ipol) 2!this subroutine builds up the absorption spectrum 3!and prints it on file 4 5USE exciton 6USE io_global, ONLY : stdout,ionode 7USE bse_wannier, ONLY : spectra_e_min,spectra_e_max,n_eig,spectra_nstep,spectra_broad 8USE constants, ONLY : RYTOEV, PI 9USE cell_base, ONLY : omega 10USE io_files, ONLY : tmp_dir,prefix 11 12implicit none 13 14REAL(DP), INTENT(in) :: ampl(n_eig), en(n_eig) 15 16INTEGER, EXTERNAL :: find_free_unit 17 18REAL(DP), ALLOCATABLE :: omega_g(:),absorption(:),broad_abs(:),excdos(:) 19COMPLEX(DP), ALLOCATABLE :: den(:,:),cspectrum(:),campl(:) 20REAL(DP) :: eta,step,prefac,sumdos,norm 21COMPLEX(DP) :: lambda_sum 22INTEGER :: i,j,iun,ipol 23 24LOGICAL :: debug 25 26call start_clock('build_spectrum') 27debug=.true. 28eta=0.001 29 30if(debug) then 31 if(ionode) then 32 do i=1,n_eig 33 write(stdout,*) '#',i,'E=',en(i),'A=',ampl(i) 34 enddo 35 endif 36endif 37 38allocate(omega_g(spectra_nstep)) 39allocate(absorption(spectra_nstep)) 40allocate(excdos(spectra_nstep)) 41allocate(broad_abs(spectra_nstep)) 42allocate(cspectrum(spectra_nstep)) 43allocate(den(spectra_nstep,n_eig)) 44allocate(campl(n_eig)) 45 46if (ipol==1) then 47!convert energy range in Ry 48 spectra_e_min=spectra_e_min/RYTOEV 49 spectra_e_max=spectra_e_max/RYTOEV 50endif 51 52!build the omega grid 53step=(spectra_e_max-spectra_e_min)/dble(spectra_nstep-1) 54 55do i=0, spectra_nstep-1 56 den(i+1,1:n_eig)=dcmplx(1.d0,0.d0)/(dcmplx(en(1:n_eig),0.d0)-dcmplx(spectra_e_min+dble(i)*step,eta)) 57 omega_g(i+1)=(spectra_e_min+dble(i)*step)*RYTOEV 58enddo 59 60!compute the absorption spectrum 61campl(1:n_eig)=dcmplx(ampl(1:n_eig),0.d0) 62!campl(1:n_eig)=dcmplx(1.d0,0.d0) 63 64cspectrum(1:spectra_nstep)=(0.d0,0.d0) 65 66call zgemm('N','N',spectra_nstep,1,n_eig,(1.d0,0.d0),den,spectra_nstep,campl,n_eig,(0.d0,0.d0),cspectrum,spectra_nstep) 67 68prefac=8.d0*PI/omega 69!prefac=1.d0 70absorption(1:spectra_nstep)=prefac*aimag(cspectrum(1:spectra_nstep)) 71 72!add gaussian broadening 73broad_abs(1:spectra_nstep)=0.d0 74do i=1,spectra_nstep 75 norm=0.d0 76 do j=1,spectra_nstep 77 broad_abs(i)=broad_abs(i)+& 78 absorption(j)*exp(-((omega_g(i)-omega_g(j))**2)/(2.d0*spectra_broad**2)) 79 norm=norm+exp(-((omega_g(i)-omega_g(j))**2)/(2.d0*spectra_broad**2)) 80 enddo 81 broad_abs(i)=broad_abs(i)/norm 82enddo 83 84!compute DOS using a lorentzian 85if (ipol==1) then 86 excdos(1:spectra_nstep)=0.d0 87 do i=0,spectra_nstep-1 88 do j=1,n_eig 89 excdos(i+1)=excdos(i+1)+2.d0*eta/(PI*((en(j)-spectra_e_min-dble(i)*step)**2+eta**2)) 90 enddo 91 enddo 92 excdos(1:spectra_nstep)=excdos(1:spectra_nstep)/(2.d0*n_eig) 93endif 94 95 96 97write(*,*) 'Absorption' 98write(*,*) 'Energy(eV) Eps2 Eps2(Nogaussbroad)' 99 100!print absorption spectrum on file 101if(ionode) then 102 iun = find_free_unit() 103 104 if (ipol==1) then 105 open(unit=iun, file=trim(tmp_dir)//trim(prefix)//'.eps2x.dat', status='unknown', form='formatted') 106 elseif (ipol==2) then 107 open(unit=iun, file=trim(tmp_dir)//trim(prefix)//'.eps2y.dat', status='unknown', form='formatted') 108 elseif (ipol==3) then 109 open(unit=iun, file=trim(tmp_dir)//trim(prefix)//'.eps2z.dat', status='unknown', form='formatted') 110 endif 111 112 do i=1,spectra_nstep 113 write(iun,*) omega_g(i),broad_abs(i), absorption (i) 114 write(*,*) omega_g(i),broad_abs(i), absorption(i) 115 enddo 116 close(iun) 117endif 118 119!print excdos spectrum on file 120if(ionode.and.(ipol==1)) then 121 sumdos=0.d0 122 iun = find_free_unit() 123 open(unit=iun, file=trim(tmp_dir)//trim(prefix)//'.excdos.dat', status='unknown', form='formatted') 124 do i=1,spectra_nstep 125 write(iun,*) omega_g(i),excdos(i) 126 sumdos=sumdos+excdos(i) 127 enddo 128 close(iun) 129 write(*,*) 'sumdos=',sumdos/dble(spectra_nstep) 130endif 131 132deallocate(omega_g,absorption,broad_abs,cspectrum,den,campl,excdos) 133 134call stop_clock('build_spectrum') 135end subroutine 136 137