Branch data 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
|
