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