Line data Source code
1 : !!****m* ABINIT/m_optic_tools
2 : !! NAME
3 : !! m_optic_tools
4 : !!
5 : !! FUNCTION
6 : !! Helper functions used in the optic code
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2002-2026 ABINIT group (SSharma,MVer,VRecoules,TD,YG, NAP,VT)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
14 : !!
15 : !! COMMENTS
16 : !!
17 : !! Right now the routine sums over the k-points. In future linear tetrahedron method might be useful.
18 : !!
19 : !! Reference articles:
20 : !!
21 : !! 1. S. Sharma, J. K. Dewhurst and C. Ambrosch-Draxl, Phys. Rev. B {\bf 67} 165332 2003 [[cite:Sharma2003]]
22 : !! 2. J. L. P. Hughes and J. E. Sipe, Phys. Rev. B {\bf 53} 10 751 1996 [[cite:Hughes1996]]
23 : !! 3. S. Sharma and C. Ambrosch-Draxl, Physica Scripta T 109 2004 [[cite:Sharma2004]]
24 : !! 4. J. E. Sipe and Ed. Ghahramani, Phys. Rev. B {\bf 48} 11 705 1993 [[cite:Sipe1993]]
25 : !!
26 : !! SOURCE
27 :
28 : #if defined HAVE_CONFIG_H
29 : #include "config.h"
30 : #endif
31 :
32 : #include "abi_common.h"
33 :
34 : module m_optic_tools
35 :
36 : use defs_basis
37 : use m_errors
38 : use m_abicore
39 : use m_linalg_interfaces
40 : use m_xmpi
41 : use m_nctk
42 : use netcdf
43 : use m_ebands
44 :
45 : use m_numeric_tools, only : c2r
46 : use m_io_tools, only : flush_unit, open_file
47 : use m_crystal, only : crystal_t
48 :
49 : implicit none
50 :
51 : private
52 :
53 : public :: pmat2cart
54 : public :: pmat_renorm
55 : public :: linopt ! Compute dielectric function for semiconductors
56 : public :: nlinopt ! Second harmonic generation susceptibility for semiconductors
57 : public :: linelop ! Linear electro-optic susceptibility for semiconductors
58 : public :: nonlinopt ! nonlinear electro-optic susceptibility for semiconductors
59 :
60 : contains
61 : !!***
62 :
63 : !----------------------------------------------------------------------
64 :
65 : !!****f* m_optic_tools/pmat2cart
66 : !! NAME
67 : !! pmat2cart
68 : !!
69 : !! FUNCTION
70 : !! turn momentum matrix elements to cartesian axes. To be used in optic calculation of linear
71 : !! and non-linear RPA dielectric matrices
72 : !!
73 : !! INPUTS
74 : !! eigen11,eigen12,eigen13 = first order ddk eigen values = d eig_i,k / dk for 3 reduced directions
75 : !! mband=maximum number of bands
76 : !! nkpt = number of k-points
77 : !! nsppol=1 for unpolarized, 2 for spin-polarized
78 : !! rprimd(3,3)=dimensional real space primitive translations (bohr)
79 : !!
80 : !! OUTPUT
81 : !! pmat(mband,mband,nkpt,3,nsppol) = matrix elements of momentum operator, in cartesian coordinates
82 : !!
83 : !! SOURCE
84 :
85 14 : subroutine pmat2cart(eigen11, eigen12, eigen13, mband, nkpt, nsppol, pmat, rprimd)
86 :
87 : !Arguments -----------------------------------------------
88 : !scalars
89 : integer,intent(in) :: mband,nkpt,nsppol
90 : !arrays
91 : real(dp),intent(in) :: eigen11(2,mband,mband,nkpt,nsppol)
92 : real(dp),intent(in) :: eigen12(2,mband,mband,nkpt,nsppol)
93 : real(dp),intent(in) :: eigen13(2,mband,mband,nkpt,nsppol),rprimd(3,3)
94 : !no_abirules
95 : complex(dp),intent(out) :: pmat(mband,mband,nkpt,3,nsppol)
96 :
97 : !Local variables -----------------------------------------
98 : !scalars
99 : integer :: iband1,iband2,ikpt,isppol
100 : !arrays
101 : real(dp) :: rprim(3,3)
102 : ! *************************************************************************
103 :
104 : !rescale the rprim
105 182 : rprim(:,:) = rprimd(:,:) / two_pi
106 :
107 32 : do isppol=1,nsppol
108 984 : do ikpt=1,nkpt
109 17978 : do iband1=1,mband
110 396232 : do iband2=1,mband
111 : pmat(iband2,iband1,ikpt,:,isppol) = &
112 : rprim(:,1)*cmplx(eigen11(1,iband2,iband1,ikpt,isppol),eigen11(2,iband2,iband1,ikpt,isppol),kind=dp) &
113 : +rprim(:,2)*cmplx(eigen12(1,iband2,iband1,ikpt,isppol),eigen12(2,iband2,iband1,ikpt,isppol),kind=dp) &
114 1530096 : +rprim(:,3)*cmplx(eigen13(1,iband2,iband1,ikpt,isppol),eigen13(2,iband2,iband1,ikpt,isppol),kind=dp)
115 : end do
116 : end do
117 : end do
118 : end do
119 :
120 14 : end subroutine pmat2cart
121 : !!***
122 :
123 : !----------------------------------------------------------------------
124 :
125 : !!****f* m_optic_tools/pmat_renorm
126 : !! NAME
127 : !! pmat_renorm
128 : !!
129 : !! FUNCTION
130 : !! Renormalize the momentum matrix elements according to the scissor shift which is imposed
131 : !!
132 : !! INPUTS
133 : !! mband= number of bands
134 : !! nkpt = number of k-points
135 : !! nsppol=1 for unpolarized, 2 for spin-polarized
136 : !! fermie = Fermi level
137 : !! sc = scissor shift for conduction bands
138 : !! eig = ground state eigenvalues
139 : !!
140 : !! OUTPUT
141 : !! pmat(mband,mband,nkpt,3,nsppol) = momentum matrix elements, renormalized by denominator change with scissor shift
142 : !!
143 : !! SOURCE
144 :
145 14 : subroutine pmat_renorm(fermie, eig, mband, nkpt, nsppol, pmat, sc)
146 :
147 : !Arguments -----------------------------------------------
148 : !scalars
149 : integer, intent(in) :: nsppol
150 : integer, intent(in) :: nkpt
151 : integer, intent(in) :: mband
152 : real(dp), intent(in) :: fermie
153 : real(dp), intent(in) :: sc
154 : !arrays
155 : real(dp), intent(in) :: eig(mband,nkpt,nsppol)
156 : complex(dp), intent(inout) :: pmat(mband,mband,nkpt,3,nsppol)
157 :
158 : !Local variables -----------------------------------------
159 : !scalars
160 : integer :: iband1,iband2,ikpt,isppol
161 : real(dp) :: corec, e1, e2
162 : ! *************************************************************************
163 :
164 14 : if (abs(sc) < tol8) then
165 12 : call wrtout(std_out,' No scissor shift to be applied. Returning to main optic routine.',"COLL")
166 12 : return
167 : end if
168 :
169 4 : do isppol=1,nsppol
170 140 : do ikpt=1,nkpt
171 1402 : do iband1=1,mband ! valence states
172 1264 : e1 = eig(iband1,ikpt,isppol)
173 1264 : if (e1 > fermie) cycle
174 5736 : do iband2=1,mband ! conduction states
175 5056 : e2 = eig(iband2,ikpt,isppol)
176 5056 : if (e2 < fermie) cycle
177 2880 : corec = (e2+sc-e1)/(e2-e1)
178 11520 : pmat(iband2,iband1,ikpt,:,isppol) = corec * pmat(iband2,iband1,ikpt,:,isppol)
179 12784 : pmat(iband1,iband2,ikpt,:,isppol) = corec * pmat(iband1,iband2,ikpt,:,isppol)
180 : end do
181 : end do
182 : end do
183 : end do
184 :
185 : end subroutine pmat_renorm
186 : !!***
187 :
188 : !----------------------------------------------------------------------
189 :
190 : !!****f* m_optic_tools/linopt
191 : !! NAME
192 : !! linopt
193 : !!
194 : !! FUNCTION
195 : !! Compute optical frequency dependent dielectric function for semiconductors
196 : !!
197 : !! INPUTS
198 : !! icomp=Sequential index associated to computed tensor components (used for netcdf output)
199 : !! itemp=Temperature index (used for netcdf output)
200 : !! nband_sum=Number of bands included in the sum. Must be <= mband
201 : !! pmat(mband,mband,nkpt,3,nsppol)=momentum matrix elements in cartesian coordinates(complex)
202 : !! v1,v2=desired component of the dielectric function(integer) 1=x,2=y,3=z
203 : !! nmesh=desired number of energy mesh points(integer)
204 : !! de=desired step in energy(real); nmesh*de=maximum energy
205 : !! sc=scissors shift in Ha(real)
206 : !! brod=broadening in Ha(real)
207 : !! fnam=root for filename that will contain the output filename will be trim(fnam)//'-linopt.out'
208 : !! ncid=Netcdf id to save output data.
209 : !! prtlincompmatrixelements=if set to 1, the matrix elements are dumped in the _OPTIC.nc file for post processing.
210 : !!
211 : !! SIDE EFFECTS
212 : !! Dielectric function for semiconductors, on a desired energy mesh and for a desired
213 : !! direction of polarisation is written to file.
214 : !! The output is in a file named trim(fnam)//'-linopt.out' and contains
215 : !! Im(\epsilon_{v1v2}(\omega), Re(\epsilon_{v1v2}(\omega) and abs(\epsilon_{v1v2}(\omega).
216 : !!
217 : !! If 'prtlincompmatrixelements' is set to 1, the matrix elements and other quantities used to build
218 : !! the chi tensor are stored in the _OPTIC.nc file as well. This includes the matrix elements,
219 : !! the occupations, the renormalized but unshifted eigenvalues and the kpts weights.
220 : !!
221 : !! Comment:
222 : !! Right now the routine sums over the kpoints. In future linear tetrahedron method should be useful.
223 : !!
224 : !! SOURCE
225 :
226 34 : subroutine linopt(icomp, itemp, nband_sum, cryst, ks_ebands, EPBSt, pmat, &
227 : v1, v2, nmesh, de, sc, brod, fnam, ncid, prtlincompmatrixelements, comm)
228 :
229 : !Arguments ------------------------------------
230 : integer, intent(in) :: icomp,itemp,nband_sum, ncid
231 : type(crystal_t), intent(in) :: cryst
232 : type(ebands_t),intent(in) :: ks_ebands,EPBSt
233 : complex(dp), intent(in) :: pmat(ks_ebands%mband, ks_ebands%mband, ks_ebands%nkpt, 3, ks_ebands%nsppol)
234 : integer, intent(in) :: v1, v2, nmesh
235 : real(dp), intent(in) :: de, sc, brod
236 : character(len=*), intent(in) :: fnam
237 : integer, intent(in) :: comm
238 : integer, intent(in) :: prtlincompmatrixelements
239 :
240 : !Local variables -------------------------
241 : integer,parameter :: master=0
242 : integer :: isp,i,j,isym,lx,ly,ik,ist1,ist2,iw,nkpt
243 : integer :: my_rank, nproc, my_k1, my_k2, ierr, fout1, mband, nsppol
244 : integer :: ncerr
245 : logical :: do_linewidth
246 : real(dp) :: deltav1v2, tmpabs, renorm_factor,emin,emax
247 : real(dp) :: ene,abs_eps,re_eps
248 : complex(dp) :: e1,e2,e12, e1_ep,e2_ep,e12_ep, b11,b12, ieta, w
249 : character(len=fnlen) :: fnam1
250 : character(len=500) :: msg
251 : ! allocatable arrays
252 : real(dp) :: s(3,3),sym(3,3)
253 34 : real(dp), allocatable :: im_refract(:),re_refract(:)
254 34 : complex(dp), allocatable :: chi(:,:), matrix_elements(:,:,:,:), renorm_eigs(:,:,:), eps(:)
255 :
256 : ! *********************************************************************
257 :
258 34 : my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
259 34 : nkpt = ks_ebands%nkpt
260 34 : nsppol = ks_ebands%nsppol
261 34 : mband = ks_ebands%mband
262 34 : ABI_CHECK(nband_sum <= mband, "nband_sum <= mband")
263 :
264 34 : if (my_rank == master) then
265 : ! check polarisation
266 34 : if (v1.le.0.or.v2.le.0.or.v1.gt.3.or.v2.gt.3) then
267 0 : write(std_out,*) '---------------------------------------------'
268 0 : write(std_out,*) ' Error in linopt: '
269 0 : write(std_out,*) ' the polarisation directions incorrect '
270 0 : write(std_out,*) ' 1=x and 2=y and 3=z '
271 0 : write(std_out,*) '---------------------------------------------'
272 0 : ABI_ERROR("Aborting now")
273 : end if
274 : ! number of energy mesh points
275 34 : if (nmesh.le.0) then
276 0 : write(std_out,*) '---------------------------------------------'
277 0 : write(std_out,*) ' Error in linopt: '
278 0 : write(std_out,*) ' number of energy mesh points incorrect '
279 0 : write(std_out,*) ' number has to integer greater than 0 '
280 0 : write(std_out,*) ' nmesh*de = max energy for calculation '
281 0 : write(std_out,*) '---------------------------------------------'
282 0 : ABI_ERROR("Aborting now")
283 : end if
284 : ! step in energy
285 34 : if (de.le.zero) then
286 0 : write(std_out,*) '---------------------------------------------'
287 0 : write(std_out,*) ' Error in linopt: '
288 0 : write(std_out,*) ' energy step is incorrect '
289 0 : write(std_out,*) ' number has to real greater than 0.0 '
290 0 : write(std_out,*) ' nmesh*de = max energy for calculation '
291 0 : write(std_out,*) '---------------------------------------------'
292 0 : ABI_ERROR("Aborting now")
293 : end if
294 : ! broadening
295 34 : if (brod.gt.0.009) then
296 0 : write(std_out,*) '---------------------------------------------'
297 0 : write(std_out,*) ' ATTENTION: broadening is quite high '
298 0 : write(std_out,*) ' ideally should be less than 0.005 '
299 0 : write(std_out,*) '---------------------------------------------'
300 : else if (brod.gt.0.015) then
301 : write(std_out,*) '----------------------------------------'
302 : write(std_out,*) ' ATTENTION: broadening is too high '
303 : write(std_out,*) ' ideally should be less than 0.005 '
304 : write(std_out,*) '----------------------------------------'
305 : end if
306 : ! fermi energy
307 34 : if(ks_ebands%fermie<-1.0d4) then
308 0 : write(std_out,*) '---------------------------------------------'
309 0 : write(std_out,*) ' ATTENTION: Fermi energy seems extremely '
310 0 : write(std_out,*) ' low '
311 0 : write(std_out,*) '---------------------------------------------'
312 : end if
313 : ! scissors operator
314 34 : if (sc.lt.zero) then
315 0 : write(std_out,*) '---------------------------------------------'
316 0 : write(std_out,*) ' Error in linopt: '
317 0 : write(std_out,*) ' scissors shift is incorrect '
318 0 : write(std_out,*) ' number has to be greater than 0.0 '
319 0 : write(std_out,*) '---------------------------------------------'
320 0 : ABI_ERROR("Aborting now")
321 : end if
322 : end if
323 :
324 34 : do_linewidth = allocated(EPBSt%linewidth)
325 : ! TODO: activate this, and remove do_linewidth - always add it in even if 0.
326 : ! if (.not. allocated(EPBSt%linewidth)) then
327 : ! ABI_MALLOC(EPBSt%linewidth, (1, mband, my_k2-my_k1+1, nsppol))
328 : ! EPBSt%linewidth = zero
329 : ! end if
330 :
331 136 : ABI_MALLOC(chi, (nmesh, nsppol))
332 102 : ABI_MALLOC(eps, (nmesh))
333 102 : ABI_MALLOC(im_refract, (nmesh))
334 68 : ABI_MALLOC(re_refract, (nmesh))
335 34 : ieta=(zero, 1._dp)*brod
336 34 : renorm_factor=1._dp/(cryst%ucvol*dble(cryst%nsym))
337 :
338 : ! output file names
339 34 : fnam1=trim(fnam)//'-linopt.out'
340 :
341 : ! construct symmetrisation tensor
342 34 : sym = zero
343 1210 : do isym=1,cryst%nsym
344 15288 : s(:,:)=cryst%symrel_cart(:,:,isym)
345 4738 : do i=1,3
346 15288 : do j=1,3
347 14112 : sym(i,j)=sym(i,j)+s(i,v1)*s(j,v2)
348 : end do
349 : end do
350 : end do
351 :
352 : ! calculate the energy window
353 34 : emin=zero
354 34 : emax=zero
355 1122 : do ik=1,nkpt
356 2402 : do isp=1,nsppol
357 29632 : do ist1=1,nband_sum
358 27264 : emin=min(emin,EPBSt%eig(ist1,ik,isp))
359 28544 : emax=max(emax,EPBSt%eig(ist1,ik,isp))
360 : end do
361 : end do
362 : end do
363 :
364 : ! Split work
365 34 : call xmpi_split_work(nkpt,comm,my_k1,my_k2)
366 : ! if we print matrix elements, allocate full arrays for each process
367 : ! this is not optimized memory-wise since we could just allocate what is needed
368 : ! however we would need to write all data using mpi-io.
369 34 : if (prtlincompmatrixelements == 1) then
370 119184 : ABI_CALLOC(matrix_elements, (mband, mband, nkpt, nsppol))
371 3982 : ABI_CALLOC(renorm_eigs, (mband, nkpt, nsppol))
372 : endif
373 :
374 : ! start calculating linear optical response
375 17276 : chi(:,:)=zero
376 76 : do isp=1,nsppol
377 1356 : do ik=my_k1,my_k2
378 1280 : write(std_out,*) "P-",my_rank,": ",ik,'of',nkpt
379 28586 : do ist1=1,nband_sum
380 27264 : e1=ks_ebands%eig(ist1,ik,isp)
381 27264 : e1_ep=EPBSt%eig(ist1,ik,isp)
382 : ! TODO: unless memory is a real issue, should set lifetimes to 0 and do this sum systematically
383 : ! instead of putting an if statement in a loop! See above
384 27264 : if(do_linewidth) then
385 1120 : e1_ep = e1_ep + EPBSt%linewidth(1,ist1,ik,isp)*(0.0_dp,1.0_dp)
386 : end if
387 719360 : do ist2=1,nband_sum
388 690816 : e2=ks_ebands%eig(ist2,ik,isp)
389 690816 : e2_ep=EPBSt%eig(ist2,ik,isp)
390 690816 : if(do_linewidth) then
391 15680 : e2_ep = e2_ep - EPBSt%linewidth(1,ist2,ik,isp)*(0.0_dp,1.0_dp)
392 : end if
393 718080 : if (ist1.ne.ist2) then
394 : ! scissors correction of momentum matrix
395 663552 : if(REAL(e1) > REAL(e2)) then
396 331776 : e12 = e1-e2+sc
397 : else
398 331776 : e12 = e1-e2-sc
399 : end if
400 663552 : if(REAL(e1_ep) > REAL(e2_ep)) then
401 331776 : e12_ep = e1_ep-e2_ep+sc
402 : else
403 331776 : e12_ep = e1_ep-e2_ep-sc
404 : end if
405 : ! e12=e1-e2-sc
406 663552 : b11=zero
407 : ! symmetrization of momentum matrix
408 2654208 : do lx=1,3
409 8626176 : do ly=1,3
410 : b11=b11+(sym(lx,ly)*pmat(ist1,ist2,ik,lx,isp)* &
411 7962624 : conjg(pmat(ist1,ist2,ik,ly,isp)))
412 : end do
413 : end do
414 663552 : b12=b11*renorm_factor*(1._dp/(e12**2))
415 : ! store data for printing if necessary
416 663552 : if (prtlincompmatrixelements == 1) then
417 111360 : matrix_elements(ist1,ist2,ik,isp) = b12
418 111360 : renorm_eigs(ist1,ik,isp) = e1_ep
419 111360 : renorm_eigs(ist2,ik,isp) = e2_ep
420 : endif
421 : ! calculate on the desired energy grid
422 204384000 : do iw=2,nmesh
423 203720448 : w=(iw-1)*de+ieta
424 : chi(iw,isp)=chi(iw,isp)+(ks_ebands%wtk(ik)*(ks_ebands%occ(ist1,ik,isp)-ks_ebands%occ(ist2,ik,isp))* &
425 204384000 : (b12/(-e12_ep-w)))
426 : end do ! frequencies
427 : end if
428 : end do ! states 2
429 : end do ! states 1
430 : end do ! k points
431 : end do ! spin
432 :
433 34 : call xmpi_sum(chi,comm,ierr)
434 34 : if (prtlincompmatrixelements == 1) then
435 : ! gather all data to main process in order to write them using a single process
436 : ! in the netcdf file. This could be avoided by doing mpiio.
437 2 : call xmpi_sum(matrix_elements,comm,ierr)
438 2 : call xmpi_sum(renorm_eigs,comm,ierr)
439 : endif
440 :
441 : ! calculate epsilon
442 34 : eps(1) = zero
443 34 : deltav1v2=zero; if (v1 == v2) deltav1v2=one
444 16000 : do iw=2,nmesh
445 33158 : eps(iw)=deltav1v2+four*pi*sum(chi(iw,:))
446 : end do
447 :
448 34 : if (my_rank == master) then
449 : ! open the output files
450 34 : if (open_file(fnam1,msg,newunit=fout1,action='WRITE',form='FORMATTED') /= 0) then
451 0 : ABI_ERROR(msg)
452 : end if
453 : ! write output
454 34 : write(fout1, '(a,2i3,a)' )' #calculated the component:',v1,v2,' of dielectric function'
455 34 : write(std_out,*) 'calculated the component:',v1,v2,' of dielectric function'
456 34 : write(fout1, '(a,2es16.6)' ) ' #broadening:', real(ieta),aimag(ieta)
457 34 : write(std_out,*) ' with broadening:',ieta
458 34 : write(fout1, '(a,es16.6)' ) ' #scissors shift:',sc
459 34 : write(std_out,*) 'and scissors shift:',sc
460 34 : write(fout1, '(a,es16.6,a,es16.6,a)' ) ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
461 34 : write(std_out,*) 'energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
462 34 : write(fout1,*)
463 34 : if(nsppol==1)write(fout1, '(a)' ) ' # Energy(eV) Im(eps(w))'
464 34 : if(nsppol==2)write(fout1, '(a)' ) ' # Energy(eV) Im(eps(w)) Spin up Spin down '
465 16000 : do iw=2,nmesh
466 15966 : ene=(iw-1)*de*Ha_eV
467 15966 : if(nsppol==1)write(fout1, '(2es16.6)' ) ene,aimag(eps(iw))
468 16000 : if(nsppol==2)write(fout1, '(4es16.6)' ) ene,aimag(eps(iw)),4._dp*pi*aimag(chi(iw,1)),4._dp*pi*aimag(chi(iw,2))
469 : end do
470 34 : write(fout1,*)
471 34 : write(fout1,*)
472 34 : if(nsppol==1)write(fout1, '(a)' ) ' # Energy(eV) Re(eps(w))'
473 34 : if(nsppol==2)write(fout1, '(a)' ) ' # Energy(eV) Re(eps(w)) Spin up Spin down +delta(diag) '
474 16000 : do iw=2,nmesh
475 15966 : ene=(iw-1)*de*Ha_eV
476 15966 : if(nsppol==1)write(fout1, '(2es16.6)' ) ene,dble(eps(iw))
477 16000 : if(nsppol==2)write(fout1, '(5es16.6)' ) ene,dble(eps(iw)),4._dp*pi*dble(chi(iw,1)),4._dp*pi*dble(chi(iw,2)),deltav1v2
478 : end do
479 34 : write(fout1,*)
480 34 : write(fout1,*)
481 34 : write(fout1, '(a)' )' # Energy(eV) abs(eps(w))'
482 16000 : do iw=2,nmesh
483 15966 : ene=(iw-1)*de*Ha_eV
484 15966 : abs_eps=abs(eps(iw))
485 15966 : re_eps=dble(eps(iw))
486 15966 : write(fout1, '(2es16.6)' ) ene,abs_eps
487 15966 : re_refract(iw)=sqrt(half*(abs_eps+re_eps))
488 16000 : im_refract(iw)=sqrt(half*(abs_eps-re_eps))
489 : end do
490 34 : write(fout1,*)
491 34 : write(fout1,*)
492 34 : write(fout1, '(a)' )' # Energy(eV) Im(refractive index(w)) aka kappa'
493 16000 : do iw=2,nmesh
494 15966 : ene=(iw-1)*de*Ha_eV
495 16000 : write(fout1, '(2es16.6)' ) ene,im_refract(iw)
496 : end do
497 34 : write(fout1,*)
498 34 : write(fout1,*)
499 34 : write(fout1, '(a)' )' # Energy(eV) Re(refractive index(w)) aka n'
500 16000 : do iw=2,nmesh
501 15966 : ene=(iw-1)*de*Ha_eV
502 16000 : write(fout1, '(2es16.6)' ) ene,re_refract(iw)
503 : end do
504 34 : write(fout1,*)
505 34 : write(fout1,*)
506 34 : write(fout1, '(a)' )' # Energy(eV) Reflectivity(w) from vacuum, at normal incidence'
507 16000 : do iw=2,nmesh
508 15966 : ene=(iw-1)*de*Ha_eV
509 16000 : write(fout1, '(2es16.6)' ) ene, ((re_refract(iw)-one)**2+im_refract(iw)**2)/((re_refract(iw)+one)**2+im_refract(iw)**2)
510 : end do
511 34 : write(fout1,*)
512 34 : write(fout1,*)
513 34 : write(fout1, '(a)' )' # Energy(eV) absorption coeff (in 10^6 m-1) = omega Im(eps) / c n(eps)'
514 16000 : do iw=2,nmesh
515 15966 : ene=(iw-1)*de
516 15966 : tmpabs=zero
517 15966 : if ( re_refract(iw) > tol10 ) then
518 14675 : tmpabs = aimag(eps(iw))*ene / re_refract(iw) / Speed_Light / Bohr_meter * 1.0d-6
519 : end if
520 16000 : write(fout1, '(2es16.6)' ) Ha_eV*ene, tmpabs
521 : end do
522 :
523 : ! close output file
524 34 : close(fout1)
525 :
526 34 : if (ncid /= nctk_noid) then
527 170 : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "linopt_epsilon"), c2r(eps), start=[1, 1, icomp, itemp])
528 34 : NCF_CHECK(ncerr)
529 : end if
530 34 : if (prtlincompmatrixelements == 1) then
531 : ! write matrix elements and other quantities used to build the chi tensor.
532 2 : write(std_out, '(a)') 'Writing linopt matrix elements in _OPTIC.nc file.'
533 : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "linopt_matrix_elements"), c2r(matrix_elements),&
534 16 : start=[1, 1, 1, 1, 1, icomp, itemp])
535 2 : NCF_CHECK(ncerr)
536 2 : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "linopt_renorm_eigs"), c2r(renorm_eigs), start=[1, 1, 1, 1])
537 2 : NCF_CHECK(ncerr)
538 :
539 : ! write occupations and kpt weights
540 2 : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "linopt_occupations"), ks_ebands%occ, start=[1, 1, 1])
541 2 : NCF_CHECK(ncerr)
542 2 : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "linopt_wkpts"), ks_ebands%wtk, start=[1])
543 2 : NCF_CHECK(ncerr)
544 2 : write(std_out, '(a)') 'Writing linopt matrix elements done.'
545 : endif
546 : end if ! rank == master
547 :
548 34 : ABI_FREE(chi)
549 34 : ABI_FREE(eps)
550 34 : ABI_FREE(im_refract)
551 34 : ABI_FREE(re_refract)
552 :
553 34 : ABI_SFREE(matrix_elements)
554 34 : ABI_SFREE(renorm_eigs)
555 :
556 102 : end subroutine linopt
557 : !!***
558 :
559 : !----------------------------------------------------------------------
560 :
561 : !!****f* m_optic_tools/nlinopt
562 : !! NAME
563 : !! nlinopt
564 : !!
565 : !! FUNCTION
566 : !! Compute second harmonic generation susceptibility for semiconductors
567 : !!
568 : !! INPUTS
569 : !! icomp=Sequential index associated to computed tensor components (used for netcdf output)
570 : !! itemp=Temperature index (used for netcdf output)
571 : !! nband_sum=Number of bands included in the sum. Must be <= mband
572 : !! fermie = Fermi energy in Ha(real)
573 : !! pmat(mband,mband,nkpt,3,nsppol) = momentum matrix elements in cartesian coordinates(complex)
574 : !! v1,v2,v3 = desired component of the dielectric function(integer) 1=x,2=y,3=z
575 : !! nmesh = desired number of energy mesh points(integer)
576 : !! de = desired step in energy(real); nmesh*de=maximum energy for plotting
577 : !! sc = scissors shift in Ha(real)
578 : !! brod = broadening in Ha(real)
579 : !! tol = tolerance:how close to the singularity exact exact is calculated(real)
580 : !! fnam=root for filenames that will contain the output :
581 : !! fnam1=trim(fnam)//'-ChiTotIm.out'
582 : !! fnam2=trim(fnam)//'-ChiTotRe.out'
583 : !! fnam3=trim(fnam)//'-ChiIm.out'
584 : !! fnam4=trim(fnam)//'-ChiRe.out'
585 : !! fnam5=trim(fnam)//'-ChiAbs.out'
586 : !! ncid=Netcdf id to save output data.
587 : !!
588 : !! OUTPUT
589 : !! Calculates the second harmonic generation susceptibility on a desired energy mesh and
590 : !! for desired direction of polarisation. The output is in files named
591 : !! ChiTot.out : Im\chi_{v1v2v3}(2\omega,\omega,-\omega) and Re\chi_{v1v2v3}(2\omega,\omega,-\omega)
592 : !! ChiIm.out : contributions to the Im\chi_{v1v2v3}(2\omega,\omega,-\omega) from various terms
593 : !! ChiRe.out : contributions to Re\chi_{v1v2v3}(2\omega,\omega,-\omega) from various terms
594 : !! ChiAbs.out : abs\chi_{v1v2v3}(2\omega,\omega,-\omega). The headers in these files contain
595 : !! information about the calculation.
596 : !! See eqs. (A4)-A(11) of S. Sharma et al Phys. Rev. B 67, 165332 (2003)
597 : !!
598 : !! SOURCE
599 :
600 8 : subroutine nlinopt(icomp, itemp, nband_sum, cryst, ks_ebands, pmat, &
601 : v1, v2, v3, nmesh, de, sc, brod, tol, w_decompo, fnam, contrib_decompo, do_decompo, do_antiresonant, &
602 : ncid, comm)
603 :
604 : !Arguments ------------------------------------
605 : integer, intent(in) :: icomp, itemp, nband_sum, ncid
606 : type(crystal_t),intent(in) :: cryst
607 : type(ebands_t),intent(in) :: ks_ebands
608 : complex(dp), intent(in) :: pmat(ks_ebands%mband, ks_ebands%mband, ks_ebands%nkpt, 3, ks_ebands%nsppol)
609 : integer, intent(in) :: v1, v2, v3, nmesh, comm
610 : real(dp), intent(in) :: de, sc, brod, tol
611 : character(len=*), intent(in) :: fnam
612 : logical, intent(in) :: do_decompo ! .TRUE. to perform the bands decomposition
613 : real(dp), intent(in) :: w_decompo ! energy in eV for bands decompo
614 : integer, intent(in) :: contrib_decompo ! 0=all, 12=inter2w, 22=intra2w, 11=inter1w, 21=intra1w, 1=intra1wS, 2=2bands
615 : logical, intent(in) :: do_antiresonant ! .FALSE. to include the AR terms, TRUE by default
616 :
617 : !Local variables -------------------------
618 : integer,parameter :: master=0
619 : integer :: iw, mband,i,j,k,lx,ly,lz
620 : integer :: isp,isym,ik,ist1,ist2,istl,istn,istm
621 : integer :: my_rank, nproc, my_k1, my_k2, ierr
622 : integer :: fout1,fout2,fout3,fout4,fout5,fout6,fout7
623 : real(dp) :: f1,f2,f3
624 : real(dp) :: ene,totre,totabs,totim
625 : real(dp) :: el,en,em,emin,emax,my_emin,my_emax
626 : real(dp) :: const_esu,const_au,au2esu,wmn,wnm,wln,wnl,wml,wlm, t1
627 : complex(dp) :: idel,w,zi
628 : complex(dp) :: mat2w,mat1w1,mat1w2,mat2w_tra,mat1w3_tra
629 : complex(dp) :: b111,b121,b131,b112,b122,b132,b113,b123,b133
630 : complex(dp) :: b241,b242,b243,b221,b222,b223,b211,b212,b213,b231
631 : complex(dp) :: b311,b312,b313,b331
632 : complex(dp) :: b24,b21_22,b11,b12_13,b31_32
633 : character(len=fnlen) :: fnam1,fnam2,fnam3,fnam4,fnam5,fnam6,fnam7
634 : character(500) :: msg
635 : ! local allocatable arrays
636 : integer :: start4(4),count4(4)
637 : real(dp) :: s(3,3),sym(3,3,3)
638 8 : complex(dp), allocatable :: px(:,:,:,:,:), py(:,:,:,:,:), pz(:,:,:,:,:)
639 8 : complex(dp), allocatable :: delta(:,:,:), inter2w(:), inter1w(:)
640 8 : complex(dp), allocatable :: intra2w(:), intra1w(:), intra1wS(:),chi2tot(:)
641 : ! Addition antiresonant (AR)
642 : ! Products of momentum matrix elem for the AR terms
643 : real(dp) :: t2
644 : complex(dp) :: mat2wa,mat1w1a,mat1w2a,mat2wa_tra,mat1w3a_tra
645 : complex(dp) :: a111,a112,a113,a11 ! AR inter2w term
646 8 : complex(dp), allocatable :: inter2wa(:)
647 : complex(dp) :: a121,a131,a122,a132,a123,a133,a12_13 ! AR inter1w term
648 8 : complex(dp), allocatable :: inter1wa(:)
649 : complex(dp) :: a241,a242,a243,a231,a24 ! AR intra2w term
650 8 : complex(dp), allocatable :: intra2wa(:)
651 : complex(dp) :: a211,a221,a212,a222,a213,a223,a21_22 ! AR intra1w term
652 8 : complex(dp), allocatable :: intra1wa(:)
653 : complex(dp) :: a311,a312,a313,a331,a31_32 ! AR intra1wS term
654 8 : complex(dp), allocatable :: intra1wSa(:)
655 8 : complex(dp), allocatable :: chi2tota(:) ! AR total sum
656 8 : complex(dp), allocatable :: chi2full(:) ! AR+R total sum
657 : ! Addition bands decomposition
658 : character(len=fnlen) :: fnam8
659 : character(len=fnlen) :: fnam9
660 : character(len=fnlen) :: fnam10
661 : character(len=fnlen) :: fnam11
662 : character(len=fnlen) :: fnam12
663 : integer :: fout8
664 : integer :: fout9
665 : integer :: fout10
666 : integer :: fout11
667 : integer :: fout12
668 : integer :: iw_tgt
669 : real(dp) :: iw_real
670 : real(dp) :: ev2ha
671 8 : complex(dp), allocatable :: inter2w_bands(:,:,:), inter2w_bands_ik(:,:,:)
672 8 : complex(dp), allocatable :: inter1w_bands(:,:,:), inter1w_bands_ik(:,:,:)
673 8 : complex(dp), allocatable :: intra2w_bands(:,:,:), intra2w_bands_ik(:,:,:)
674 8 : complex(dp), allocatable :: intra1w_bands(:,:,:), intra1w_bands_ik(:,:,:)
675 8 : complex(dp), allocatable :: intra1wS_bands(:,:,:), intra1wS_bands_ik(:,:,:)
676 : ! Addition 2bands interactions decomposition
677 : character(len=fnlen) :: fnam13
678 : character(len=fnlen) :: fnam14
679 : integer :: fout13
680 : integer :: fout14
681 8 : complex(dp), allocatable :: intra2w_2bands(:,:), intra2w_2bands_ik(:,:)
682 8 : complex(dp), allocatable :: intra1wS_2bands(:,:), intra1wS_2bands_ik(:,:)
683 : ! *********************************************************************
684 :
685 : !DEBUG
686 : !write(std_out,*)' nlinopt : enter '
687 : !write(std_out,*)' nlinopt : tol=',tol
688 : !call flush_unit(std_out)
689 : !ENDDEBUG
690 :
691 8 : my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
692 8 : mband = ks_ebands%mband
693 :
694 : !calculate the constant
695 8 : zi=(0._dp,1._dp)
696 8 : idel=zi*brod
697 : !const_au=-1._dp/(cryst%ucvol*dble(cryst%nsym)) ! VT: 1 instead of 2 bc 2 is unexplained
698 8 : const_au=-2._dp/(cryst%ucvol*dble(cryst%nsym))
699 8 : au2esu=5.8300348177d-8 ! REPLACE WITH DATA FROM DEFS_BASIS
700 8 : const_esu=const_au*au2esu
701 : if (do_decompo) then
702 : ev2ha = 1._dp/Ha_eV
703 : end if ! do_decompo
704 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
705 : !5.8300348177d-8 : au2esu : bohr*c*10^4/4pi*2*ry2ev
706 : !bohr: 5.2917ifc nlinopt.f907E-11
707 : !c: 2.99792458 velocity of light
708 : !au2esu=(5.29177E-11*2.99792458*1.0E4)/Ha_eV
709 : !this const includes (e^3*hbar^3*hbar^3)/(vol*hbar^5*m_e^3)
710 : !mass comes from converting P_mn to r_mn
711 : !hbar^3 comes from converting all frequencies to energies in denominator
712 : !hbar^3 comes from operator for momentum (hbar/i nabla)
713 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
714 : !output file names
715 8 : fnam1=trim(fnam)//'-ChiTotIm.out'
716 8 : fnam2=trim(fnam)//'-ChiTotRe.out'
717 8 : fnam3=trim(fnam)//'-ChiIm.out'
718 8 : fnam4=trim(fnam)//'-ChiRe.out'
719 8 : fnam5=trim(fnam)//'-ChiAbs.out'
720 8 : fnam6=trim(fnam)//'-ChiImDec.out'
721 8 : fnam7=trim(fnam)//'-ChiReDec.out'
722 : ! Addition bands decomposition
723 8 : if (do_decompo) then
724 0 : fnam8=trim(fnam)//'-ChiInter2wBands.out'
725 0 : fnam9=trim(fnam)//'-ChiInter1wBands.out'
726 0 : fnam10=trim(fnam)//'-ChiIntra2wBands.out'
727 0 : fnam11=trim(fnam)//'-ChiIntra1wBands.out'
728 0 : fnam12=trim(fnam)//'-ChiIntra1wSBands.out'
729 : ! Addition 2bands interactions decomposition
730 0 : fnam13=trim(fnam)//'-ChiIntra2w2Bands.out'
731 0 : fnam14=trim(fnam)//'-ChiIntra1wS2Bands.out'
732 : end if ! do_decompo
733 :
734 8 : if(my_rank == master) then
735 : ! If there exists inversion symmetry exit with a message.
736 8 : if (cryst%idx_spatial_inversion() /= 0) then
737 0 : write(std_out,*) '-------------------------------------------'
738 0 : write(std_out,*) ' The crystal has inversion symmetry '
739 0 : write(std_out,*) ' The SHG susceptibility is zero '
740 0 : write(std_out,*) ' Action : set num_nonlin_comp to zero '
741 0 : write(std_out,*) '-------------------------------------------'
742 0 : ABI_ERROR("Aborting now")
743 : end if
744 : ! check polarisation
745 8 : if (v1.le.0.or.v2.le.0.or.v3.le.0.or.v1.gt.3.or.v2.gt.3.or.v3.gt.3) then
746 0 : write(std_out,*) '---------------------------------------------'
747 0 : write(std_out,*) ' Error in nlinopt: '
748 0 : write(std_out,*) ' Incorrect polarisation directions '
749 0 : write(std_out,*) ' 1=x, 2=y and 3=z '
750 0 : write(std_out,*) ' Action : check your input file, '
751 0 : write(std_out,*) ' use only 1, 2 or 3 to define directions '
752 0 : write(std_out,*) '---------------------------------------------'
753 0 : ABI_ERROR("Aborting now")
754 : end if
755 : !number of energy mesh points
756 8 : if (nmesh.le.0) then
757 0 : write(std_out,*) '---------------------------------------------'
758 0 : write(std_out,*) ' Error in nlinopt: '
759 0 : write(std_out,*) ' number of energy mesh points incorrect '
760 0 : write(std_out,*) ' number has to be integer greater than 0 '
761 0 : write(std_out,*) ' nmesh*de = max energy for calculation '
762 0 : write(std_out,*) '---------------------------------------------'
763 0 : ABI_ERROR("Aborting now")
764 : end if
765 : !step in energy
766 8 : if (de.le.zero) then
767 0 : write(std_out,*) '---------------------------------------------'
768 0 : write(std_out,*) ' Error in nlinopt: '
769 0 : write(std_out,*) ' energy step is incorrect '
770 0 : write(std_out,*) ' number has to real greater than 0.0 '
771 0 : write(std_out,*) ' nmesh*de = max energy for calculation '
772 0 : write(std_out,*) '---------------------------------------------'
773 0 : ABI_ERROR("Aborting now")
774 : end if
775 : !broadening
776 8 : if (brod.gt.0.009) then
777 0 : write(std_out,*) '---------------------------------------------'
778 0 : write(std_out,*) ' WARNING : broadening is quite high '
779 0 : write(std_out,*) ' ideally should be less than 0.005 '
780 0 : write(std_out,*) '---------------------------------------------'
781 : else if (brod.gt.0.015) then
782 : write(std_out,*) '----------------------------------------'
783 : write(std_out,*) ' WARNING : broadening is too high '
784 : write(std_out,*) ' ideally should be less than 0.005 '
785 : write(std_out,*) '----------------------------------------'
786 : end if
787 : !tolerance
788 8 : if (tol.gt.0.006) then
789 0 : write(std_out,*) '----------------------------------------'
790 0 : write(std_out,*) ' WARNING : tolerance is too high '
791 0 : write(std_out,*) ' ideally should be less than 0.004 '
792 0 : write(std_out,*) '----------------------------------------'
793 : end if
794 : end if
795 :
796 : !allocate local arrays
797 56 : ABI_MALLOC(px, (mband, mband, 3, 3, 3))
798 24 : ABI_MALLOC(py, (mband, mband, 3, 3, 3))
799 24 : ABI_MALLOC(pz, (mband, mband, 3, 3, 3))
800 24 : ABI_MALLOC(inter2w, (nmesh))
801 16 : ABI_MALLOC(inter1w, (nmesh))
802 16 : ABI_MALLOC(intra2w, (nmesh))
803 16 : ABI_MALLOC(intra1w, (nmesh))
804 16 : ABI_MALLOC(intra1wS, (nmesh))
805 32 : ABI_MALLOC(delta, (mband, mband, 3))
806 : ! Addition antiresonant (AR)
807 8 : if (.not.do_antiresonant) then
808 8 : ABI_MALLOC(inter2wa, (nmesh))
809 8 : ABI_MALLOC(inter1wa, (nmesh))
810 8 : ABI_MALLOC(intra2wa, (nmesh))
811 8 : ABI_MALLOC(intra1wa, (nmesh))
812 8 : ABI_MALLOC(intra1wSa, (nmesh))
813 : end if ! if .not.do_antiresonant
814 : ! Addition bands decomposition
815 8 : if (do_decompo) then
816 0 : if (contrib_decompo==0 .or. contrib_decompo==12) then
817 0 : ABI_MALLOC(inter2w_bands, (nband_sum,nband_sum,nband_sum)) ! (l,m,n)
818 0 : ABI_MALLOC(inter2w_bands_ik, (nband_sum,nband_sum,nband_sum)) ! (l,m,n)
819 : end if
820 0 : if (contrib_decompo==0 .or. contrib_decompo==11) then
821 0 : ABI_MALLOC(inter1w_bands, (nband_sum,nband_sum,nband_sum)) ! (l,m,n)
822 0 : ABI_MALLOC(inter1w_bands_ik, (nband_sum,nband_sum,nband_sum)) ! (l,m,n)
823 : end if
824 0 : if (contrib_decompo==0 .or. contrib_decompo==22) then
825 0 : ABI_MALLOC(intra2w_bands, (nband_sum,nband_sum,nband_sum)) ! (l,m,n)
826 0 : ABI_MALLOC(intra2w_bands_ik, (nband_sum,nband_sum,nband_sum)) ! (l,m,n)
827 : end if
828 0 : if (contrib_decompo==0 .or. contrib_decompo==21) then
829 0 : ABI_MALLOC(intra1w_bands, (nband_sum,nband_sum,nband_sum)) ! (l,m,n)
830 0 : ABI_MALLOC(intra1w_bands_ik, (nband_sum,nband_sum,nband_sum)) ! (l,m,n)
831 : end if
832 0 : if (contrib_decompo==0 .or. contrib_decompo==1) then
833 0 : ABI_MALLOC(intra1wS_bands, (nband_sum,nband_sum,nband_sum)) ! (l,m,n)
834 0 : ABI_MALLOC(intra1wS_bands_ik, (nband_sum,nband_sum,nband_sum)) ! (l,m,n)
835 : end if
836 0 : if (contrib_decompo==0 .or. contrib_decompo==2) then
837 : ! Addition 2bands interactions decomposition
838 0 : ABI_MALLOC(intra2w_2bands, (nband_sum,nband_sum)) ! (m,n)
839 0 : ABI_MALLOC(intra2w_2bands_ik, (nband_sum,nband_sum)) ! (m,n)
840 0 : ABI_MALLOC(intra1wS_2bands, (nband_sum,nband_sum)) ! (m,n)
841 0 : ABI_MALLOC(intra1wS_2bands_ik, (nband_sum,nband_sum)) ! (m,n)
842 : end if
843 : end if ! do_decompo
844 :
845 : !generate the symmetrizing tensor
846 8 : sym = zero
847 200 : do isym=1,cryst%nsym
848 2496 : s(:,:)=cryst%symrel_cart(:,:,isym)
849 776 : do i=1,3
850 2496 : do j=1,3
851 7488 : do k=1,3
852 6912 : sym(i,j,k)=sym(i,j,k)+(s(i,v1)*s(j,v2)*s(k,v3))
853 : end do
854 : end do
855 : end do
856 : end do
857 : ! Disable symmetries for now
858 : !sym(:,:,:) = zero
859 : !sym(v1,v2,v3) = nsym
860 :
861 : ! Split work
862 8 : call xmpi_split_work(ks_ebands%nkpt, comm, my_k1, my_k2)
863 :
864 : ! initialise
865 2808 : inter2w(:)=zero
866 2808 : inter1w(:)=zero
867 2808 : intra2w(:)=zero
868 2808 : intra1w(:)=zero
869 2808 : intra1wS(:)=zero
870 22892 : delta(:,:,:)=zero
871 : ! Addition antiresonant (AR)
872 8 : if (.not.do_antiresonant) then
873 604 : inter2wa(:)=zero
874 604 : inter1wa(:)=zero
875 604 : intra2wa(:)=zero
876 604 : intra1wa(:)=zero
877 604 : intra1wSa(:)=zero
878 : end if ! if .not.do_antiresonant
879 : ! Addition bands decomposition
880 8 : if (do_decompo) then
881 0 : if (contrib_decompo==0 .or. contrib_decompo==12) then
882 0 : inter2w_bands(:,:,:) = zero
883 0 : inter2w_bands_ik(:,:,:) = zero
884 : end if
885 0 : if (contrib_decompo==0 .or. contrib_decompo==11) then
886 0 : inter1w_bands(:,:,:) = zero
887 0 : inter1w_bands_ik(:,:,:) = zero
888 : end if
889 0 : if (contrib_decompo==0 .or. contrib_decompo==22) then
890 0 : intra2w_bands(:,:,:) = zero
891 0 : intra2w_bands_ik(:,:,:) = zero
892 : end if
893 0 : if (contrib_decompo==0 .or. contrib_decompo==21) then
894 0 : intra1w_bands(:,:,:) = zero
895 0 : intra1w_bands_ik(:,:,:) = zero
896 : end if
897 0 : if (contrib_decompo==0 .or. contrib_decompo==1) then
898 0 : intra1wS_bands(:,:,:) = zero
899 0 : intra1wS_bands_ik(:,:,:) = zero
900 : end if
901 0 : if (contrib_decompo==0 .or. contrib_decompo==2) then
902 : ! Addition 2bands interactions decomposition
903 0 : intra2w_2bands(:,:) = zero
904 0 : intra2w_2bands_ik(:,:) = zero
905 0 : intra1wS_2bands(:,:) = zero
906 0 : intra1wS_2bands_ik(:,:) = zero
907 : end if
908 : end if ! do_decompo
909 :
910 8 : my_emin=HUGE(zero)
911 8 : my_emax=-HUGE(zero)
912 :
913 : ! loop over kpts
914 680 : do ik=my_k1,my_k2
915 672 : write(std_out,*) "P-",my_rank,": ",ik,'of',ks_ebands%nkpt
916 : ! loop over spins
917 1352 : do isp=1,ks_ebands%nsppol
918 : ! loop over states
919 15072 : do ist1=1,nband_sum
920 : ! Addition antiresonant (AR) commented
921 : !e1 = ks_ebands%eig(ist1,ik,isp)
922 : !if (e1.lt.ks_ebands%fermie) then ! ist1 is a valence state
923 351072 : do ist2=1,nband_sum
924 : ! Addition antiresonant (AR) commented
925 : !e2 = ks_ebands%eig(ist2,ik,isp)
926 : !if (e2.gt.ks_ebands%fermie) then ! ist2 is a conduction state
927 : ! symmetrize the momentum matrix elements
928 1358400 : do lx=1,3
929 4368000 : do ly=1,3
930 13104000 : do lz=1,3
931 9072000 : f1=sym(lx,ly,lz)+sym(lx,lz,ly)
932 9072000 : f2=sym(ly,lx,lz)+sym(ly,lz,lx)
933 9072000 : f3=sym(lz,lx,ly)+sym(lz,ly,lx)
934 9072000 : px(ist1,ist2,lx,ly,lz)=f1*pmat(ist1,ist2,ik,lx,isp)
935 9072000 : py(ist2,ist1,lx,ly,lz)=f2*pmat(ist2,ist1,ik,lx,isp)
936 12096000 : pz(ist2,ist1,lx,ly,lz)=f3*pmat(ist2,ist1,ik,lx,isp)
937 : end do !lz
938 : end do !ly
939 : end do ! lx end loop over states
940 : !end if ! e2.gt.fermie
941 : end do ! ist2
942 : !end if ! e1.lt.fermie
943 : end do ! ist1
944 :
945 : ! calculate the energy window and \Delta_nm
946 15072 : do ist1=1,nband_sum
947 14400 : my_emin=min(my_emin, ks_ebands%eig(ist1,ik,isp))
948 14400 : my_emax=max(my_emax, ks_ebands%eig(ist1,ik,isp))
949 351072 : do ist2=1,nband_sum
950 1358400 : delta(ist1,ist2,1:3)=pmat(ist1,ist1,ik,1:3,isp)-pmat(ist2,ist2,ik,1:3,isp)
951 : end do
952 : end do
953 : ! initialise the factors
954 : ! factors are named according to the Ref. article 2.
955 672 : b111=zero
956 672 : b121=zero
957 672 : b131=zero
958 672 : b112=zero
959 672 : b122=zero
960 672 : b132=zero
961 672 : b113=zero
962 672 : b123=zero
963 672 : b133=zero
964 672 : b211=zero
965 672 : b221=zero
966 672 : b212=zero
967 672 : b222=zero
968 672 : b213=zero
969 672 : b223=zero
970 672 : b231=zero
971 672 : b241=zero
972 672 : b242=zero
973 672 : b243=zero
974 672 : b311=zero
975 672 : b312=zero
976 672 : b313=zero
977 672 : b331=zero
978 : ! Addition AR
979 : ! factors are named as above but b->a for AR
980 672 : if (.not.do_antiresonant) then
981 128 : a111=zero
982 128 : a112=zero
983 128 : a113=zero
984 128 : a121=zero
985 128 : a131=zero
986 128 : a122=zero
987 128 : a132=zero
988 128 : a123=zero
989 128 : a133=zero
990 128 : a241=zero
991 128 : a242=zero
992 128 : a243=zero
993 128 : a231=zero
994 128 : a211=zero
995 128 : a221=zero
996 128 : a212=zero
997 128 : a222=zero
998 128 : a213=zero
999 128 : a223=zero
1000 128 : a311=zero
1001 128 : a312=zero
1002 128 : a313=zero
1003 128 : a331=zero
1004 : end if ! if .not.do_antiresonant
1005 : ! start the calculation
1006 15744 : do istn=1,nband_sum
1007 14400 : en=ks_ebands%eig(istn,ik,isp)
1008 15072 : if (en.lt.ks_ebands%fermie) then ! istn is a valence state
1009 60288 : do istm=1,nband_sum
1010 57600 : em=ks_ebands%eig(istm,ik,isp)
1011 60288 : if (em.gt.ks_ebands%fermie) then ! istm is a conduction state
1012 46848 : em = em + sc ! Should add the scissor to conduction energies
1013 46848 : wmn=em-en
1014 46848 : wnm=-wmn
1015 : ! calculate the matrix elements for two band intraband term
1016 46848 : mat2w_tra=zero
1017 46848 : mat1w3_tra=zero
1018 : ! Addition AR
1019 46848 : if (.not.do_antiresonant) then
1020 8192 : mat2wa_tra = zero
1021 8192 : mat1w3a_tra = zero
1022 : end if ! if .not.do_antiresonant
1023 187392 : do lx=1,3
1024 609024 : do ly=1,3
1025 1827072 : do lz=1,3
1026 : ! See Sharma03, A10, last term of first line, corrected ! Should be Delta^b_mn r c_mn, see A2.
1027 : mat2w_tra=mat2w_tra+px(istn,istm,lx,ly,lz)*pmat(istm,istn,ik,lz,isp) &
1028 1264896 : *delta(istm,istn,ly)
1029 : ! See Sharma03, A11, last term.
1030 : mat1w3_tra=mat1w3_tra+px(istn,istm,lx,ly,lz)*pmat(istm,istn,ik,ly,isp) &
1031 1264896 : *delta(istm,istn,lz)
1032 : ! NOTE:: lx to ly m to n in pmat matrices respectively
1033 : ! Changes are made so that this (b3) term is according to paper
1034 : ! [[cite:Sipe1993]] (Ref. 4) rather than [[cite:Hughes1996]] (Ref 2) in which this term is incorrect
1035 :
1036 : ! Addition AR
1037 1686528 : if (.not.do_antiresonant) then
1038 : mat2wa_tra=mat2wa_tra+px(istm,istn,lx,ly,lz)*pmat(istn,istm,ik,lz,isp) &
1039 221184 : *delta(istn,istm,ly)
1040 :
1041 : mat1w3a_tra=mat1w3a_tra+px(istm,istn,lx,ly,lz)*pmat(istn,istm,ik,ly,isp) &
1042 221184 : *delta(istn,istm,lz)
1043 : end if ! if .not.do_antiresonant
1044 :
1045 : end do
1046 : end do
1047 : end do
1048 46848 : b331=mat1w3_tra/wnm
1049 46848 : b231=8._dp*mat2w_tra/wmn
1050 : ! Addition AR
1051 46848 : if (.not.do_antiresonant) then
1052 8192 : a331=mat1w3a_tra/wmn
1053 8192 : a231=8._dp*mat2wa_tra/wnm
1054 : end if ! if .not.do_antiresonant
1055 : ! Addition 2bands interactions decomposition
1056 46848 : if (do_decompo) then
1057 0 : if (contrib_decompo==0 .or. contrib_decompo==2) then
1058 0 : intra2w_2bands_ik(istm, istn) = b231
1059 0 : intra1wS_2bands_ik(istm, istn) = b331
1060 : end if
1061 : end if ! do_decompo
1062 :
1063 46848 : b11=zero
1064 46848 : b12_13=zero
1065 46848 : b24=zero
1066 46848 : b31_32=zero
1067 46848 : b21_22=zero
1068 : ! Addition AR
1069 46848 : if (.not.do_antiresonant) then
1070 8192 : a11 = zero
1071 8192 : a12_13 = zero
1072 8192 : a24 = zero
1073 8192 : a21_22 = zero
1074 8192 : a31_32 = zero
1075 : end if ! if .not.do_antiresonant
1076 : ! istl < istn
1077 117120 : do istl=1,istn-1 ! istl is a valence state below istn
1078 70272 : el=ks_ebands%eig(istl,ik,isp)
1079 70272 : wln=el-en ! do not add sc to the valence bands!
1080 70272 : wml=em-el
1081 70272 : wnl=-wln
1082 70272 : wlm=-wml
1083 : ! calculate the matrix elements for three band terms
1084 70272 : mat2w=zero
1085 70272 : mat1w1=zero
1086 70272 : mat1w2=zero
1087 : ! Addition AR
1088 70272 : if (.not.do_antiresonant) then
1089 12288 : mat2wa = zero
1090 12288 : mat1w1a = zero
1091 12288 : mat1w2a = zero
1092 : end if ! if .not.do_antiresonant
1093 281088 : do lx=1,3
1094 913536 : do ly=1,3
1095 2740608 : do lz=1,3
1096 :
1097 : mat2w=mat2w+(px(istn,istm,lx,ly,lz)*pmat(istm,istl,ik,ly,isp) &
1098 1897344 : *pmat(istl,istn,ik,lz,isp))
1099 :
1100 : mat1w1=mat1w1+(py(istm,istn,lx,ly,lz)*pmat(istl,istm,ik,lz,isp) &
1101 1897344 : *pmat(istn,istl,ik,ly,isp))
1102 :
1103 : mat1w2=mat1w2+(pz(istm,istn,lx,ly,lz)*pmat(istl,istm,ik,lz,isp) &
1104 1897344 : *pmat(istn,istl,ik,ly,isp))
1105 :
1106 : ! Addition AR
1107 2529792 : if (.not.do_antiresonant) then
1108 : mat2wa=mat2wa+(px(istm,istn,lx,ly,lz)*pmat(istn,istl,ik,ly,isp) &
1109 331776 : *pmat(istl,istm,ik,lz,isp))
1110 :
1111 : mat1w1a=mat1w1a+(py(istn,istm,lx,ly,lz)*pmat(istl,istn,ik,lz,isp) &
1112 331776 : *pmat(istm,istl,ik,ly,isp))
1113 :
1114 : mat1w2a=mat1w2a+(pz(istn,istm,lx,ly,lz)*pmat(istl,istn,ik,lz,isp) &
1115 331776 : *pmat(istm,istl,ik,ly,isp))
1116 : end if ! if .not.do_antiresonant
1117 :
1118 : end do
1119 : end do
1120 : end do
1121 70272 : b111=mat2w*(1._dp/(wln+wlm))*(1._dp/wlm)
1122 70272 : b121=mat1w1*(1._dp/(wnm+wlm))*(1._dp/wlm)
1123 70272 : b131=mat1w2*(1._dp/wlm)
1124 70272 : b221=zero
1125 70272 : b211=mat1w1/wml
1126 70272 : b241=-mat2w/wml
1127 70272 : b311=mat1w2/wlm
1128 : ! Addition AR
1129 70272 : if (.not.do_antiresonant) then
1130 12288 : a111=mat2wa*(1._dp/(wln+wlm))*(1._dp/wml)
1131 12288 : a121=mat1w1a*(1._dp/wml) !*(1._dp/(wmn+wln))
1132 12288 : a131=mat1w2a*(1._dp/(wlm-wmn))*(1._dp/wml)
1133 12288 : a241=-mat2wa/wlm
1134 12288 : a221=mat1w2a/wlm
1135 12288 : a211=zero
1136 12288 : a311=mat1w1a/wlm
1137 : end if ! if .not.do_antiresonant
1138 70272 : if (abs(wln).gt.tol) then
1139 68880 : b111=b111/wln
1140 68880 : b121=b121/wln
1141 68880 : b131=b131/wln
1142 68880 : b221=mat1w2/wln
1143 68880 : b241=b241+(mat2w/wln)
1144 68880 : b311=b311+(mat1w1/wln)
1145 : ! Addition AR
1146 68880 : if (.not.do_antiresonant) then
1147 11776 : a111=a111/wnl
1148 11776 : a121=a121/wnl
1149 11776 : a131=a131/wnl
1150 11776 : a241=a241+(mat2wa/wnl)
1151 11776 : a211=mat1w1a/wnl
1152 11776 : a311=a311+(mat1w2a/wln)
1153 : end if ! if .not.do_antiresonant
1154 : else
1155 1392 : b111=zero
1156 1392 : b121=zero
1157 1392 : b131=zero
1158 1392 : b221=zero
1159 : ! Addition AR
1160 1392 : if (.not.do_antiresonant) then
1161 512 : a111=zero
1162 512 : a121=zero
1163 512 : a131=zero
1164 512 : a211=zero
1165 : end if ! if .not.do_antiresonant
1166 : end if
1167 :
1168 70272 : t1=wln-wnm
1169 70272 : if (abs(t1).gt.tol) then
1170 70176 : b131=b131/t1
1171 : ! Addition AR
1172 70176 : if (.not.do_antiresonant) then
1173 12288 : a121=a121/t1
1174 : end if ! if .not.do_antiresonant
1175 : else
1176 96 : b131=zero
1177 : ! Addition AR
1178 96 : if (.not.do_antiresonant) then
1179 0 : a121=zero
1180 : end if ! if .not.do_antiresonant
1181 : end if
1182 70272 : b11=b11-2._dp*b111
1183 70272 : b12_13=b12_13+b121+b131
1184 70272 : b21_22=b21_22-b211+b221
1185 70272 : b24=b24+2._dp*b241
1186 70272 : b31_32=b31_32+b311
1187 : ! Addition AR
1188 70272 : if (.not.do_antiresonant) then
1189 12288 : a11 = a11+2._dp*a111
1190 12288 : a12_13 = a12_13-a121-a131
1191 12288 : a24 = a24+2._dp*a241
1192 12288 : a21_22 = a21_22+a211-a221
1193 12288 : a31_32 = a31_32-a311
1194 : end if ! if .not.do_antiresonant
1195 : ! Addition bands decomposition
1196 117120 : if (do_decompo) then
1197 0 : if (contrib_decompo==0 .or. contrib_decompo==12) then
1198 0 : inter2w_bands_ik(istl, istm, istn) = -2._dp*b111
1199 : end if
1200 0 : if (contrib_decompo==0 .or. contrib_decompo==11) then
1201 0 : inter1w_bands_ik(istl, istm, istn) = b121+b131
1202 : end if
1203 0 : if (contrib_decompo==0 .or. contrib_decompo==22) then
1204 0 : intra2w_bands_ik(istl, istm, istn) = 2._dp*b241
1205 : end if
1206 0 : if (contrib_decompo==0 .or. contrib_decompo==21) then
1207 0 : intra1w_bands_ik(istl, istm, istn) = -b211+b221
1208 : end if
1209 0 : if (contrib_decompo==0 .or. contrib_decompo==1) then
1210 0 : intra1wS_bands_ik(istl, istm, istn) = b311
1211 : end if
1212 : end if ! do_decompo
1213 : end do ! istl
1214 :
1215 : ! istn < istl < istm
1216 556800 : do istl=istn+1,istm-1
1217 509952 : el=ks_ebands%eig(istl,ik,isp)
1218 : ! calculate the matrix elements for three band terms
1219 509952 : mat2w=zero
1220 509952 : mat1w1=zero
1221 509952 : mat1w2=zero
1222 : ! Addition AR
1223 509952 : if (.not.do_antiresonant) then
1224 73728 : mat2wa = zero
1225 73728 : mat1w1a = zero
1226 73728 : mat1w2a = zero
1227 : end if ! if .not.do_antiresonant
1228 2039808 : do lx=1,3
1229 6629376 : do ly=1,3
1230 19888128 : do lz=1,3
1231 :
1232 : mat2w=mat2w+(px(istn,istm,lx,ly,lz)*pmat(istm,istl,ik,ly,isp) &
1233 13768704 : *pmat(istl,istn,ik,lz,isp))
1234 :
1235 : mat1w1=mat1w1+(py(istm,istn,lx,ly,lz)*pmat(istl,istm,ik,lz,isp) &
1236 13768704 : *pmat(istn,istl,ik,ly,isp))
1237 :
1238 : mat1w2=mat1w2+(pz(istm,istn,lx,ly,lz)*pmat(istl,istm,ik,lz,isp) &
1239 13768704 : *pmat(istn,istl,ik,ly,isp))
1240 :
1241 : ! Addition AR
1242 18358272 : if (.not.do_antiresonant) then
1243 : mat2wa=mat2wa+(px(istm,istn,lx,ly,lz)*pmat(istn,istl,ik,ly,isp) &
1244 1990656 : *pmat(istl,istm,ik,lz,isp))
1245 :
1246 : mat1w1a=mat1w1a+(py(istn,istm,lx,ly,lz)*pmat(istl,istn,ik,lz,isp) &
1247 1990656 : *pmat(istm,istl,ik,ly,isp))
1248 :
1249 : mat1w2a=mat1w2a+(pz(istn,istm,lx,ly,lz)*pmat(istl,istn,ik,lz,isp) &
1250 1990656 : *pmat(istm,istl,ik,ly,isp))
1251 : end if ! if .not.do_antiresonant
1252 :
1253 : end do
1254 : end do
1255 : end do
1256 509952 : if (el.lt.ks_ebands%fermie) then
1257 70272 : wln=el-en
1258 70272 : wnl=-wln
1259 70272 : wml=em-el
1260 70272 : wlm=-wml
1261 : else
1262 439680 : el=el+sc
1263 439680 : wln=el-en
1264 439680 : wnl=-wln
1265 439680 : wml=em-el
1266 439680 : wlm=-wml
1267 : end if
1268 : !
1269 509952 : b112=zero
1270 509952 : b122=mat1w1*(1._dp/(wnm+wlm))
1271 509952 : b132=mat1w2*(1._dp/(wnm+wnl))
1272 509952 : b242=zero
1273 509952 : b222=zero
1274 509952 : b212=zero
1275 509952 : b312=zero
1276 : ! Addition AR
1277 509952 : if (.not.do_antiresonant) then
1278 73728 : a112=zero
1279 73728 : a122=mat1w1a*(1._dp/(wmn+wln))
1280 73728 : a132=mat1w2a*(1._dp/(wmn+wml))
1281 73728 : a242=zero
1282 73728 : a212=zero
1283 73728 : a222=zero
1284 73728 : a312=zero
1285 : end if ! .not.do_antiresonant
1286 509952 : if (abs(wnl).gt.tol) then
1287 508560 : b112=mat2w/wln
1288 508560 : b122=b122/wnl
1289 508560 : b132=b132/wnl
1290 508560 : b242=mat2w/wln
1291 508560 : b222=mat1w2/wln
1292 508560 : b312=mat1w1/wln
1293 : ! Addition AR
1294 508560 : if (.not.do_antiresonant) then
1295 73216 : a112=mat2wa/wnl
1296 73216 : a122=a122/wln
1297 73216 : a132=a132/wln
1298 73216 : a242=mat2wa/wnl
1299 73216 : a212=mat1w1a/wnl
1300 73216 : a312=mat1w2a/wln
1301 : end if ! .not.do_antiresonant
1302 : else
1303 1392 : b122=zero
1304 1392 : b132=zero
1305 : ! Addition AR
1306 1392 : if (.not.do_antiresonant) then
1307 512 : a122=zero
1308 512 : a132=zero
1309 : end if ! .not.do_antiresonant
1310 : end if
1311 509952 : if (abs(wlm).gt.tol) then
1312 507872 : b112=b112/wml
1313 507872 : b122=b122/wlm
1314 507872 : b132=b132/wlm
1315 507872 : b242=b242-(mat2w/wml)
1316 507872 : b212=mat1w1/wml
1317 507872 : b312=b312+(mat1w2/wlm)
1318 : ! Addition AR
1319 507872 : if (.not.do_antiresonant) then
1320 73088 : a112=a112/wlm
1321 73088 : a122=a122/wml
1322 73088 : a132=a132/wml
1323 73088 : a242=a242-(mat2wa/wlm)
1324 73088 : a222=mat1w2a/wlm
1325 73088 : a312=a312+(mat1w1a/wlm)
1326 : end if ! .not.do_antiresonant
1327 : else
1328 2080 : b112=zero
1329 2080 : b122=zero
1330 2080 : b132=zero
1331 2080 : b212=zero
1332 : ! Addition AR
1333 2080 : if (.not.do_antiresonant) then
1334 640 : a112=zero
1335 640 : a122=zero
1336 640 : a132=zero
1337 640 : a222=zero
1338 : end if ! .not.do_antiresonant
1339 : end if
1340 509952 : t1=wlm-wnl
1341 509952 : if (abs(t1).gt.tol) then
1342 509184 : b112=b112/t1
1343 : ! Addition AR
1344 509184 : if (.not.do_antiresonant) then
1345 73728 : a112=a112/t1
1346 : end if ! .not.do_antiresonant
1347 : else
1348 768 : b112=zero
1349 : ! Addition AR
1350 768 : if (.not.do_antiresonant) then
1351 0 : a112=zero
1352 : end if ! .not.do_antiresonant
1353 : end if
1354 509952 : b11=b11+2._dp*b112
1355 509952 : b12_13=b12_13-b122+b132
1356 509952 : b24=b24+2._dp*b242
1357 509952 : b21_22=b21_22-b212+b222
1358 509952 : b31_32=b31_32+b312
1359 : ! Addition AR
1360 509952 : if (.not.do_antiresonant) then
1361 73728 : a11 = a11-2._dp*a112
1362 73728 : a12_13 = a12_13+a122-a132
1363 73728 : a24 = a24+2._dp*a242
1364 73728 : a21_22 = a21_22+a212-a222
1365 73728 : a31_32 = a31_32-a312
1366 : end if ! .not.do_antiresonant
1367 : ! Addition bands decomposition
1368 556800 : if (do_decompo) then
1369 0 : if (contrib_decompo==0 .or. contrib_decompo==12) then
1370 0 : inter2w_bands_ik(istl, istm, istn) = 2._dp*b112
1371 : end if
1372 0 : if (contrib_decompo==0 .or. contrib_decompo==11) then
1373 0 : inter1w_bands_ik(istl, istm, istn) = -b122+b132
1374 : end if
1375 0 : if (contrib_decompo==0 .or. contrib_decompo==22) then
1376 0 : intra2w_bands_ik(istl, istm, istn) = 2._dp*b242
1377 : end if
1378 0 : if (contrib_decompo==0 .or. contrib_decompo==21) then
1379 0 : intra1w_bands_ik(istl, istm, istn) = -b212+b222
1380 : end if
1381 0 : if (contrib_decompo==0 .or. contrib_decompo==1) then
1382 0 : intra1wS_bands_ik(istl, istm, istn) = b312
1383 : end if
1384 : end if ! do_decompo
1385 : end do ! istl
1386 :
1387 : ! istl > istm !
1388 486528 : do istl=istm+1,nband_sum
1389 439680 : el=ks_ebands%eig(istl,ik,isp)+sc
1390 439680 : wln=el-en
1391 439680 : wnl=-wln
1392 439680 : wml=em-el
1393 439680 : wlm=-wml
1394 : ! calculate the matrix elements for three band terms
1395 439680 : mat2w=zero
1396 439680 : mat1w1=zero
1397 439680 : mat1w2=zero
1398 : ! Addition AR
1399 439680 : if (.not.do_antiresonant) then
1400 61440 : mat2wa = zero
1401 61440 : mat1w1a = zero
1402 61440 : mat1w2a = zero
1403 : end if ! .not.do_antiresonant
1404 1758720 : do lx=1,3
1405 5715840 : do ly=1,3
1406 17147520 : do lz=1,3
1407 :
1408 : mat2w=mat2w+px(istn,istm,lx,ly,lz)*pmat(istm,istl,ik,ly,isp) &
1409 11871360 : *pmat(istl,istn,ik,lz,isp)
1410 :
1411 : mat1w1=mat1w1+(py(istm,istn,lx,ly,lz)*pmat(istl,istm,ik,lz,isp) &
1412 11871360 : *pmat(istn,istl,ik,ly,isp))
1413 :
1414 : mat1w2=mat1w2+(pz(istm,istn,lx,ly,lz)*pmat(istl,istm,ik,lz,isp) &
1415 11871360 : *pmat(istn,istl,ik,ly,isp))
1416 :
1417 : ! Addition AR
1418 15828480 : if (.not.do_antiresonant) then
1419 : mat2wa=mat2wa+px(istm,istn,lx,ly,lz)*pmat(istn,istl,ik,ly,isp) &
1420 1658880 : *pmat(istl,istm,ik,lz,isp)
1421 :
1422 : mat1w1a=mat1w1a+(py(istn,istm,lx,ly,lz)*pmat(istl,istn,ik,lz,isp) &
1423 1658880 : *pmat(istm,istl,ik,ly,isp))
1424 :
1425 : mat1w2a=mat1w2a+(pz(istn,istm,lx,ly,lz)*pmat(istl,istn,ik,lz,isp) &
1426 1658880 : *pmat(istm,istl,ik,ly,isp))
1427 : end if ! .not.do_antiresonant
1428 :
1429 : end do
1430 : end do
1431 : end do
1432 :
1433 439680 : b113=mat2w*(1._dp/(wnl+wml))*(1._dp/wnl)
1434 439680 : b123=mat1w1*(1._dp/wnl)
1435 439680 : b133=mat1w2*(1._dp/wnl)*(1._dp/(wnl+wnm))
1436 439680 : b243=mat2w/wln
1437 439680 : b223=mat1w2/wln
1438 439680 : b213=zero
1439 439680 : b313=-1._dp*mat1w1/wnl
1440 : ! Addition AR
1441 439680 : if (.not.do_antiresonant) then
1442 61440 : a113=mat2wa*(1._dp/(wnl+wml))*(1._dp/wln)
1443 61440 : a123=mat1w1a*(1._dp/(wmn-wnl))*(1._dp/wln)
1444 61440 : a133=mat1w2a*(1._dp/wln)!*(1._dp/(wml+wmn))
1445 61440 : a243=mat2wa/wnl
1446 61440 : a223=zero
1447 61440 : a213=mat1w1a/wnl
1448 61440 : a313=mat1w2a/wln
1449 : end if ! .not.do_antiresonant
1450 439680 : if (abs(wml).gt.tol) then
1451 437600 : b113=b113/wml
1452 437600 : b123=b123/wml
1453 437600 : b133=b133/wml
1454 437600 : b243=b243-(mat2w/wml)
1455 437600 : b213=mat1w1/wml
1456 437600 : b313=b313+(mat1w2/wlm)
1457 : ! Addition AR
1458 437600 : if (.not.do_antiresonant) then
1459 60800 : a113=a113/wlm
1460 60800 : a123=a123/wlm
1461 60800 : a133=a133/wlm
1462 60800 : a243=a243-(mat2wa/wlm)
1463 60800 : a223=mat1w2a/wlm
1464 60800 : a313=a313-(mat1w1a/wml)
1465 : end if ! .not.do_antiresonant
1466 : else
1467 2080 : b113=zero
1468 2080 : b123=zero
1469 2080 : b133=zero
1470 : ! Addition AR
1471 2080 : if (.not.do_antiresonant) then
1472 640 : a113=zero
1473 640 : a123=zero
1474 640 : a133=zero
1475 : end if ! .not.do_antiresonant
1476 : end if
1477 :
1478 439680 : t1=wnm-wml
1479 439680 : if (abs(t1).gt.tol) then
1480 439008 : b123=b123/t1
1481 : else
1482 : b123=zero
1483 : end if
1484 :
1485 439680 : if (.not.do_antiresonant) then
1486 61440 : t2=(wml+wmn)
1487 61440 : if (abs(t2).gt.tol) then
1488 61440 : a133=a133/t2
1489 : else
1490 : a133=zero
1491 : end if
1492 : end if ! .not.do_antiresonant
1493 439680 : b11=b11+2._dp*b113
1494 439680 : b12_13=b12_13+b123-b133
1495 439680 : b21_22=b21_22-b213+b223
1496 439680 : b24=b24+2._dp*b243
1497 439680 : b31_32=b31_32+b313
1498 : ! Addition AR
1499 439680 : if (.not.do_antiresonant) then
1500 61440 : a11 = a11-2._dp*a113
1501 61440 : a12_13 = a12_13-a123+a133
1502 61440 : a24 = a24+2._dp*a243
1503 61440 : a21_22 = a21_22+a213-a223
1504 61440 : a31_32 = a31_32-a313
1505 : end if ! .not.do_antiresonant
1506 : ! Addition bands decomposition
1507 486528 : if (do_decompo) then
1508 0 : if (contrib_decompo==0 .or. contrib_decompo==12) then
1509 0 : inter2w_bands_ik(istl, istm, istn) = 2._dp*b113
1510 : end if
1511 0 : if (contrib_decompo==0 .or. contrib_decompo==11) then
1512 0 : inter1w_bands_ik(istl, istm, istn) = b123-b133
1513 : end if
1514 0 : if (contrib_decompo==0 .or. contrib_decompo==22) then
1515 0 : intra2w_bands_ik(istl, istm, istn) = 2._dp*b243
1516 : end if
1517 0 : if (contrib_decompo==0 .or. contrib_decompo==21) then
1518 0 : intra1w_bands_ik(istl, istm, istn) = -b213+b223
1519 : end if
1520 0 : if (contrib_decompo==0 .or. contrib_decompo==1) then
1521 0 : intra1wS_bands_ik(istl, istm, istn) = b313
1522 : end if
1523 : end if ! do_decompo
1524 : end do ! istl
1525 :
1526 46848 : b11 = b11*zi*(1._dp/wnm)*const_esu
1527 46848 : b12_13 = b12_13*zi*(1._dp/wnm)*const_esu
1528 46848 : b24 = (b24+b231)*zi*(1._dp/(wnm**3))*const_esu
1529 46848 : b21_22 = (b21_22)*zi*(1._dp/(wnm**3))*const_esu
1530 46848 : b31_32 = (b31_32-b331)*zi*(1._dp/(wmn**3))*const_esu*0.5_dp
1531 : ! Addition AR
1532 46848 : if (.not.do_antiresonant) then
1533 8192 : a11 = a11*zi*(1._dp/wmn)*const_esu
1534 8192 : a12_13 = a12_13*zi*(1._dp/wmn)*const_esu
1535 8192 : a24 = (a24-a231)*zi*(1._dp/(wmn**3))*const_esu
1536 8192 : a21_22 = (a21_22)*zi*(1._dp/(wmn**3))*const_esu
1537 8192 : a31_32 = (a31_32+a331)*zi*(1._dp/(wnm**3))*const_esu*0.5_dp
1538 : end if ! .not.do_antiresonant
1539 : ! Addition bands decomposition
1540 46848 : if (do_decompo) then
1541 0 : if (contrib_decompo==0 .or. contrib_decompo==12) then
1542 0 : inter2w_bands_ik(:, istm, istn) = inter2w_bands_ik(:, istm, istn)*zi*(1._dp/wnm)*const_esu
1543 : end if
1544 0 : if (contrib_decompo==0 .or. contrib_decompo==11) then
1545 0 : inter1w_bands_ik(:, istm, istn) = inter1w_bands_ik(:, istm, istn)*zi*(1._dp/wnm)*const_esu
1546 : end if
1547 0 : if (contrib_decompo==0 .or. contrib_decompo==22) then
1548 0 : intra2w_bands_ik(:, istm, istn) = intra2w_bands_ik(:, istm, istn)*zi*(1._dp/(wnm**3))*const_esu
1549 : end if
1550 0 : if (contrib_decompo==0 .or. contrib_decompo==21) then
1551 0 : intra1w_bands_ik(:, istm, istn) = intra1w_bands_ik(:, istm, istn)*zi*(1._dp/(wnm**3))*const_esu
1552 : end if
1553 0 : if (contrib_decompo==0 .or. contrib_decompo==1) then
1554 0 : intra1wS_bands_ik(:, istm, istn) = intra1wS_bands_ik(:, istm, istn)*zi*(1._dp/(wmn**3))*const_esu*0.5_dp
1555 : end if
1556 0 : if (contrib_decompo==0 .or. contrib_decompo==2) then
1557 : ! Addition 2bands interactions decomposition
1558 0 : intra2w_2bands_ik(istm, istn) = intra2w_2bands_ik(istm, istn)*zi*(1._dp/(wnm**3))*const_esu
1559 0 : intra1wS_2bands_ik(istm, istn) = -intra1wS_2bands_ik(istm, istn)*zi*(1._dp/(wmn**3))*const_esu*0.5_dp
1560 : end if
1561 : end if ! do_decompo
1562 :
1563 : ! calculate over the desired energy mesh and sum over k-points
1564 34632448 : do iw=1,nmesh
1565 34585600 : w=(iw-1)*de+idel
1566 34585600 : inter2w(iw)=inter2w(iw)+(ks_ebands%wtk(ik)*(b11/(wmn-2._dp*w))) ! Inter(2w) from chi
1567 34585600 : inter1w(iw)=inter1w(iw)+(ks_ebands%wtk(ik)*(b12_13/(wmn-w))) ! Inter(1w) from chi
1568 34585600 : intra2w(iw)=intra2w(iw)+(ks_ebands%wtk(ik)*(b24/(wmn-2._dp*w))) ! Intra(2w) from eta
1569 34585600 : intra1w(iw)=intra1w(iw)+(ks_ebands%wtk(ik)*((b21_22)/(wmn-w))) ! Intra(1w) from eta
1570 34585600 : intra1wS(iw)=intra1wS(iw)+(ks_ebands%wtk(ik)*((b31_32)/(wmn-w))) ! Intra(1w) from sigma
1571 : ! Addition AR
1572 34585600 : if (.not.do_antiresonant) then
1573 1228800 : inter2wa(iw)=inter2wa(iw)+(ks_ebands%wtk(ik)*(a11/(wnm-2._dp*dble(w)))) ! Inter(2w) from chi AR
1574 1228800 : inter1wa(iw)=inter1wa(iw)+(ks_ebands%wtk(ik)*(a12_13/(wnm-dble(w)))) ! Inter(1w) from chi AR
1575 1228800 : intra2wa(iw)=intra2wa(iw)+(ks_ebands%wtk(ik)*(a24/(wnm-2._dp*dble(w)))) ! Intra(2w) from eta AR
1576 1228800 : intra1wa(iw)=intra1wa(iw)+(ks_ebands%wtk(ik)*((a21_22)/(wnm-dble(w)))) ! Intra(1w) from eta AR
1577 1228800 : intra1wSa(iw)=intra1wSa(iw)+(ks_ebands%wtk(ik)*((a31_32)/(wnm-dble(w)))) ! Intra(1w) from sigma AR
1578 : end if ! .not.do_antiresonant
1579 : ! Addition bands decomposition
1580 34632448 : if (do_decompo) then
1581 0 : iw_real = (w_decompo*ev2ha)/de + 1
1582 0 : iw_tgt = idnint(iw_real)
1583 0 : if (iw==iw_tgt) then ! need to restrict to one frequency due to memory issues
1584 0 : if (contrib_decompo==0 .or. contrib_decompo==12) then
1585 0 : inter2w_bands(:, istm, istn) = inter2w_bands(:,istm,istn)+(ks_ebands%wtk(ik)*(inter2w_bands_ik(:, istm, istn)/(wmn-2._dp*w)))
1586 : end if
1587 0 : if (contrib_decompo==0 .or. contrib_decompo==11) then
1588 0 : inter1w_bands(:, istm, istn) = inter1w_bands(:,istm,istn)+(ks_ebands%wtk(ik)*(inter1w_bands_ik(:, istm, istn)/(wmn-w)))
1589 : end if
1590 0 : if (contrib_decompo==0 .or. contrib_decompo==22) then
1591 0 : intra2w_bands(:, istm, istn) = intra2w_bands(:,istm,istn)+(ks_ebands%wtk(ik)*(intra2w_bands_ik(:, istm, istn)/(wmn-2._dp*w)))
1592 : end if
1593 0 : if (contrib_decompo==0 .or. contrib_decompo==21) then
1594 0 : intra1w_bands(:, istm, istn) = intra1w_bands(:,istm,istn)+(ks_ebands%wtk(ik)*((intra1w_bands_ik(:, istm, istn))/(wmn-w)))
1595 : end if
1596 0 : if (contrib_decompo==0 .or. contrib_decompo==1) then
1597 0 : intra1wS_bands(:, istm, istn) = intra1wS_bands(:,istm,istn)+(ks_ebands%wtk(ik)*((intra1wS_bands_ik(:, istm, istn))/(wmn-w)))
1598 : end if
1599 0 : if (contrib_decompo==0 .or. contrib_decompo==2) then
1600 : ! Addition 2bands interaction decomposition
1601 0 : intra2w_2bands(istm, istn) = intra2w_2bands(istm,istn)+(ks_ebands%wtk(ik)*(intra2w_2bands_ik(istm, istn)/(wmn-2._dp*w)))
1602 0 : intra1wS_2bands(istm, istn) = intra1wS_2bands(istm,istn)+(ks_ebands%wtk(ik)*((intra1wS_2bands_ik(istm, istn))/(wmn-w)))
1603 : end if
1604 : end if ! Chosen frequency for bands decomposition
1605 : end if ! do_decompo
1606 : end do ! iw
1607 : end if ! istm is CB
1608 : end do ! istm
1609 : end if ! istn is VB
1610 : end do ! istn
1611 : end do ! spins
1612 : end do ! k-points
1613 :
1614 8 : call xmpi_sum(inter2w,comm,ierr)
1615 8 : call xmpi_sum(inter1w,comm,ierr)
1616 8 : call xmpi_sum(intra2w,comm,ierr)
1617 8 : call xmpi_sum(intra1w,comm,ierr)
1618 8 : call xmpi_sum(intra1wS,comm,ierr)
1619 8 : call xmpi_min(my_emin,emin,comm,ierr)
1620 8 : call xmpi_max(my_emax,emax,comm,ierr)
1621 : ! Addition AR
1622 8 : if (.not.do_antiresonant) then
1623 4 : call xmpi_sum(inter2wa,comm,ierr)
1624 4 : call xmpi_sum(inter1wa,comm,ierr)
1625 4 : call xmpi_sum(intra2wa,comm,ierr)
1626 4 : call xmpi_sum(intra1wa,comm,ierr)
1627 4 : call xmpi_sum(intra1wSa,comm,ierr)
1628 : end if ! .not.do_antiresonant
1629 : ! Addition bands decomposition
1630 8 : if (do_decompo) then
1631 0 : if (contrib_decompo==0 .or. contrib_decompo==12) then
1632 0 : call xmpi_sum(inter2w_bands,comm,ierr)
1633 : end if
1634 0 : if (contrib_decompo==0 .or. contrib_decompo==11) then
1635 0 : call xmpi_sum(inter1w_bands,comm,ierr)
1636 : end if
1637 0 : if (contrib_decompo==0 .or. contrib_decompo==22) then
1638 0 : call xmpi_sum(intra2w_bands,comm,ierr)
1639 : end if
1640 0 : if (contrib_decompo==0 .or. contrib_decompo==21) then
1641 0 : call xmpi_sum(intra1w_bands,comm,ierr)
1642 : end if
1643 0 : if (contrib_decompo==0 .or. contrib_decompo==1) then
1644 0 : call xmpi_sum(intra1wS_bands,comm,ierr)
1645 : end if
1646 0 : if (contrib_decompo==0 .or. contrib_decompo==2) then
1647 : ! Addition 2bands interactions decomposition
1648 0 : call xmpi_sum(intra2w_2bands,comm,ierr)
1649 0 : call xmpi_sum(intra1wS_2bands,comm,ierr)
1650 : end if
1651 : end if ! do_decompo
1652 :
1653 8 : if (my_rank == master) then
1654 : ! write output in SI units and esu (esu to SI(m/v)=(value_esu)*(4xpi)/30000)
1655 :
1656 8 : if (ncid /= nctk_noid) then
1657 40 : start4 = [1, 1, icomp, itemp]
1658 40 : count4 = [2, nmesh, 1, 1]
1659 24 : ABI_MALLOC(chi2tot, (nmesh))
1660 2816 : chi2tot = inter2w + inter1w + intra2w + intra1w + intra1wS
1661 : ! Addition AR
1662 8 : if (.not.do_antiresonant) then
1663 8 : ABI_MALLOC(chi2tota, (nmesh))
1664 8 : ABI_MALLOC(chi2full, (nmesh))
1665 608 : chi2tota = inter2wa + inter1wa + intra2wa + intra1wa + intra1wSa
1666 608 : chi2full = chi2tot + chi2tota
1667 : end if ! .not.do_antiresonant
1668 8 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_inter2w"), c2r(inter2w), start=start4, count=count4))
1669 8 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_inter1w"), c2r(inter1w), start=start4, count=count4))
1670 8 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_intra2w"), c2r(intra2w), start=start4, count=count4))
1671 8 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_intra1w"), c2r(intra1w), start=start4, count=count4))
1672 8 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_intra1wS"), c2r(intra1wS), start=start4, count=count4))
1673 8 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_chi2tot"), c2r(chi2tot), start=start4, count=count4))
1674 : ! Addition AR
1675 8 : if (.not.do_antiresonant) then
1676 4 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_inter2w_AR"), c2r(inter2wa), start=start4, count=count4))
1677 4 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_inter1w_AR"), c2r(inter1wa), start=start4, count=count4))
1678 4 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_intra2w_AR"), c2r(intra2wa), start=start4, count=count4))
1679 4 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_intra1w_AR"), c2r(intra1wa), start=start4, count=count4))
1680 4 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_intra1wS_AR"), c2r(intra1wSa), start=start4, count=count4))
1681 4 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_chi2tot_AR"), c2r(chi2tota), start=start4, count=count4))
1682 4 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_chi2full"), c2r(chi2full), start=start4, count=count4))
1683 : end if ! .not.do_antiresonant
1684 :
1685 8 : ABI_FREE(chi2tot)
1686 : ! Addition AR
1687 8 : if (.not.do_antiresonant) then
1688 4 : ABI_FREE(chi2tota)
1689 4 : ABI_FREE(chi2full)
1690 : end if ! .not.do_antiresonant
1691 : end if
1692 :
1693 8 : if (open_file(fnam1,msg,newunit=fout1,action='WRITE',form='FORMATTED') /= 0) then
1694 0 : ABI_ERROR(msg)
1695 : end if
1696 8 : if (open_file(fnam2,msg,newunit=fout2,action='WRITE',form='FORMATTED') /= 0) then
1697 0 : ABI_ERROR(msg)
1698 : end if
1699 8 : if (open_file(fnam3,msg,newunit=fout3,action='WRITE',form='FORMATTED') /= 0) then
1700 0 : ABI_ERROR(msg)
1701 : end if
1702 8 : if (open_file(fnam4,msg,newunit=fout4,action='WRITE',form='FORMATTED') /= 0) then
1703 0 : ABI_ERROR(msg)
1704 : end if
1705 8 : if (open_file(fnam5,msg,newunit=fout5,action='WRITE',form='FORMATTED') /= 0) then
1706 0 : ABI_ERROR(msg)
1707 : end if
1708 8 : if (open_file(fnam6,msg,newunit=fout6,action='WRITE',form='FORMATTED') /= 0) then
1709 0 : ABI_ERROR(msg)
1710 : end if
1711 8 : if (open_file(fnam7,msg,newunit=fout7,action='WRITE',form='FORMATTED') /= 0) then
1712 0 : ABI_ERROR(msg)
1713 : end if
1714 : ! Addition bands decomposition
1715 8 : if (do_decompo) then
1716 0 : if (contrib_decompo==0 .or. contrib_decompo==12) then
1717 0 : if (open_file(fnam8,msg,newunit=fout8,action='WRITE',form='FORMATTED') /= 0) then
1718 0 : ABI_ERROR(msg)
1719 : end if
1720 : end if
1721 0 : if (contrib_decompo==0 .or. contrib_decompo==11) then
1722 0 : if (open_file(fnam9,msg,newunit=fout9,action='WRITE',form='FORMATTED') /= 0) then
1723 0 : ABI_ERROR(msg)
1724 : end if
1725 : end if
1726 0 : if (contrib_decompo==0 .or. contrib_decompo==22) then
1727 0 : if (open_file(fnam10,msg,newunit=fout10,action='WRITE',form='FORMATTED') /= 0) then
1728 0 : ABI_ERROR(msg)
1729 : end if
1730 : end if
1731 0 : if (contrib_decompo==0 .or. contrib_decompo==21) then
1732 0 : if (open_file(fnam11,msg,newunit=fout11,action='WRITE',form='FORMATTED') /= 0) then
1733 0 : ABI_ERROR(msg)
1734 : end if
1735 : end if
1736 0 : if (contrib_decompo==0 .or. contrib_decompo==1) then
1737 0 : if (open_file(fnam12,msg,newunit=fout12,action='WRITE',form='FORMATTED') /= 0) then
1738 0 : ABI_ERROR(msg)
1739 : end if
1740 : end if
1741 0 : if (contrib_decompo==0 .or. contrib_decompo==2) then
1742 : ! Addition 2bands interactions decomposition
1743 0 : if (open_file(fnam13,msg,newunit=fout13,action='WRITE',form='FORMATTED') /= 0) then
1744 0 : ABI_ERROR(msg)
1745 : end if
1746 0 : if (open_file(fnam14,msg,newunit=fout14,action='WRITE',form='FORMATTED') /= 0) then
1747 0 : ABI_ERROR(msg)
1748 : end if
1749 : end if
1750 : end if ! do_decompo
1751 :
1752 : ! write headers
1753 8 : write(fout1, '(a,3i3)' ) ' #calculated the component:',v1,v2,v3
1754 8 : write(fout1, '(a,es16.6)' ) ' #tolerance:',tol
1755 8 : write(fout1, '(a,es16.6,a)' ) ' #broadening:',brod,'Ha'
1756 8 : write(fout1, '(a,es16.6,a)' ) ' #scissors shift:',sc,'Ha'
1757 8 : write(fout1, '(a,es16.6,a,es16.6,a)' ) ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
1758 8 : write(fout1, '(a)' )' # Energy Tot-Im Chi(-2w,w,w) Tot-Im Chi(-2w,w,w)'
1759 8 : write(fout1, '(a)' )' # eV *10^-7 esu *10^-12 m/V SI units '
1760 8 : write(fout1, '(a)' )' # '
1761 :
1762 8 : write(fout2, '(a,3i3)' ) ' #calculated the component:',v1,v2,v3
1763 8 : write(fout2, '(a,es16.6)') ' #tolerance:',tol
1764 8 : write(fout2, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1765 8 : write(fout2, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1766 8 : write(fout2, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
1767 8 : write(fout2, '(a)')' # Energy Tot-Re Chi(-2w,w,w) Tot-Re Chi(-2w,w,w)'
1768 8 : write(fout2, '(a)')' # eV *10^-7 esu *10^-12 m/V SI units '
1769 8 : write(fout2, '(a)')' # '
1770 :
1771 8 : write(fout3, '(a,3i3)') ' #calculated the component:',v1,v2,v3
1772 8 : write(fout3, '(a,es16.6)') ' #tolerance:',tol
1773 8 : write(fout3, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1774 8 : write(fout3, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1775 8 : write(fout3, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
1776 8 : write(fout3, '(a)')' # Energy(eV) Inter(2w) inter(1w) intra(2w) intra(1w)'
1777 8 : write(fout3, '(a)')' # in esu'
1778 8 : write(fout3, '(a)')' # '
1779 :
1780 8 : write(fout4, '(a,3i3)') ' #calculated the component:',v1,v2,v3
1781 8 : write(fout4, '(a,es16.6)') ' #tolerance:',tol
1782 8 : write(fout4, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1783 8 : write(fout4, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1784 8 : write(fout4, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
1785 8 : write(fout4, '(a)')' # Energy(eV) Inter(2w) inter(1w) intra(2w) intra(1w)'
1786 8 : write(fout4, '(a)')' # in esu'
1787 8 : write(fout4, '(a)')' # '
1788 :
1789 8 : write(fout5, '(a,3i3)') ' #calculated the component:',v1,v2,v3
1790 8 : write(fout5, '(a,es16.6)') ' #tolerance:',tol
1791 8 : write(fout5, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1792 8 : write(fout5, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1793 8 : write(fout5, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
1794 8 : write(fout5, '(a)')' # Energy(eV) |TotChi(-2w,w,w)| |Tot Chi(-2w,w,w)|'
1795 8 : write(fout5, '(a)')' # eV *10^-7 esu *10^-12 m/V SI units '
1796 8 : write(fout5, '(a)')' # '
1797 :
1798 8 : write(fout6, '(a,3i3)') ' #calculated the component:',v1,v2,v3
1799 8 : write(fout6, '(a,es16.6)') ' #tolerance:',tol
1800 8 : write(fout6, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1801 8 : write(fout6, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1802 8 : write(fout6, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
1803 8 : write(fout6, '(a)')' # Energy(eV) Chi(w) Eta(w) Sigma(w)'
1804 8 : write(fout6, '(a)')' # in esu'
1805 8 : write(fout6, '(a)')' # '
1806 :
1807 8 : write(fout7, '(a,3i3)') ' #calculated the component:',v1,v2,v3
1808 8 : write(fout7, '(a,es16.6)') ' #tolerance:',tol
1809 8 : write(fout7, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1810 8 : write(fout7, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1811 8 : write(fout7, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
1812 8 : write(fout7, '(a)')' # Energy(eV) Chi(w) Eta(w) Sigma(w)'
1813 8 : write(fout7, '(a)')' # in esu'
1814 8 : write(fout7, '(a)')' # '
1815 : ! Addition bands decomposition
1816 8 : if (do_decompo) then
1817 0 : if (contrib_decompo==0 .or. contrib_decompo==12) then
1818 0 : write(fout8, '(a,3i3)') ' #calculated the component:',v1,v2,v3
1819 0 : write(fout8, '(a,es16.6)') ' #tolerance:',tol
1820 0 : write(fout8, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1821 0 : write(fout8, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1822 0 : write(fout8, '(a,es16.6,a)') ' #energy decompo input:',w_decompo,'eV'
1823 0 : write(fout8, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
1824 0 : write(fout8, '(a)')' # Energy(eV) n valence m conduction l both Re Inter2w Im Inter2w'
1825 0 : write(fout8, '(a)')' # eV *10^-12 m/V SI units *10^-12 m/V SI units '
1826 0 : write(fout8, '(a)')' # '
1827 : end if
1828 :
1829 0 : if (contrib_decompo==0 .or. contrib_decompo==11) then
1830 0 : write(fout9, '(a,3i3)') ' #calculated the component:',v1,v2,v3
1831 0 : write(fout9, '(a,es16.6)') ' #tolerance:',tol
1832 0 : write(fout9, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1833 0 : write(fout9, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1834 0 : write(fout9, '(a,es16.6,a)') ' #energy decompo input:',w_decompo,'eV'
1835 0 : write(fout9, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
1836 0 : write(fout9, '(a)')' # Energy(eV) n valence m conduction l both Re Inter1w Im Inter1w'
1837 0 : write(fout9, '(a)')' # eV *10^-12 m/V SI units *10^-12 m/V SI units '
1838 0 : write(fout9, '(a)')' # '
1839 : end if
1840 :
1841 0 : if (contrib_decompo==0 .or. contrib_decompo==22) then
1842 0 : write(fout10, '(a,3i3)') ' #calculated the component:',v1,v2,v3
1843 0 : write(fout10, '(a,es16.6)') ' #tolerance:',tol
1844 0 : write(fout10, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1845 0 : write(fout10, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1846 0 : write(fout10, '(a,es16.6,a)') ' #energy decompo input:',w_decompo,'eV'
1847 0 : write(fout10, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
1848 0 : write(fout10, '(a)')' # Energy(eV) n valence m conduction l both Re Intra2w Im Intra2w'
1849 0 : write(fout10, '(a)')' # eV *10^-12 m/V SI units *10^-12 m/V SI units '
1850 0 : write(fout10, '(a)')' # '
1851 : end if
1852 :
1853 0 : if (contrib_decompo==0 .or. contrib_decompo==21) then
1854 0 : write(fout11, '(a,3i3)') ' #calculated the component:',v1,v2,v3
1855 0 : write(fout11, '(a,es16.6)') ' #tolerance:',tol
1856 0 : write(fout11, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1857 0 : write(fout11, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1858 0 : write(fout11, '(a,es16.6,a)') ' #energy decompo input:',w_decompo,'eV'
1859 0 : write(fout11, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
1860 0 : write(fout11, '(a)')' # Energy(eV) n valence m conduction l both Re Intra1w Im Intra1w'
1861 0 : write(fout11, '(a)')' # eV *10^-12 m/V SI units *10^-12 m/V SI units '
1862 0 : write(fout11, '(a)')' # '
1863 : end if
1864 :
1865 0 : if (contrib_decompo==0 .or. contrib_decompo==1) then
1866 0 : write(fout12, '(a,3i3)') ' #calculated the component:',v1,v2,v3
1867 0 : write(fout12, '(a,es16.6)') ' #tolerance:',tol
1868 0 : write(fout12, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1869 0 : write(fout12, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1870 0 : write(fout12, '(a,es16.6,a)') ' #energy decompo input:',w_decompo,'eV'
1871 0 : write(fout12, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
1872 0 : write(fout12, '(a)')' # Energy(eV) n valence m conduction l both Re Intra1wS Im Intra1wS'
1873 0 : write(fout12, '(a)')' # eV *10^-12 m/V SI units *10^-12 m/V SI units '
1874 0 : write(fout12, '(a)')' # '
1875 : end if
1876 :
1877 0 : if (contrib_decompo==0 .or. contrib_decompo==2) then
1878 : ! Addition 2bands interactions decomposition
1879 0 : write(fout13, '(a,3i3)') ' #calculated the component:',v1,v2,v3
1880 0 : write(fout13, '(a,es16.6)') ' #tolerance:',tol
1881 0 : write(fout13, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1882 0 : write(fout13, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1883 0 : write(fout13, '(a,es16.6,a)') ' #energy decompo input:',w_decompo,'eV'
1884 0 : write(fout13, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
1885 0 : write(fout13, '(a)')' # Energy(eV) n valence m conduction Re Intra2w Im Intra2w'
1886 0 : write(fout13, '(a)')' # eV *10^-12 m/V SI units *10^-12 m/V SI units '
1887 0 : write(fout13, '(a)')' # '
1888 :
1889 0 : write(fout14, '(a,3i3)') ' #calculated the component:',v1,v2,v3
1890 0 : write(fout14, '(a,es16.6)') ' #tolerance:',tol
1891 0 : write(fout14, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1892 0 : write(fout14, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1893 0 : write(fout14, '(a,es16.6,a)') ' #energy decompo input:',w_decompo,'eV'
1894 0 : write(fout14, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
1895 0 : write(fout14, '(a)')' # Energy(eV) n valence m conduction Re Intra2w Im Intra2w'
1896 0 : write(fout14, '(a)')' # eV *10^-12 m/V SI units *10^-12 m/V SI units '
1897 0 : write(fout14, '(a)')' # '
1898 : end if
1899 : end if ! (do_decompo)
1900 :
1901 8 : totim=zero
1902 8 : totre=zero
1903 8 : totabs=zero
1904 8 : if (.not.do_antiresonant) then
1905 600 : do iw=2,nmesh
1906 596 : ene=(iw-1)*de
1907 596 : ene=ene*Ha_eV
1908 :
1909 : totim=aimag(inter2w(iw)+inter1w(iw)+intra2w(iw)+intra1w(iw)+intra1wS(iw) + &
1910 596 : inter2wa(iw)+inter1wa(iw)+intra2wa(iw)+intra1wa(iw)+intra1wSa(iw))/1.d-7
1911 596 : write(fout1,'(f15.6,2es15.6)') ene,totim,totim*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
1912 596 : totim=zero
1913 :
1914 : totre=dble(inter2w(iw)+inter1w(iw)+intra2w(iw)+intra1w(iw)+intra1wS(iw) + &
1915 596 : inter2wa(iw)+inter1wa(iw)+intra2wa(iw)+intra1wa(iw)+intra1wSa(iw))/1.d-7
1916 596 : write(fout2,'(f15.6,2es15.6)') ene,totre,totre*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
1917 596 : totre=zero
1918 :
1919 596 : write(fout3,'(f15.6,4es15.6)') ene,aimag(inter2w(iw)+inter2wa(iw))/1.d-7, &
1920 596 : aimag(inter1w(iw)+inter1wa(iw))/1.d-7,aimag(intra2w(iw)+intra2wa(iw))/1.d-7, &
1921 1192 : aimag(intra1w(iw)+intra1wS(iw)+intra1wa(iw)+intra1wSa(iw))/1.d-7
1922 :
1923 596 : write(fout4,'(f15.6,4es15.6)') ene,dble(inter2w(iw)+inter2wa(iw))/1.d-7, &
1924 596 : dble(inter1w(iw)+inter1wa(iw))/1.d-7,dble(intra2w(iw)+intra2wa(iw))/1.d-7, &
1925 1192 : dble(intra1w(iw)+intra1wS(iw)+intra1wa(iw)+intra1wSa(iw))/1.d-7
1926 :
1927 : totabs=abs(inter2w(iw)+inter1w(iw)+intra2w(iw)+intra1w(iw)+intra1wS(iw) + &
1928 596 : inter2wa(iw)+inter1wa(iw)+intra2wa(iw)+intra1wa(iw)+intra1wSa(iw))/1.d-7
1929 596 : write(fout5,'(f15.6,2es15.6)') ene,totabs,totabs*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
1930 596 : totabs=zero
1931 :
1932 596 : write(fout6,'(f15.6,4es15.6)') ene,aimag(inter2w(iw)+inter1w(iw)+ &
1933 596 : inter2wa(iw)+inter1wa(iw))/1.d-7, &
1934 596 : aimag(intra2w(iw)+intra1w(iw)+intra2wa(iw)+intra1wa(iw))/1.d-7, &
1935 1192 : aimag(intra1wS(iw)+intra1wSa(iw))/1.d-7
1936 :
1937 596 : write(fout7,'(f15.6,4es15.6)') ene,dble(inter2w(iw)+inter1w(iw) + &
1938 596 : inter2wa(iw)+inter1wa(iw))/1.d-7, &
1939 596 : dble(intra2w(iw)+intra1w(iw)+intra2wa(iw)+intra1wa(iw))/1.d-7, &
1940 1196 : dble(intra1wS(iw)+intra1wSa(iw))/1.d-7
1941 : end do ! iw=2, nmesh
1942 : else
1943 2200 : do iw=2,nmesh
1944 2196 : ene=(iw-1)*de
1945 2196 : ene=ene*Ha_eV
1946 :
1947 2196 : totim=aimag(inter2w(iw)+inter1w(iw)+intra2w(iw)+intra1w(iw)+intra1wS(iw))/1.d-7
1948 2196 : write(fout1,'(f15.6,2es15.6)') ene,totim,totim*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
1949 2196 : totim=zero
1950 :
1951 2196 : totre=dble(inter2w(iw)+inter1w(iw)+intra2w(iw)+intra1w(iw)+intra1wS(iw))/1.d-7
1952 2196 : write(fout2,'(f15.6,2es15.6)') ene,totre,totre*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
1953 2196 : totre=zero
1954 :
1955 2196 : write(fout3,'(f15.6,4es15.6)') ene,aimag(inter2w(iw))/1.d-7, &
1956 4392 : aimag(inter1w(iw))/1.d-7,aimag(intra2w(iw))/1.d-7, aimag(intra1w(iw)+intra1wS(iw))/1.d-7
1957 :
1958 2196 : write(fout4,'(f15.6,4es15.6)') ene,dble(inter2w(iw))/1.d-7, &
1959 4392 : dble(inter1w(iw))/1.d-7,dble(intra2w(iw))/1.d-7,dble(intra1w(iw)+intra1wS(iw))/1.d-7
1960 :
1961 2196 : totabs=abs(inter2w(iw)+inter1w(iw)+intra2w(iw)+intra1w(iw)+intra1wS(iw))/1.d-7
1962 2196 : write(fout5,'(f15.6,2es15.6)') ene,totabs,totabs*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
1963 2196 : totabs=zero
1964 :
1965 2196 : write(fout6,'(f15.6,4es15.6)') ene,aimag(inter2w(iw)+inter1w(iw))/1.d-7, &
1966 4392 : aimag(intra2w(iw)+intra1w(iw))/1.d-7,aimag(intra1wS(iw))/1.d-7
1967 :
1968 2196 : write(fout7,'(f15.6,4es15.6)') ene,dble(inter2w(iw)+inter1w(iw))/1.d-7, &
1969 4396 : dble(intra2w(iw)+intra1w(iw))/1.d-7,dble(intra1wS(iw))/1.d-7
1970 : end do ! iw=2, nmesh
1971 : end if ! .not.do_antiresonant, else
1972 : ! Addition bands decomposition
1973 8 : if (do_decompo) then
1974 0 : do iw=2,nmesh
1975 0 : if (iw==iw_tgt) then
1976 0 : ene=(iw-1)*de
1977 0 : ene=ene*Ha_eV
1978 :
1979 0 : do istn=1,nband_sum
1980 0 : do istm=1,nband_sum
1981 0 : do istl=1,nband_sum
1982 :
1983 0 : if (contrib_decompo==0 .or. contrib_decompo==12) then
1984 0 : write(fout8,'(f15.6,3i15,2es16.6)') ene,istn,istm,istl, &
1985 0 : dble(inter2w_bands(istl,istm,istn)) *4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-12), &
1986 0 : aimag(inter2w_bands(istl,istm,istn))*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-12)
1987 : end if
1988 :
1989 0 : if (contrib_decompo==0 .or. contrib_decompo==11) then
1990 0 : write(fout9,'(f15.6,3i15,2es16.6)') ene,istn,istm,istl, &
1991 0 : dble(inter1w_bands(istl,istm,istn)) *4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-12), &
1992 0 : aimag(inter1w_bands(istl,istm,istn))*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-12)
1993 : end if
1994 :
1995 0 : if (contrib_decompo==0 .or. contrib_decompo==22) then
1996 0 : write(fout10,'(f15.6,3i15,2es16.6)') ene,istn,istm,istl, &
1997 0 : dble(intra2w_bands(istl,istm,istn)) *4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-12), &
1998 0 : aimag(intra2w_bands(istl,istm,istn))*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-12)
1999 : end if
2000 :
2001 0 : if (contrib_decompo==0 .or. contrib_decompo==21) then
2002 0 : write(fout11,'(f15.6,3i15,2es16.6)') ene,istn,istm,istl, &
2003 0 : dble(intra1w_bands(istl,istm,istn)) *4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-12), &
2004 0 : aimag(intra1w_bands(istl,istm,istn))*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-12)
2005 : end if
2006 :
2007 0 : if (contrib_decompo==0 .or. contrib_decompo==1) then
2008 0 : write(fout12,'(f15.6,3i15,2es16.6)') ene,istn,istm,istl, &
2009 0 : dble(intra1wS_bands(istl,istm,istn)) *4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-12), &
2010 0 : aimag(intra1wS_bands(istl,istm,istn))*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-12)
2011 : end if
2012 :
2013 : end do ! istl
2014 0 : if (contrib_decompo==0 .or. contrib_decompo==2) then
2015 : ! Addition 2bands interactions decomposition
2016 0 : write(fout13,'(f15.6,2i15,2es16.6)') ene,istn,istm, &
2017 0 : dble(intra2w_2bands(istm,istn)) *4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-12), &
2018 0 : aimag(intra2w_2bands(istm,istn))*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-12)
2019 :
2020 0 : write(fout14,'(f15.6,2i15,2es16.6)') ene,istn,istm, &
2021 0 : dble(intra1wS_2bands(istm,istn)) *4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-12), &
2022 0 : aimag(intra1wS_2bands(istm,istn))*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-12)
2023 : end if
2024 :
2025 : end do ! istm
2026 : end do ! istn
2027 : end if ! iw==2
2028 : end do ! iw
2029 : end if ! do_decompo
2030 :
2031 :
2032 8 : close(fout1)
2033 8 : close(fout2)
2034 8 : close(fout3)
2035 8 : close(fout4)
2036 8 : close(fout5)
2037 8 : close(fout6)
2038 8 : close(fout7)
2039 : ! Addition bands decomposition
2040 8 : if (do_decompo) then
2041 0 : if (contrib_decompo==0 .or. contrib_decompo==12) then
2042 0 : close(fout8)
2043 : end if
2044 0 : if (contrib_decompo==0 .or. contrib_decompo==11) then
2045 0 : close(fout9)
2046 : end if
2047 0 : if (contrib_decompo==0 .or. contrib_decompo==22) then
2048 0 : close(fout10)
2049 : end if
2050 0 : if (contrib_decompo==0 .or. contrib_decompo==21) then
2051 0 : close(fout11)
2052 : end if
2053 0 : if (contrib_decompo==0 .or. contrib_decompo==1) then
2054 0 : close(fout12)
2055 : end if
2056 0 : if (contrib_decompo==0 .or. contrib_decompo==2) then
2057 : ! Addition 2bands interactions decomposition
2058 0 : close(fout13)
2059 0 : close(fout14)
2060 : end if
2061 : end if ! (do_decompo)
2062 :
2063 : ! print information
2064 8 : write(std_out,*) ' '
2065 8 : write(std_out,*) 'information about calculation just performed:'
2066 8 : write(std_out,*) ' '
2067 8 : write(std_out,*) 'calculated the component:',v1,v2,v3 ,'of second order susceptibility'
2068 8 : write(std_out,*) 'tolerance:',tol
2069 8 : if (tol.gt.0.008) write(std_out,*) 'ATTENTION: tolerance is too high'
2070 8 : write(std_out,*) 'broadening:',brod,'Hartree'
2071 8 : if (brod.gt.0.009) then
2072 0 : write(std_out,*) ' '
2073 0 : write(std_out,*) 'ATTENTION: broadening is quite high'
2074 0 : write(std_out,*) ' '
2075 : else if (brod.gt.0.015) then
2076 : write(std_out,*) ' '
2077 : write(std_out,*) 'ATTENTION: broadening is too high'
2078 : write(std_out,*) ' '
2079 : end if
2080 8 : write(std_out,*) 'scissors shift:',sc,'Hartree'
2081 8 : write(std_out,*) 'energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Hartree'
2082 : end if
2083 :
2084 : ! deallocate local arrays
2085 8 : ABI_FREE(px)
2086 8 : ABI_FREE(py)
2087 8 : ABI_FREE(pz)
2088 8 : ABI_FREE(inter2w)
2089 8 : ABI_FREE(inter1w)
2090 8 : ABI_FREE(intra2w)
2091 8 : ABI_FREE(intra1w)
2092 8 : ABI_FREE(intra1wS)
2093 8 : ABI_FREE(delta)
2094 : ! Addition AR
2095 8 : if (.not.do_antiresonant) then
2096 4 : ABI_FREE(inter2wa)
2097 4 : ABI_FREE(inter1wa)
2098 4 : ABI_FREE(intra2wa)
2099 4 : ABI_FREE(intra1wa)
2100 4 : ABI_FREE(intra1wSa)
2101 : end if ! .not.do_antiresonant
2102 : ! Addition bands decomposition
2103 8 : if (do_decompo) then
2104 0 : if (contrib_decompo==0 .or. contrib_decompo==12) then
2105 0 : ABI_FREE(inter2w_bands)
2106 0 : ABI_FREE(inter2w_bands_ik)
2107 : end if
2108 0 : if (contrib_decompo==0 .or. contrib_decompo==11) then
2109 0 : ABI_FREE(inter1w_bands)
2110 0 : ABI_FREE(inter1w_bands_ik)
2111 : end if
2112 0 : if (contrib_decompo==0 .or. contrib_decompo==22) then
2113 0 : ABI_FREE(intra2w_bands)
2114 0 : ABI_FREE(intra2w_bands_ik)
2115 : end if
2116 0 : if (contrib_decompo==0 .or. contrib_decompo==21) then
2117 0 : ABI_FREE(intra1w_bands)
2118 0 : ABI_FREE(intra1w_bands_ik)
2119 : end if
2120 0 : if (contrib_decompo==0 .or. contrib_decompo==1) then
2121 0 : ABI_FREE(intra1wS_bands)
2122 0 : ABI_FREE(intra1wS_bands_ik)
2123 : end if
2124 0 : if (contrib_decompo==0 .or. contrib_decompo==2) then
2125 : ! Addition 2bands interactions decomposition
2126 0 : ABI_FREE(intra2w_2bands)
2127 0 : ABI_FREE(intra2w_2bands_ik)
2128 0 : ABI_FREE(intra1wS_2bands)
2129 0 : ABI_FREE(intra1wS_2bands_ik)
2130 : end if
2131 : end if ! (do_decompo)
2132 :
2133 24 : end subroutine nlinopt
2134 : !!***
2135 :
2136 : !----------------------------------------------------------------------
2137 :
2138 : !!****f* m_optic_tools/linelop
2139 : !! NAME
2140 : !! linelop
2141 : !!
2142 : !! FUNCTION
2143 : !! Compute optical frequency dependent linear electro-optic susceptibility for semiconductors
2144 : !!
2145 : !! INPUTS
2146 : !! icomp=Sequential index associated to computed tensor components (used for netcdf output)
2147 : !! itemp=Temperature index (used for netcdf output)
2148 : !! nband_sum=Number of bands included in the sum. Must be <= mband
2149 : !! pmat(mband,mband,nkpt,3,nsppol) = momentum matrix elements in cartesian coordinates(complex)
2150 : !! v1,v2,v3 = desired component of the dielectric function(integer) 1=x,2=y,3=z
2151 : !! nmesh = desired number of energy mesh points(integer)
2152 : !! de = desired step in energy(real); nmesh*de=maximum energy for plotting
2153 : !! sc = scissors shift in Ha(real)
2154 : !! brod = broadening in Ha(real)
2155 : !! tol = tolerance:how close to the singularity exact exact is calculated(real)
2156 : !! fnam=root for filenames that will contain the output :
2157 : !! fnam1=trim(fnam)//'-ChiTotIm.out'
2158 : !! fnam2=trim(fnam)//'-ChiTotRe.out'
2159 : !! fnam3=trim(fnam)//'-ChiIm.out'
2160 : !! fnam4=trim(fnam)//'-ChiRe.out'
2161 : !! fnam5=trim(fnam)//'-ChiAbs.out'
2162 : !! ncid=Netcdf id to save output data.
2163 : !!
2164 : !! OUTPUT
2165 : !! Calculates the second harmonic generation susceptibility on a desired energy mesh and
2166 : !! for desired direction of polarisation. The output is in files named
2167 : !! ChiEOTot.out : Im\chi_{v1v2v3}(\omega,\omega,0) and Re\chi_{v1v2v3}(\omega,\omega,0)
2168 : !! ChiEOIm.out : contributions to the Im\chi_{v1v2v3}(\omega,\omega,0) from various terms
2169 : !! ChiEORe.out : contributions to Re\chi_{v1v2v3}(\omega,\omega,-0) from various terms
2170 : !! ChiEOAbs.out : abs\chi_{v1v2v3}(\omega,\omega,0). The headers in these files contain
2171 : !! information about the calculation.
2172 : !!
2173 : !! NOTES:
2174 : !! - The routine has been written using notations of Ref. 2
2175 : !! - This routine does not symmetrize the tensor (up to now)
2176 : !! - Sum over all the states and use occupation factors instead of looping only on resonant contributions
2177 : !!
2178 : !! SOURCE
2179 :
2180 5 : subroutine linelop(icomp, itemp, nband_sum, cryst, ks_ebands, &
2181 5 : pmat,v1,v2,v3,nmesh,de,sc,brod,tol,fnam,do_antiresonant,ncid,comm)
2182 :
2183 : !Arguments ------------------------------------
2184 : integer, intent(in) :: icomp, itemp, nband_sum, ncid
2185 : type(crystal_t),intent(in) :: cryst
2186 : type(ebands_t),intent(in) :: ks_ebands
2187 : complex(dp), intent(in) :: pmat(ks_ebands%mband, ks_ebands%mband, ks_ebands%nkpt, 3, ks_ebands%nsppol)
2188 : integer, intent(in) :: v1, v2, v3
2189 : integer, intent(in) :: nmesh
2190 : integer, intent(in) :: comm
2191 : real(dp), intent(in) :: de
2192 : real(dp), intent(in) :: sc
2193 : real(dp), intent(in) :: brod
2194 : real(dp), intent(in) :: tol
2195 : character(len=*), intent(in) :: fnam
2196 : logical, intent(in) :: do_antiresonant
2197 :
2198 : !Local variables -------------------------
2199 : integer,parameter :: master = 0
2200 : integer :: iw
2201 : integer :: i,j,k,lx,ly,lz
2202 : integer :: isp,isym,ik
2203 : integer :: ist1,istl,istn,istm, mband
2204 : real(dp) :: ene,totre,totabs,totim
2205 : real(dp) :: el,en,em
2206 : real(dp) :: emin,emax,my_emin,my_emax
2207 : real(dp) :: const_esu,const_au,au2esu
2208 : real(dp) :: wmn,wnm,wln,wnl,wml,wlm
2209 : complex(dp) :: idel,w,zi
2210 : character(len=fnlen) :: fnam1,fnam2,fnam3,fnam4,fnam5
2211 : ! local allocatable arrays
2212 5 : real(dp), allocatable :: s(:,:), sym(:,:,:)
2213 : integer :: start4(4),count4(4)
2214 : integer :: istp
2215 : real(dp) :: ep, wmp, wpn
2216 5 : real(dp), allocatable :: enk(:) ! (n) = \omega_n(k), with scissor included !
2217 : real(dp) :: fn, fm, fl, fnm, fnl, fml, fln, fmn
2218 5 : complex(dp), allocatable :: delta(:,:,:) ! (m,n,a) = \Delta_{mn}^{a}
2219 5 : complex(dp), allocatable :: rmna(:,:,:) ! (m,n,a) = r_{mn}^{a}
2220 5 : complex(dp), allocatable :: rmnbc(:,:,:,:) ! (m,n,b,c) = r^b_{mn;c}(k)
2221 5 : complex(dp), allocatable :: roverw(:,:,:,:) ! (m,n,b,c) = [r^b_{mn}(k)/w_{mn(k)];c
2222 5 : complex(dp), allocatable :: chi(:) ! \chi_{II}^{abc}(-\omega,\omega,0)
2223 5 : complex(dp), allocatable :: eta(:) ! \eta_{II}^{abc}(-\omega,\omega,0)
2224 5 : complex(dp), allocatable :: sigma(:) ! \frac{i}{\omega} \sigma_{II}^{abc}(-\omega,\omega,0)
2225 5 : complex(dp), allocatable :: chi2tot(:)
2226 : complex(dp) :: num1, num2, den1, den2, term1, term2
2227 : complex(dp) :: chi1, chi1_1, chi1_2, chi2_1b, chi2_2b
2228 5 : complex(dp), allocatable :: chi2(:) ! Second term that depends on the frequency ! (omega)
2229 : complex(dp) :: eta1, eta2, eta2_1, eta2_2
2230 : complex(dp) :: sigma1, sigma1_1, sigma1_2, sigma2
2231 : !Parallelism
2232 : integer :: my_rank, nproc, ierr, my_k1, my_k2
2233 : integer :: fout1,fout2,fout3,fout4,fout5
2234 : character(500) :: msg
2235 : ! *********************************************************************
2236 :
2237 5 : my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
2238 :
2239 : !calculate the constant
2240 5 : zi=(0._dp,1._dp)
2241 5 : idel=zi*brod
2242 : ! Disable symmetries for now
2243 5 : const_au=-2._dp/(cryst%ucvol*dble(cryst%nsym))
2244 5 : au2esu=5.8300348177d-8 ! REPLACE WITH DATA FROM DEFS_BASIS
2245 5 : const_esu=const_au*au2esu
2246 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2247 : !5.8300348177d-8 : au2esu : bohr*c*10^4/4pi*2*ry2ev
2248 : !bohr: 5.2917ifc nlinopt.f907E-11
2249 : !c: 2.99792458 velocity of light
2250 : !au2esu=(5.29177E-11*2.99792458*1.0E4)/Ha_eV
2251 : !this const includes (e^3*hbar^3*hbar^3)/(vol*hbar^5*m_e^3)
2252 : !mass comes from converting P_mn to r_mn
2253 : !hbar^3 comes from converting all frequencies to energies in denominator
2254 : !hbar^3 comes from operator for momentum (hbar/i nabla)
2255 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2256 : !output file names
2257 5 : fnam1=trim(fnam)//'-ChiEOTotIm.out'
2258 5 : fnam2=trim(fnam)//'-ChiEOTotRe.out'
2259 5 : fnam3=trim(fnam)//'-ChiEOIm.out'
2260 5 : fnam4=trim(fnam)//'-ChiEORe.out'
2261 5 : fnam5=trim(fnam)//'-ChiEOAbs.out'
2262 :
2263 : ! If there exists inversion symmetry exit with a mesg.
2264 5 : if (cryst%idx_spatial_inversion() /= 0) then
2265 0 : write(std_out,*) '-----------------------------------------'
2266 0 : write(std_out,*) ' the crystal has inversion symmetry '
2267 0 : write(std_out,*) ' the LEO susceptibility is zero '
2268 0 : write(std_out,*) '-----------------------------------------'
2269 0 : ABI_ERROR("Aborting now")
2270 : end if
2271 :
2272 : ! check polarisation
2273 5 : if (v1.le.0.or.v2.le.0.or.v3.le.0.or.v1.gt.3.or.v2.gt.3.or.v3.gt.3) then
2274 0 : write(std_out,*) '---------------------------------------------'
2275 0 : write(std_out,*) ' Error in linelop: '
2276 0 : write(std_out,*) ' the polarisation directions incorrect '
2277 0 : write(std_out,*) ' 1=x, 2=y and 3=z '
2278 0 : write(std_out,*) '---------------------------------------------'
2279 0 : ABI_ERROR("Aborting now")
2280 : end if
2281 :
2282 : ! number of energy mesh points
2283 5 : if (nmesh.le.0) then
2284 0 : write(std_out,*) '---------------------------------------------'
2285 0 : write(std_out,*) ' Error in linelop: '
2286 0 : write(std_out,*) ' number of energy mesh points incorrect '
2287 0 : write(std_out,*) ' number has to be integer greater than 0 '
2288 0 : write(std_out,*) ' nmesh*de = max energy for calculation '
2289 0 : write(std_out,*) '---------------------------------------------'
2290 0 : ABI_ERROR("Aborting now")
2291 : end if
2292 :
2293 : ! step in energy
2294 5 : if (de.le.zero) then
2295 0 : write(std_out,*) '---------------------------------------------'
2296 0 : write(std_out,*) ' Error in linelop: '
2297 0 : write(std_out,*) ' energy step is incorrect '
2298 0 : write(std_out,*) ' number has to real greater than 0.0 '
2299 0 : write(std_out,*) ' nmesh*de = max energy for calculation '
2300 0 : write(std_out,*) '---------------------------------------------'
2301 0 : ABI_ERROR("Aborting now")
2302 : end if
2303 :
2304 : ! broadening
2305 5 : if (brod.gt.0.009) then
2306 0 : write(std_out,*) '---------------------------------------------'
2307 0 : write(std_out,*) ' ATTENTION: broadening is quite high '
2308 0 : write(std_out,*) ' ideally should be less than 0.005 '
2309 0 : write(std_out,*) '---------------------------------------------'
2310 : else if (brod.gt.0.015) then
2311 : write(std_out,*) '----------------------------------------'
2312 : write(std_out,*) ' ATTENTION: broadening is too high '
2313 : write(std_out,*) ' ideally should be less than 0.005 '
2314 : write(std_out,*) '----------------------------------------'
2315 : end if
2316 :
2317 : ! tolerance
2318 5 : if (tol.gt.0.006) then
2319 0 : write(std_out,*) '----------------------------------------'
2320 0 : write(std_out,*) ' ATTENTION: tolerance is too high '
2321 0 : write(std_out,*) ' ideally should be less than 0.004 '
2322 0 : write(std_out,*) '----------------------------------------'
2323 : end if
2324 :
2325 5 : mband = ks_ebands%mband
2326 15 : ABI_MALLOC(enk, (mband))
2327 25 : ABI_MALLOC(delta, (mband, mband, 3))
2328 25 : ABI_MALLOC(rmnbc,(mband,mband, 3, 3))
2329 15 : ABI_MALLOC(roverw,(mband, mband, 3, 3))
2330 15 : ABI_MALLOC(rmna, (mband, mband, 3))
2331 15 : ABI_MALLOC(chi, (nmesh))
2332 10 : ABI_MALLOC(eta, (nmesh))
2333 10 : ABI_MALLOC(sigma, (nmesh))
2334 10 : ABI_MALLOC(chi2, (nmesh))
2335 5 : ABI_MALLOC(sym, (3, 3, 3))
2336 5 : ABI_MALLOC(s, (3, 3))
2337 :
2338 5 : ABI_CHECK(nband_sum <= mband, "nband_sum <= mband")
2339 :
2340 : ! generate the symmetrizing tensor
2341 200 : sym(:,:,:)=zero
2342 125 : do isym=1,cryst%nsym
2343 1560 : s(:,:)=cryst%symrel_cart(:,:,isym)
2344 485 : do i=1,3
2345 1560 : do j=1,3
2346 4680 : do k=1,3
2347 4320 : sym(i,j,k)=sym(i,j,k)+(s(i,v1)*s(j,v2)*s(k,v3))
2348 : end do
2349 : end do
2350 : end do
2351 : end do
2352 :
2353 : ! initialise
2354 5330 : delta(:,:,:)=zero
2355 15995 : rmnbc(:,:,:,:)=zero
2356 1605 : chi(:)=zero
2357 1605 : chi2(:) = zero
2358 1605 : eta(:)=zero
2359 1605 : sigma(:)=zero
2360 5 : my_emin=HUGE(zero)
2361 5 : my_emax=-HUGE(zero)
2362 :
2363 : ! Split work
2364 5 : call xmpi_split_work(ks_ebands%nkpt,comm,my_k1,my_k2)
2365 :
2366 : ! loop over kpts
2367 261 : do ik=my_k1,my_k2
2368 256 : write(std_out,*) "P-",my_rank,": ",ik,'of',ks_ebands%nkpt
2369 517 : do isp=1,ks_ebands%nsppol
2370 : ! Calculate the scissor corrected energies and the energy window
2371 3968 : do ist1=1,nband_sum
2372 3712 : en = ks_ebands%eig(ist1,ik,isp)
2373 3712 : my_emin=min(my_emin,en)
2374 3712 : my_emax=max(my_emax,en)
2375 3712 : if(en > ks_ebands%fermie) then
2376 2688 : en = en + sc
2377 : end if
2378 3968 : enk(ist1) = en
2379 : end do
2380 :
2381 : ! calculate \Delta_nm and r_mn^a
2382 3968 : do istn=1,nband_sum
2383 3712 : en = enk(istn)
2384 65536 : do istm=1,nband_sum
2385 61568 : em = enk(istm)
2386 61568 : wmn = em - en
2387 246272 : delta(istn,istm,1:3)=pmat(istn,istn,ik,1:3,isp)-pmat(istm,istm,ik,1:3,isp)
2388 65280 : if(abs(wmn) < tol) then
2389 16576 : rmna(istm,istn,1:3) = zero
2390 : else
2391 229696 : rmna(istm,istn,1:3)=-zi*pmat(istm,istn,ik,1:3,isp)/wmn
2392 : end if
2393 : end do
2394 : end do
2395 :
2396 : ! calculate \r^b_mn;c
2397 3968 : do istm=1,nband_sum
2398 3712 : em = enk(istm)
2399 65536 : do istn=1,nband_sum
2400 61568 : en = enk(istn)
2401 61568 : wmn = em - en
2402 65280 : if(abs(wmn) > tol) then
2403 229696 : do ly = 1,3
2404 746512 : do lz = 1,3
2405 516816 : num1 = (rmna(istm,istn,ly)*delta(istm,istn,lz))+(rmna(istm,istn,lz)*delta(istm,istn,ly))
2406 516816 : den1 = wmn
2407 516816 : term1 = num1/den1
2408 516816 : term2 = zero
2409 9945504 : do istp=1,nband_sum
2410 9428688 : ep = enk(istp)
2411 9428688 : wmp = em - ep
2412 9428688 : wpn = ep - en
2413 9428688 : num2 = (wmp*rmna(istm,istp,ly)*rmna(istp,istn,lz))-(wpn*rmna(istm,istp,lz)*rmna(istp,istn,ly))
2414 9428688 : den2 = wmn
2415 9945504 : term2 = term2 + (num2/den2)
2416 : end do
2417 516816 : rmnbc(istm,istn,ly,lz) = -term1-(zi*term2)
2418 689088 : roverw(istm,istn,ly,lz) = (rmnbc(istm,istn,ly,lz)/wmn) - (rmna(istm,istn,ly)/(wmn**2))*delta(istm,istn,lz)
2419 : end do
2420 : end do
2421 : end if
2422 : end do
2423 : end do
2424 :
2425 : ! initialise the factors
2426 : ! start the calculation
2427 4224 : do istn=1,nband_sum
2428 3712 : en=enk(istn)
2429 3712 : if (do_antiresonant .and. en .ge. ks_ebands%fermie) then
2430 : cycle
2431 : end if
2432 3072 : fn=ks_ebands%occ(istn,ik,isp)
2433 59136 : do istm=1,nband_sum
2434 55808 : em=enk(istm)
2435 55808 : if (do_antiresonant .and. em .le. ks_ebands%fermie) then
2436 : cycle
2437 : end if
2438 53760 : wmn=em-en
2439 53760 : wnm=-wmn
2440 53760 : fm = ks_ebands%occ(istm,ik,isp)
2441 53760 : fnm = fn - fm
2442 53760 : fmn = fm - fn
2443 53760 : eta1 = zero
2444 53760 : eta2_1 = zero
2445 53760 : eta2_2 = zero
2446 53760 : sigma1_1 = zero
2447 53760 : sigma1_2 = zero
2448 53760 : sigma2 = zero
2449 53760 : if(abs(wmn) > tol) then
2450 203264 : do lx = 1,3
2451 660608 : do ly = 1,3
2452 1981824 : do lz = 1,3
2453 1372032 : eta1 = eta1 + sym(lx,ly,lz)*(fnm*rmna(istn,istm,lx)*(roverw(istm,istn,lz,ly)))
2454 1372032 : eta2_1 = eta2_1 + sym(lx,ly,lz)*(fnm*(rmna(istn,istm,lx)*rmnbc(istm,istn,ly,lz)))
2455 1372032 : eta2_2 = eta2_2 + sym(lx,ly,lz)*(fnm*(rmnbc(istn,istm,lx,lz)*rmna(istm,istn,ly)))
2456 1372032 : sigma1_1 = sigma1_1 + sym(lx,ly,lz)*(fnm*delta(istn,istm,lx)*rmna(istn,istm,ly)*rmna(istm,istn,lz))/(wmn**2)
2457 1372032 : sigma1_2 = sigma1_2 + sym(lx,ly,lz)*(fnm*delta(istn,istm,lx)*rmna(istn,istm,lz)*rmna(istm,istn,ly))/(wmn**2)
2458 1829376 : sigma2 = sigma2 + sym(lx,ly,lz)*(fnm*rmnbc(istn,istm,lz,lx)*rmna(istm,istn,ly))/wmn
2459 : end do
2460 : end do
2461 : end do
2462 : end if
2463 53760 : chi1_1 = zero
2464 53760 : chi1_2 = zero
2465 53760 : chi2_1b = zero
2466 53760 : chi2_2b = zero
2467 10293760 : chi2(:) = zero
2468 : ! Three band terms
2469 1100800 : do istl=1,nband_sum
2470 1047040 : el=enk(istl)
2471 1047040 : fl = ks_ebands%occ(istl,ik,isp)
2472 1047040 : wlm = el-em
2473 1047040 : wln = el-en
2474 1047040 : wnl = en-el
2475 1047040 : wml = em-el
2476 1047040 : fnl = fn-fl
2477 1047040 : fln = fl-fn
2478 1047040 : fml = fm-fl
2479 4241920 : do lx = 1,3
2480 13611520 : do ly = 1,3
2481 40834560 : do lz = 1,3
2482 28270080 : if(abs(wlm) > tol) then
2483 26607744 : chi1_1 = chi1_1 + sym(lx,ly,lz)*(fnm*rmna(istn,istm,lx)*rmna(istm,istl,lz)*rmna(istl,istn,ly))/(wlm)
2484 26607744 : chi2_1b = chi2_1b + sym(lx,ly,lz)*(fnm*rmna(istn,istl,lx)*rmna(istl,istm,lz)*rmna(istm,istn,ly))/(wlm)
2485 : end if
2486 37693440 : if(abs(wln) > tol) then
2487 26609040 : chi1_2 = chi1_2 + sym(lx,ly,lz)*(fnm*rmna(istn,istm,lx)*rmna(istm,istl,ly)*rmna(istl,istn,lz))/(wln)
2488 26609040 : chi2_2b = chi2_2b + sym(lx,ly,lz)*(fmn*rmna(istl,istm,lx)*rmna(istm,istn,ly)*rmna(istn,istl,lz))/(wnl)
2489 : end if
2490 : end do
2491 : end do
2492 : end do
2493 : end do
2494 :
2495 53760 : sigma1 = 0.5_dp*(sigma1_1-sigma1_2)
2496 53760 : eta2 = 0.5_dp*(eta2_1-eta2_2)
2497 53760 : chi1 = chi1_1 + chi1_2
2498 :
2499 : ! calculate over the desired energy mesh and sum over k-points
2500 10297472 : do iw=1,nmesh
2501 10240000 : w=(iw-1)*de+idel
2502 : ! Better way to compute it
2503 10240000 : chi(iw) = chi(iw) + 0.5_dp*ks_ebands%wtk(ik)*((chi1/(wmn-w)) + ((chi2_1b+chi2_2b)/(wmn-w)))*const_esu
2504 10240000 : eta(iw) = eta(iw) + 0.5_dp*zi*ks_ebands%wtk(ik)*((eta1/(wmn-w)) + (eta2/((wmn-w)**2)))*const_esu
2505 10295808 : sigma(iw) = sigma(iw) + 0.5_dp*zi*ks_ebands%wtk(ik)*((sigma1/(wmn-w))- (sigma2/(wmn-w)))*const_esu
2506 : end do
2507 : end do ! istn and istm
2508 : end do
2509 : end do ! spins
2510 : end do ! k-points
2511 :
2512 5 : call xmpi_sum(chi,comm,ierr)
2513 5 : call xmpi_sum(eta,comm,ierr)
2514 5 : call xmpi_sum(sigma,comm,ierr)
2515 5 : call xmpi_min(my_emin,emin,comm,ierr)
2516 5 : call xmpi_max(my_emax,emax,comm,ierr)
2517 :
2518 5 : if (my_rank == master) then
2519 :
2520 5 : if (ncid /= nctk_noid) then
2521 25 : start4 = [1, 1, icomp, itemp]
2522 25 : count4 = [2, nmesh, 1, 1]
2523 10 : ABI_MALLOC(chi2tot, (nmesh))
2524 1610 : chi2tot = chi + eta + sigma
2525 5 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo_chi"), c2r(chi), start=start4, count=count4))
2526 5 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo_eta"), c2r(eta), start=start4, count=count4))
2527 5 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo_sigma"), c2r(sigma), start=start4, count=count4))
2528 5 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo_chi2tot"), c2r(chi2tot), start=start4, count=count4))
2529 5 : ABI_FREE(chi2tot)
2530 : end if
2531 :
2532 : ! write output in SI units and esu (esu to SI(m/v)=(value_esu)*(4xpi)/30000)
2533 5 : if (open_file(fnam1,msg,newunit=fout1,action='WRITE',form='FORMATTED') /= 0) then
2534 0 : ABI_ERROR(msg)
2535 : end if
2536 5 : if (open_file(fnam2,msg,newunit=fout2,action='WRITE',form='FORMATTED') /= 0) then
2537 0 : ABI_ERROR(msg)
2538 : end if
2539 5 : if (open_file(fnam3,msg,newunit=fout3,action='WRITE',form='FORMATTED') /= 0) then
2540 0 : ABI_ERROR(msg)
2541 : end if
2542 5 : if (open_file(fnam4,msg,newunit=fout4,action='WRITE',form='FORMATTED') /= 0) then
2543 0 : ABI_ERROR(msg)
2544 : end if
2545 5 : if (open_file(fnam5,msg,newunit=fout5,action='WRITE',form='FORMATTED') /= 0) then
2546 0 : ABI_ERROR(msg)
2547 : end if
2548 : ! write headers
2549 5 : write(fout1, '(a,3i3)' ) ' #calculated the component:',v1,v2,v3
2550 5 : write(fout1, '(a,es16.6)' ) ' #tolerance:',tol
2551 5 : write(fout1, '(a,es16.6,a)' ) ' #broadening:',brod,'Ha'
2552 5 : write(fout1, '(a,es16.6,a)' ) ' #scissors shift:',sc,'Ha'
2553 5 : write(fout1, '(a,es16.6,a,es16.6,a)' ) ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
2554 5 : write(fout1, '(a)' )' # Energy Tot-Im Chi(-w,w,0) Tot-Im Chi(-w,w,0)'
2555 5 : write(fout1, '(a)' )' # eV *10^-7 esu *10^-12 m/V SI units '
2556 5 : write(fout1, '(a)' )' # '
2557 :
2558 5 : write(fout2, '(a,3i3)' ) ' #calculated the component:',v1,v2,v3
2559 5 : write(fout2, '(a,es16.6)') ' #tolerance:',tol
2560 5 : write(fout2, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
2561 5 : write(fout2, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
2562 5 : write(fout2, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
2563 5 : write(fout2, '(a)')' # Energy Tot-Re Chi(-w,w,0) Tot-Re Chi(-w,w,0)'
2564 5 : write(fout2, '(a)')' # eV *10^-7 esu *10^-12 m/V SI units '
2565 5 : write(fout2, '(a)')' # '
2566 :
2567 5 : write(fout3, '(a,3i3)') ' #calculated the component:',v1,v2,v3
2568 5 : write(fout3, '(a,es16.6)') ' #tolerance:',tol
2569 5 : write(fout3, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
2570 5 : write(fout3, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
2571 5 : write(fout3, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
2572 5 : write(fout3, '(a)')' # Energy(eV) Chi(w) Eta(w) Sigma(w)'
2573 5 : write(fout3, '(a)')' # in esu'
2574 5 : write(fout3, '(a)')' # '
2575 :
2576 5 : write(fout4, '(a,3i3)') ' #calculated the component:',v1,v2,v3
2577 5 : write(fout4, '(a,es16.6)') ' #tolerance:',tol
2578 5 : write(fout4, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
2579 5 : write(fout4, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
2580 5 : write(fout4, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
2581 5 : write(fout4, '(a)')' # Energy(eV) Chi(w) Eta(w) Sigma(w)'
2582 5 : write(fout4, '(a)')' # in esu'
2583 5 : write(fout4, '(a)')' # '
2584 :
2585 5 : write(fout5, '(a,3i3)') ' #calculated the component:',v1,v2,v3
2586 5 : write(fout5, '(a,es16.6)') ' #tolerance:',tol
2587 5 : write(fout5, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
2588 5 : write(fout5, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
2589 5 : write(fout5, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
2590 5 : write(fout5, '(a)')' # Energy(eV) |TotChi(-w,w,0)| |Tot Chi(-w,w,0)|'
2591 5 : write(fout5, '(a)')' # eV *10^-7 esu *10^-12 m/V SI units '
2592 5 : write(fout5, '(a)')' # '
2593 :
2594 5 : totim=zero
2595 5 : totre=zero
2596 5 : totabs=zero
2597 1600 : do iw=2,nmesh
2598 1595 : ene=(iw-1)*de
2599 1595 : ene=ene*Ha_eV
2600 1595 : totim=aimag(chi(iw)+eta(iw)+sigma(iw))/1.d-7
2601 1595 : write(fout1,'(f15.6,2es15.6)') ene,totim,totim*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
2602 1595 : totim=zero
2603 1595 : totre=dble(chi(iw)+eta(iw)+sigma(iw))/1.d-7
2604 1595 : write(fout2,'(f15.6,2es15.6)') ene,totre,totre*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
2605 1595 : totre=zero
2606 1595 : write(fout3,'(f15.6,3es15.6)') ene,aimag(chi(iw))/1.d-7, &
2607 3190 : aimag(eta(iw))/1.d-7,aimag(sigma(iw))/1.d-7
2608 1595 : write(fout4,'(f15.6,3es15.6)') ene,dble(chi(iw))/1.d-7, &
2609 3190 : dble(eta(iw))/1.d-7,dble(sigma(iw))/1.d-7
2610 1595 : totabs=abs(chi(iw)+eta(iw)+sigma(iw))/1.d-7
2611 1595 : write(fout5,'(f15.6,2es15.6)') ene,totabs,totabs*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
2612 1600 : totabs=zero
2613 : end do
2614 :
2615 5 : close(fout1)
2616 5 : close(fout2)
2617 5 : close(fout3)
2618 5 : close(fout4)
2619 5 : close(fout5)
2620 : ! print information
2621 5 : write(std_out,*) ' '
2622 5 : write(std_out,*) 'information about calculation just performed:'
2623 5 : write(std_out,*) ' '
2624 5 : write(std_out,*) 'calculated the component:',v1,v2,v3 ,'of LEO susceptibility'
2625 5 : write(std_out,*) 'tolerance:',tol
2626 5 : if (tol.gt.0.008) write(std_out,*) 'ATTENTION: tolerance is too high'
2627 5 : write(std_out,*) 'broadening:',brod,'Hartree'
2628 5 : if (brod.gt.0.009) then
2629 0 : write(std_out,*) ' '
2630 0 : write(std_out,*) 'ATTENTION: broadening is quite high'
2631 0 : write(std_out,*) ' '
2632 : else if (brod.gt.0.015) then
2633 : write(std_out,*) ' '
2634 : write(std_out,*) 'ATTENTION: broadening is too high'
2635 : write(std_out,*) ' '
2636 : end if
2637 5 : write(std_out,*) 'scissors shift:',sc,'Hartree'
2638 5 : write(std_out,*) 'energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Hartree'
2639 :
2640 : end if
2641 :
2642 : ! deallocate local arrays
2643 5 : ABI_FREE(enk)
2644 5 : ABI_FREE(delta)
2645 5 : ABI_FREE(rmnbc)
2646 5 : ABI_FREE(roverw)
2647 5 : ABI_FREE(rmna)
2648 5 : ABI_FREE(chi)
2649 5 : ABI_FREE(chi2)
2650 5 : ABI_FREE(eta)
2651 5 : ABI_FREE(sigma)
2652 5 : ABI_FREE(s)
2653 5 : ABI_FREE(sym)
2654 :
2655 15 : end subroutine linelop
2656 : !!***
2657 :
2658 : !----------------------------------------------------------------------
2659 :
2660 : !!****f* m_optic_tools/nonlinopt
2661 : !! NAME
2662 : !! nonlinopt
2663 : !!
2664 : !! FUNCTION
2665 : !! Compute the frequency dependent nonlinear electro-optic susceptibility for semiconductors
2666 : !!
2667 : !! INPUTS
2668 : !! icomp=Sequential index associated to computed tensor components (used for netcdf output)
2669 : !! itemp=Temperature index (used for netcdf output)
2670 : !! nband_sum=Number of bands included in the sum. Must be <= mband
2671 : !! pmat(mband,mband,nkpt,3,nsppol) = momentum matrix elements in cartesian coordinates(complex)
2672 : !! v1,v2,v3 = desired component of the dielectric function(integer) 1=x,2=y,3=z
2673 : !! nmesh = desired number of energy mesh points(integer)
2674 : !! de = desired step in energy(real); nmesh*de=maximum energy for plotting
2675 : !! sc = scissors shift in Ha(real)
2676 : !! brod = broadening in Ha(real)
2677 : !! tol = tolerance:how close to the singularity exact exact is calculated(real)
2678 : !! fnam=root for filenames that will contain the output :
2679 : !! fnam1=trim(fnam)//'-ChiTotIm.out'
2680 : !! fnam2=trim(fnam)//'-ChiTotRe.out'
2681 : !! fnam3=trim(fnam)//'-ChiIm.out'
2682 : !! fnam4=trim(fnam)//'-ChiRe.out'
2683 : !! fnam5=trim(fnam)//'-ChiAbs.out'
2684 : !!
2685 : !! OUTPUT
2686 : !! Calculates the nonlinear electro-optical susceptibility on a desired energy mesh and
2687 : !! for desired direction of polarisation. The output is in files named
2688 : !! ChiEOTot.out : Im\chi_{v1v2v3}(\omega,\omega,0) and Re\chi_{v1v2v3}(\omega,\omega,0)
2689 : !! ChiEOIm.out : contributions to the Im\chi_{v1v2v3}(\omega,\omega,0) from various terms
2690 : !! ChiEORe.out : contributions to Re\chi_{v1v2v3}(\omega,\omega,-0) from various terms
2691 : !! ChiEOAbs.out : abs\chi_{v1v2v3}(\omega,\omega,0). The headers in these files contain
2692 : !! information about the calculation.
2693 : !! ncid=Netcdf id to save output data.
2694 : !!
2695 : !! COMMENTS
2696 : !! - The routine has been written using notations of Ref. 2
2697 : !! - This routine does not symmetrize the tensor (up to now)
2698 : !! - Sum over all the states and use occupation factors instead of looping only on resonant contributions
2699 : !!
2700 : !! SOURCE
2701 :
2702 4 : subroutine nonlinopt(icomp, itemp, nband_sum, cryst, ks_ebands, &
2703 4 : pmat, v1, v2, v3, nmesh, de, sc, brod, tol, fnam, do_antiresonant, ncid, comm)
2704 :
2705 : !Arguments ------------------------------------
2706 : integer, intent(in) :: icomp, itemp, nband_sum, ncid
2707 : type(crystal_t),intent(in) :: cryst
2708 : type(ebands_t),intent(in) :: ks_ebands
2709 : complex(dp), intent(in) :: pmat(ks_ebands%mband, ks_ebands%mband, ks_ebands%nkpt, 3, ks_ebands%nsppol)
2710 : integer, intent(in) :: v1, v2, v3
2711 : integer, intent(in) :: nmesh
2712 : integer, intent(in) :: comm
2713 : real(dp), intent(in) :: de, sc, brod, tol
2714 : character(len=*), intent(in) :: fnam
2715 : logical, intent(in) :: do_antiresonant
2716 :
2717 : !Local variables -------------------------
2718 : integer :: iw,i,j,k,lx,ly,lz,mband
2719 : integer :: isp,isym,ik,ist1,istl,istn,istm
2720 : real(dp) :: ene,totre,totabs,totim
2721 : real(dp) :: el,en,em
2722 : real(dp) :: emin,emax, my_emin,my_emax
2723 : real(dp) :: const_esu,const_au,au2esu
2724 : real(dp) :: wmn,wnm,wln,wnl,wml,wlm !, t1
2725 : complex(dp) :: idel,w,zi
2726 : character(len=fnlen) :: fnam1,fnam2,fnam3,fnam4,fnam5,fnam6,fnam7
2727 : ! local allocatable arrays
2728 : integer :: start4(4),count4(4)
2729 : real(dp) :: s(3,3),sym(3,3,3)
2730 : integer :: istp
2731 : real(dp) :: ep, wmp, wpn, wtk
2732 4 : real(dp), allocatable :: enk(:) ! (n) = \omega_n(k), with scissor included !
2733 : real(dp) :: fn, fm, fl, fnm, fnl, fml, fln, flm
2734 4 : complex(dp), allocatable :: delta(:,:,:) ! (m,n,a) = \Delta_{mn}^{a}
2735 4 : complex(dp), allocatable :: rmna(:,:,:) ! (m,n,a) = r_{mn}^{a}
2736 4 : complex(dp), allocatable :: rmnbc(:,:,:,:) ! (m,n,b,c) = r^b_{mn;c}(k)
2737 4 : complex(dp), allocatable :: roverw(:,:,:,:) ! (m,n,b,c) = [r^b_{mn}(k)/w_{mn(k)];c
2738 4 : complex(dp), allocatable :: chiw(:), chi2w(:) ! \chi_{II}^{abc}(-\omega,\omega,0)
2739 4 : complex(dp), allocatable :: etaw(:), eta2w(:) ! \eta_{II}^{abc}(-\omega,\omega,0)
2740 4 : complex(dp), allocatable :: sigmaw(:) ! \frac{i}{\omega} \sigma_{II}^{abc}(-\omega,\omega,0)
2741 : complex(dp) :: num1, num2, den1, den2, term1, term2
2742 : complex(dp) :: chi1, chi2_1, chi2_2
2743 4 : complex(dp), allocatable :: chi2(:) ! Second term that depends on the frequency ! (omega)
2744 4 : complex(dp), allocatable :: eta1(:) ! Second term that depends on the frequency ! (omega)
2745 4 : complex(dp), allocatable :: chi2tot(:)
2746 : complex(dp) :: eta1_1, eta1_2, eta2_1, eta2_2
2747 : complex(dp) :: sigma2_1, sigma1
2748 4 : complex(dp), allocatable :: symrmn(:,:,:) ! (m,l,n) = 1/2*(rml^b rln^c+rml^c rln^b)
2749 : complex(dp) :: symrmnl(3,3), symrlmn(3,3), symrmln(3,3)
2750 : !Parallelism
2751 : integer :: my_rank, nproc
2752 : integer,parameter :: master = 0
2753 : integer :: ierr
2754 : integer :: my_k1, my_k2
2755 : character(500) :: msg
2756 : integer :: fout1,fout2,fout3,fout4,fout5,fout6,fout7
2757 :
2758 : ! *********************************************************************
2759 :
2760 4 : my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
2761 :
2762 : !calculate the constant
2763 4 : zi=(0._dp,1._dp)
2764 4 : idel=zi*brod
2765 4 : const_au=-2._dp/(cryst%ucvol*dble(cryst%nsym))
2766 : !const_au=-2._dp/(cryst%ucvol)
2767 4 : au2esu=5.8300348177d-8 ! REPLACE WITH DATA FROM DEFS_BASIS
2768 4 : const_esu=const_au*au2esu
2769 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2770 : !5.8300348177d-8 : au2esu : bohr*c*10^4/4pi*2*ry2ev
2771 : !bohr: 5.2917ifc nlinopt.f907E-11
2772 : !c: 2.99792458 velocity of light
2773 : !au2esu=(5.29177E-11*2.99792458*1.0E4)/Ha_eV
2774 : !this const includes (e^3*hbar^3*hbar^3)/(vol*hbar^5*m_e^3)
2775 : !mass comes from converting P_mn to r_mn
2776 : !hbar^3 comes from converting all frequencies to energies in denominator
2777 : !hbar^3 comes from operator for momentum (hbar/i nabla)
2778 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2779 : !output file names
2780 4 : fnam1=trim(fnam)//'-ChiSHGTotIm.out'
2781 4 : fnam2=trim(fnam)//'-ChiSHGTotRe.out'
2782 4 : fnam3=trim(fnam)//'-ChiSHGIm.out'
2783 4 : fnam4=trim(fnam)//'-ChiSHGRe.out'
2784 4 : fnam5=trim(fnam)//'-ChiSHGAbs.out'
2785 4 : fnam6=trim(fnam)//'-ChiSHGImDec.out'
2786 4 : fnam7=trim(fnam)//'-ChiSHGReDec.out'
2787 :
2788 : ! If there exists inversion symmetry exit with a mesg.
2789 4 : if (cryst%idx_spatial_inversion() /= 0) then
2790 0 : write(std_out,*) '-----------------------------------------'
2791 0 : write(std_out,*) ' the crystal has inversion symmetry '
2792 0 : write(std_out,*) ' the nl electro-optical susceptibility'
2793 0 : write(std_out,*) ' is zero '
2794 0 : write(std_out,*) '-----------------------------------------'
2795 0 : ABI_ERROR("Aborting now")
2796 : end if
2797 :
2798 : ! check polarisation
2799 4 : if (v1.le.0.or.v2.le.0.or.v3.le.0.or.v1.gt.3.or.v2.gt.3.or.v3.gt.3) then
2800 0 : write(std_out,*) '---------------------------------------------'
2801 0 : write(std_out,*) ' Error in nonlinopt: '
2802 0 : write(std_out,*) ' the polarisation directions incorrect '
2803 0 : write(std_out,*) ' 1=x, 2=y and 3=z '
2804 0 : write(std_out,*) '---------------------------------------------'
2805 0 : ABI_ERROR("Aborting now")
2806 : end if
2807 :
2808 : ! number of energy mesh points
2809 4 : if (nmesh.le.0) then
2810 0 : write(std_out,*) '---------------------------------------------'
2811 0 : write(std_out,*) ' Error in nonlinopt: '
2812 0 : write(std_out,*) ' number of energy mesh points incorrect '
2813 0 : write(std_out,*) ' number has to be integer greater than 0 '
2814 0 : write(std_out,*) ' nmesh*de = max energy for calculation '
2815 0 : write(std_out,*) '---------------------------------------------'
2816 0 : ABI_ERROR("Aborting now")
2817 : end if
2818 :
2819 : ! step in energy
2820 4 : if (de.le.zero) then
2821 0 : write(std_out,*) '---------------------------------------------'
2822 0 : write(std_out,*) ' Error in nonlinopt: '
2823 0 : write(std_out,*) ' energy step is incorrect '
2824 0 : write(std_out,*) ' number has to real greater than 0.0 '
2825 0 : write(std_out,*) ' nmesh*de = max energy for calculation '
2826 0 : write(std_out,*) '---------------------------------------------'
2827 0 : ABI_ERROR("Aborting now")
2828 : end if
2829 :
2830 : ! broadening
2831 4 : if (brod.gt.0.009) then
2832 0 : write(std_out,*) '---------------------------------------------'
2833 0 : write(std_out,*) ' ATTENTION: broadening is quite high '
2834 0 : write(std_out,*) ' ideally should be less than 0.005 '
2835 0 : write(std_out,*) '---------------------------------------------'
2836 : else if (brod.gt.0.015) then
2837 : write(std_out,*) '----------------------------------------'
2838 : write(std_out,*) ' ATTENTION: broadening is too high '
2839 : write(std_out,*) ' ideally should be less than 0.005 '
2840 : write(std_out,*) '----------------------------------------'
2841 : end if
2842 :
2843 : ! tolerance
2844 4 : if (tol.gt.0.006) then
2845 0 : write(std_out,*) '----------------------------------------'
2846 0 : write(std_out,*) ' ATTENTION: tolerance is too high '
2847 0 : write(std_out,*) ' ideally should be less than 0.004 '
2848 0 : write(std_out,*) '----------------------------------------'
2849 : end if
2850 :
2851 : ! allocate local arrays
2852 4 : mband = ks_ebands%mband
2853 4 : ABI_CHECK(nband_sum <= mband, "nband_sum <= mband")
2854 12 : ABI_MALLOC(enk, (mband))
2855 20 : ABI_MALLOC(delta, (mband, mband, 3))
2856 20 : ABI_MALLOC(rmnbc, (mband, mband, 3, 3))
2857 12 : ABI_MALLOC(roverw, (mband, mband, 3, 3))
2858 12 : ABI_MALLOC(rmna, (mband, mband, 3))
2859 12 : ABI_MALLOC(chiw, (nmesh))
2860 8 : ABI_MALLOC(etaw, (nmesh))
2861 8 : ABI_MALLOC(chi2w, (nmesh))
2862 8 : ABI_MALLOC(eta2w, (nmesh))
2863 8 : ABI_MALLOC(sigmaw, (nmesh))
2864 8 : ABI_MALLOC(chi2, (nmesh))
2865 8 : ABI_MALLOC(eta1, (nmesh))
2866 20 : ABI_MALLOC(symrmn, (mband, mband, mband))
2867 :
2868 : ! generate the symmetrizing tensor
2869 4 : sym = zero
2870 100 : do isym=1,cryst%nsym
2871 1248 : s(:,:)=cryst%symrel_cart(:,:,isym)
2872 388 : do i=1,3
2873 1248 : do j=1,3
2874 3744 : do k=1,3
2875 3456 : sym(i,j,k)=sym(i,j,k)+(s(i,v1)*s(j,v2)*s(k,v3))
2876 : end do
2877 : end do
2878 : end do
2879 : end do
2880 :
2881 : ! initialise
2882 5056 : delta = zero
2883 15172 : rmnbc = zero
2884 604 : chiw = zero
2885 604 : chi2w = zero
2886 604 : chi2 = zero
2887 604 : etaw = zero
2888 604 : eta2w = zero
2889 604 : sigmaw = zero
2890 4 : my_emin=HUGE(one)
2891 4 : my_emax=-HUGE(one)
2892 :
2893 : ! Split work
2894 4 : call xmpi_split_work(ks_ebands%nkpt, comm, my_k1, my_k2)
2895 :
2896 : ! loop over kpts
2897 132 : do ik=my_k1,my_k2
2898 128 : write(std_out,*) "P-",my_rank,": ",ik,'of ', ks_ebands%nkpt
2899 260 : do isp=1,ks_ebands%nsppol
2900 : ! Calculate the scissor corrected energies and the energy window
2901 2688 : do ist1=1,nband_sum
2902 2560 : en = ks_ebands%eig(ist1,ik,isp)
2903 2560 : my_emin=min(my_emin,en)
2904 2560 : my_emax=max(my_emax,en)
2905 2560 : if(en > ks_ebands%fermie) then
2906 2048 : en = en + sc
2907 : end if
2908 2688 : enk(ist1) = en
2909 : end do
2910 :
2911 : ! calculate \Delta_nm and r_mn^a
2912 2688 : do istn=1,nband_sum
2913 2560 : en = enk(istn)
2914 53888 : do istm=1,nband_sum
2915 51200 : em = enk(istm)
2916 51200 : wmn = em - en
2917 204800 : delta(istn,istm,1:3)=pmat(istn,istn,ik,1:3,isp)-pmat(istm,istm,ik,1:3,isp)
2918 53760 : if(abs(wmn) < tol) then
2919 11776 : rmna(istm,istn,1:3) = zero
2920 : else
2921 193024 : rmna(istm,istn,1:3)=pmat(istm,istn,ik,1:3,isp)/wmn
2922 : end if
2923 : end do
2924 : end do
2925 :
2926 : ! calculate \r^b_mn;c
2927 2688 : do istm=1,nband_sum
2928 2560 : em = enk(istm)
2929 53888 : do istn=1,nband_sum
2930 51200 : en = enk(istn)
2931 51200 : wmn = em - en
2932 51200 : if (abs(wmn) < tol) then ! Degenerate energies
2933 38272 : rmnbc(istm,istn,:,:) = zero
2934 38272 : roverw(istm,istn,:,:) = zero
2935 : cycle
2936 : end if
2937 195584 : do ly = 1,3
2938 627328 : do lz = 1,3
2939 434304 : num1 = (rmna(istm,istn,ly)*delta(istm,istn,lz))+(rmna(istm,istn,lz)*delta(istm,istn,ly))
2940 434304 : den1 = wmn
2941 434304 : term1 = num1/den1
2942 434304 : term2 = zero
2943 9120384 : do istp=1,nband_sum
2944 8686080 : ep = enk(istp)
2945 8686080 : wmp = em - ep
2946 8686080 : wpn = ep - en
2947 8686080 : num2 = (wmp*rmna(istm,istp,ly)*rmna(istp,istn,lz))-(wpn*rmna(istm,istp,lz)*rmna(istp,istn,ly))
2948 8686080 : den2 = wmn
2949 9120384 : term2 = term2 + (num2/den2)
2950 : end do
2951 434304 : rmnbc(istm,istn,ly,lz) = -term1-(zi*term2)
2952 579072 : roverw(istm,istn,ly,lz) = (rmnbc(istm,istn,ly,lz)/wmn) - (rmna(istm,istn,ly)/(wmn**2))*delta(istm,istn,lz)
2953 : end do
2954 : end do
2955 : end do
2956 : end do
2957 :
2958 : ! initialise the factors
2959 : ! start the calculation
2960 2816 : do istn=1,nband_sum
2961 2560 : en=enk(istn)
2962 2560 : fn=ks_ebands%occ(istn,ik,isp)
2963 2560 : if(do_antiresonant .and. en .ge. ks_ebands%fermie) then
2964 : cycle
2965 : end if
2966 53888 : do istm=1,nband_sum
2967 51200 : em=enk(istm)
2968 51200 : if (do_antiresonant .and. em .le. ks_ebands%fermie) then
2969 : cycle
2970 : end if
2971 51200 : wmn=em-en
2972 51200 : wnm=-wmn
2973 51200 : fm = ks_ebands%occ(istm,ik,isp)
2974 51200 : fnm = fn - fm
2975 53760 : if(abs(wmn) > tol) then
2976 7286656 : chi1 = zero
2977 7286656 : chi2(:) = zero
2978 7286656 : chi2_1 = zero
2979 7286656 : chi2_2 = zero
2980 7286656 : eta1(:) = zero
2981 : eta1_1 = zero
2982 : eta1_2 = zero
2983 1013376 : eta2_1 = zero
2984 : eta2_2 = zero
2985 : sigma1 = zero
2986 1013376 : sigma2_1 = zero
2987 : ! Three band terms
2988 1013376 : do istl=1,nband_sum
2989 965120 : el=enk(istl)
2990 965120 : fl = ks_ebands%occ(istl,ik,isp)
2991 965120 : wlm = el-em
2992 965120 : wln = el-en
2993 965120 : wnl = -wln
2994 965120 : wml = em-el
2995 965120 : fnl = fn-fl
2996 965120 : fml = fm-fl
2997 965120 : flm = -fml
2998 965120 : fln = -fnl
2999 3860480 : do ly = 1,3
3000 12546560 : do lz = 1,3
3001 8686080 : symrmnl(ly,lz) = 0.5_dp*(rmna(istm,istn,ly)*rmna(istn,istl,lz)+rmna(istm,istn,lz)*rmna(istn,istl,ly))
3002 8686080 : symrlmn(ly,lz) = 0.5_dp*(rmna(istl,istm,ly)*rmna(istm,istn,lz)+rmna(istl,istm,lz)*rmna(istm,istn,ly))
3003 11581440 : symrmln(ly,lz) = 0.5_dp*(rmna(istm,istl,ly)*rmna(istl,istn,lz)+rmna(istm,istl,lz)*rmna(istl,istn,ly))
3004 : end do
3005 : end do
3006 :
3007 3908736 : do lx = 1,3
3008 12546560 : do ly = 1,3
3009 37639680 : do lz = 1,3
3010 26058240 : sigma1 = sigma1 + sym(lx,ly,lz)*(wnl*rmna(istl,istm,lx)*symrmnl(ly,lz)-wlm*rmna(istn,istl,lx)*symrlmn(ly,lz))
3011 26058240 : eta2_2 = eta2_2 + sym(lx,ly,lz)*fnm*rmna(istn,istm,lx)*symrmln(ly,lz)*(wml-wln)
3012 26058240 : if(abs(wln-wml) > tol) then
3013 26042688 : chi1 = chi1 + sym(lx,ly,lz)*(rmna(istn,istm,lx)*symrmln(ly,lz))/(wln-wml)
3014 : end if
3015 26058240 : eta1_1 = eta1_1 + sym(lx,ly,lz)*wln*rmna(istn,istl,lx)*symrlmn(ly,lz)
3016 26058240 : eta1_2 = eta1_2 - sym(lx,ly,lz)*wml*rmna(istl,istm,lx)*symrmnl(ly,lz)
3017 26058240 : if(abs(wnl-wmn) > tol) then
3018 26042688 : chi2_1 = chi2_1 - sym(lx,ly,lz)*(fnm*rmna(istl,istm,lx)*symrmnl(ly,lz)/(wnl-wmn))
3019 : end if
3020 34744320 : if(abs(wmn-wlm) > tol) then
3021 26042688 : chi2_2 = chi2_2 - sym(lx,ly,lz)*(fnm*rmna(istn,istl,lx)*symrlmn(ly,lz)/(wmn-wlm))
3022 : end if
3023 : end do
3024 : end do
3025 : end do
3026 : end do
3027 :
3028 : ! Two band terms
3029 : eta2_1 = zero
3030 : sigma2_1 = zero
3031 193024 : do lx = 1,3
3032 627328 : do ly = 1,3
3033 1881984 : do lz = 1,3
3034 : eta2_1 = eta2_1 + sym(lx,ly,lz)*fnm*rmna(istn,istm,lx)*0.5_dp &
3035 1302912 : *(delta(istm,istn,ly)*rmna(istm,istn,lz)+delta(istm,istn,lz)*rmna(istm,istn,ly))
3036 : ! Correct version (Sipe 1993)
3037 : sigma2_1 = sigma2_1 + sym(lx,ly,lz)*fnm*rmna(istn,istm,lx)*0.5_dp &
3038 1737216 : *(rmna(istm,istn,ly)*delta(istn,istm,lz)+rmna(istm,istn,lz)*delta(istn,istm,ly))
3039 :
3040 : ! Incorrect version (Hughes 1996)
3041 : !sigma2_1 = fnm*delta(istn,istm,v1)*0.5_dp*(rmna(istm,istn,v2)*rmna(istn,istm,v3)+rmna(istm,istn,v3)*rmna(istn,istm,v2))
3042 : end do
3043 : end do
3044 : end do
3045 :
3046 : ! calculate over the desired energy mesh and sum over k-points
3047 48256 : wtk = ks_ebands%wtk(ik)
3048 7286656 : do iw=1,nmesh
3049 7238400 : w=(iw-1)*de+idel
3050 7238400 : chi2w(iw) = chi2w(iw) + zi*wtk*((2.0_dp*fnm*chi1/(wmn-2.0_dp*w)))*const_esu ! Inter(2w) from chi
3051 7238400 : chiw(iw) = chiw(iw) + zi*wtk*((chi2_1+chi2_2)/(wmn-w))*const_esu ! Inter(w) from chi
3052 : eta2w(iw) = eta2w(iw) + zi*wtk*(8.0_dp*(eta2_1/((wmn**2)*(wmn-2.0_dp*w))) &
3053 7238400 : + 2.0_dp*eta2_2/((wmn**2)*(wmn-2.0_dp*w)))*const_esu ! Intra(2w) from eta
3054 7238400 : etaw(iw) = etaw(iw) + zi*wtk*((eta1_1 + eta1_2)*fnm/((wmn**2)*(wmn-w)))*const_esu ! Intra(w) from eta
3055 : sigmaw(iw) = sigmaw(iw) + 0.5_dp*zi*wtk*(fnm*sigma1/((wmn**2)*(wmn-w)) &
3056 7286656 : + (sigma2_1/((wmn**2)*(wmn-w))))*const_esu ! Intra(1w) from sigma
3057 : end do
3058 : end if
3059 : end do ! end loop over istn and istm
3060 : end do
3061 : end do ! spins
3062 : end do ! k-points
3063 :
3064 : ! Collect info among the nodes
3065 4 : call xmpi_min(my_emin,emin,comm,ierr)
3066 4 : call xmpi_max(my_emax,emax,comm,ierr)
3067 :
3068 4 : call xmpi_sum(chiw,comm,ierr)
3069 4 : call xmpi_sum(etaw,comm,ierr)
3070 4 : call xmpi_sum(chi2w,comm,ierr)
3071 4 : call xmpi_sum(eta2w,comm,ierr)
3072 4 : call xmpi_sum(sigmaw,comm,ierr)
3073 :
3074 : ! Master writes the output
3075 4 : if (my_rank == master) then
3076 :
3077 4 : if (ncid /= nctk_noid) then
3078 20 : start4 = [1, 1, icomp, itemp]
3079 20 : count4 = [2, nmesh, 1, 1]
3080 8 : ABI_MALLOC(chi2tot, (nmesh))
3081 608 : chi2tot = chiw + chi2w + etaw + eta2w + sigmaw
3082 4 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo2_chi2tot"), c2r(chi2tot), start=start4, count=count4))
3083 4 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo2_chiw"), c2r(chiw), start=start4, count=count4))
3084 4 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo2_etaw"), c2r(etaw), start=start4, count=count4))
3085 4 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo2_chi2w"), c2r(chi2w), start=start4, count=count4))
3086 4 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo2_eta2w"), c2r(eta2w), start=start4, count=count4))
3087 4 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo2_sigmaw"), c2r(sigmaw), start=start4, count=count4))
3088 4 : ABI_FREE(chi2tot)
3089 : end if
3090 :
3091 : ! write output in SI units and esu (esu to SI(m/v)=(value_esu)*(4xpi)/30000)
3092 4 : if (open_file(fnam1,msg,newunit=fout1,action='WRITE',form='FORMATTED') /= 0) then
3093 0 : ABI_ERROR(msg)
3094 : end if
3095 4 : if (open_file(fnam2,msg,newunit=fout2,action='WRITE',form='FORMATTED') /= 0) then
3096 0 : ABI_ERROR(msg)
3097 : end if
3098 4 : if (open_file(fnam3,msg,newunit=fout3,action='WRITE',form='FORMATTED') /= 0) then
3099 0 : ABI_ERROR(msg)
3100 : end if
3101 4 : if (open_file(fnam4,msg,newunit=fout4,action='WRITE',form='FORMATTED') /= 0) then
3102 0 : ABI_ERROR(msg)
3103 : end if
3104 4 : if (open_file(fnam5,msg,newunit=fout5,action='WRITE',form='FORMATTED') /= 0) then
3105 0 : ABI_ERROR(msg)
3106 : end if
3107 4 : if (open_file(fnam6,msg,newunit=fout6,action='WRITE',form='FORMATTED') /= 0) then
3108 0 : ABI_ERROR(msg)
3109 : end if
3110 4 : if (open_file(fnam7,msg,newunit=fout7,action='WRITE',form='FORMATTED') /= 0) then
3111 0 : ABI_ERROR(msg)
3112 : end if
3113 :
3114 : ! write headers
3115 4 : write(fout1, '(a,3i3)' ) ' #calculated the component:',v1,v2,v3
3116 4 : write(fout1, '(a,es16.6)' ) ' #tolerance:',tol
3117 4 : write(fout1, '(a,es16.6,a)' ) ' #broadening:',brod,'Ha'
3118 4 : write(fout1, '(a,es16.6,a)' ) ' #scissors shift:',sc,'Ha'
3119 4 : write(fout1, '(a,es16.6,a,es16.6,a)' ) ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
3120 4 : write(fout1, '(a)' )' # Energy Tot-Im Chi(-w,w,0) Tot-Im Chi(-w,w,0)'
3121 4 : write(fout1, '(a)' )' # eV *10^-7 esu *10^-12 m/V SI units '
3122 4 : write(fout1, '(a)' )' # '
3123 :
3124 4 : write(fout2, '(a,3i3)' ) ' #calculated the component:',v1,v2,v3
3125 4 : write(fout2, '(a,es16.6)') ' #tolerance:',tol
3126 4 : write(fout2, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
3127 4 : write(fout2, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
3128 4 : write(fout2, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
3129 4 : write(fout2, '(a)')' # Energy Tot-Re Chi(-w,w,0) Tot-Re Chi(-w,w,0)'
3130 4 : write(fout2, '(a)')' # eV *10^-7 esu *10^-12 m/V SI units '
3131 4 : write(fout2, '(a)')' # '
3132 :
3133 4 : write(fout3, '(a,3i3)') ' #calculated the component:',v1,v2,v3
3134 4 : write(fout3, '(a,es16.6)') ' #tolerance:',tol
3135 4 : write(fout3, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
3136 4 : write(fout3, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
3137 4 : write(fout3, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
3138 4 : write(fout3, '(a)')' # Energy(eV) Inter(2w) inter(1w) intra(2w) intra(1w)'
3139 4 : write(fout3, '(a)')' # in esu'
3140 4 : write(fout3, '(a)')' # '
3141 :
3142 4 : write(fout4, '(a,3i3)') ' #calculated the component:',v1,v2,v3
3143 4 : write(fout4, '(a,es16.6)') ' #tolerance:',tol
3144 4 : write(fout4, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
3145 4 : write(fout4, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
3146 4 : write(fout4, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
3147 4 : write(fout4, '(a)')' # Energy(eV) Inter(2w) inter(1w) intra(2w) intra(1w)'
3148 4 : write(fout4, '(a)')' # in esu'
3149 4 : write(fout4, '(a)')' # '
3150 :
3151 4 : write(fout5, '(a,3i3)') ' #calculated the component:',v1,v2,v3
3152 4 : write(fout5, '(a,es16.6)') ' #tolerance:',tol
3153 4 : write(fout5, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
3154 4 : write(fout5, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
3155 4 : write(fout5, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
3156 4 : write(fout5, '(a)')' # Energy(eV) |TotChi(-w,w,0)| |Tot Chi(-w,w,0)|'
3157 4 : write(fout5, '(a)')' # eV *10^-7 esu *10^-12 m/V SI units '
3158 4 : write(fout5, '(a)')' # '
3159 :
3160 4 : write(fout6, '(a,3i3)') ' #calculated the component:',v1,v2,v3
3161 4 : write(fout6, '(a,es16.6)') ' #tolerance:',tol
3162 4 : write(fout6, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
3163 4 : write(fout6, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
3164 4 : write(fout6, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
3165 4 : write(fout6, '(a)')' # Energy(eV) Chi(w) Eta(w) Sigma(w)'
3166 4 : write(fout6, '(a)')' # in esu'
3167 4 : write(fout6, '(a)')' # '
3168 :
3169 4 : write(fout7, '(a,3i3)') ' #calculated the component:',v1,v2,v3
3170 4 : write(fout7, '(a,es16.6)') ' #tolerance:',tol
3171 4 : write(fout7, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
3172 4 : write(fout7, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
3173 4 : write(fout7, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Ha'
3174 4 : write(fout7, '(a)')' # Energy(eV) Chi(w) Eta(w) Sigma(w)'
3175 4 : write(fout7, '(a)')' # in esu'
3176 4 : write(fout7, '(a)')' # '
3177 :
3178 4 : totim=zero
3179 4 : totre=zero
3180 4 : totabs=zero
3181 600 : do iw=2,nmesh
3182 596 : ene=(iw-1)*de
3183 596 : ene=ene*Ha_eV
3184 :
3185 596 : totim=aimag(chiw(iw)+chi2w(iw)+etaw(iw)+eta2w(iw)+sigmaw(iw))/1.d-7
3186 596 : write(fout1,'(f15.6,2es15.6)') ene,totim,totim*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
3187 596 : totim=zero
3188 :
3189 596 : totre=dble(chiw(iw)+chi2w(iw)+eta2w(iw)+etaw(iw)+sigmaw(iw))/1.d-7
3190 596 : write(fout2,'(f15.6,2es15.6)') ene,totre,totre*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
3191 596 : totre=zero
3192 :
3193 596 : write(fout3,'(f15.6,4es15.6)') ene,aimag(chi2w(iw))/1.d-7,aimag(chiw(iw))/1.d-7, &
3194 1192 : aimag(eta2w(iw))/1.d-7,aimag(etaw(iw))/1.d-7+aimag(sigmaw(iw))/1.d-7
3195 :
3196 596 : write(fout4,'(f15.6,4es15.6)') ene,dble(chi2w(iw))/1.d-7,aimag(chiw(iw))/1.d-7, &
3197 1192 : dble(eta2w(iw))/1.d-7,dble(etaw(iw))/1.d-7+dble(sigmaw(iw))/1.d-7
3198 :
3199 596 : totabs=abs(chiw(iw)+chi2w(iw)+etaw(iw)+eta2w(iw)+sigmaw(iw))/1.d-7
3200 596 : write(fout5,'(f15.6,2es15.6)') ene,totabs,totabs*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
3201 596 : totabs=zero
3202 :
3203 596 : write(fout6,'(f15.6,4es15.6)') ene,aimag(chi2w(iw)+chiw(iw))/1.d-7, &
3204 1192 : aimag(eta2w(iw)+etaw(iw))/1.d-7,aimag(sigmaw(iw))/1.d-7
3205 :
3206 596 : write(fout7,'(f15.6,4es15.6)') ene,dble(chi2w(iw)+chiw(iw))/1.d-7, &
3207 1196 : dble(eta2w(iw)+etaw(iw))/1.d-7,dble(sigmaw(iw))/1.d-7
3208 : end do
3209 :
3210 4 : close(fout1)
3211 4 : close(fout2)
3212 4 : close(fout3)
3213 4 : close(fout4)
3214 4 : close(fout5)
3215 4 : close(fout6)
3216 4 : close(fout7)
3217 :
3218 : ! print information
3219 4 : write(std_out,*) ' '
3220 4 : write(std_out,*) 'information about calculation just performed:'
3221 4 : write(std_out,*) ' '
3222 4 : write(std_out,*) 'calculated the component:',v1,v2,v3 ,'of the nonlinear electro-optical susceptibility'
3223 4 : write(std_out,*) 'tolerance:',tol
3224 4 : if (tol.gt.0.008) write(std_out,*) 'ATTENTION: tolerance is too high'
3225 4 : write(std_out,*) 'broadening:',brod,'Hartree'
3226 4 : if (brod.gt.0.009) then
3227 0 : write(std_out,*) ' '
3228 0 : write(std_out,*) 'ATTENTION: broadening is quite high'
3229 0 : write(std_out,*) ' '
3230 : else if (brod.gt.0.015) then
3231 : write(std_out,*) ' '
3232 : write(std_out,*) 'ATTENTION: broadening is too high'
3233 : write(std_out,*) ' '
3234 : end if
3235 4 : write(std_out,*) 'scissors shift:',sc,'Hartree'
3236 4 : write(std_out,*) 'energy window:',(emax-emin)*Ha_eV,'eV',(emax-emin),'Hartree'
3237 :
3238 : end if
3239 :
3240 : ! deallocate local arrays
3241 4 : ABI_FREE(enk)
3242 4 : ABI_FREE(delta)
3243 4 : ABI_FREE(rmnbc)
3244 4 : ABI_FREE(roverw)
3245 4 : ABI_FREE(rmna)
3246 4 : ABI_FREE(chiw)
3247 4 : ABI_FREE(chi2w)
3248 4 : ABI_FREE(chi2)
3249 4 : ABI_FREE(etaw)
3250 4 : ABI_FREE(eta1)
3251 4 : ABI_FREE(symrmn)
3252 4 : ABI_FREE(eta2w)
3253 4 : ABI_FREE(sigmaw)
3254 :
3255 12 : end subroutine nonlinopt
3256 : !!***
3257 :
3258 : !----------------------------------------------------------------------
3259 :
3260 : end module m_optic_tools
|
