Line data Source code
1 : !!****m* ABINIT/m_screening
2 : !! NAME
3 : !! m_screening
4 : !!
5 : !! FUNCTION
6 : !! This module contains the definition of the object used to deal
7 : !! with the inverse dielectric matrix as well as related methods.
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2008-2026 ABINIT group (MG)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : module m_screening
24 :
25 : use, intrinsic :: iso_c_binding
26 : use defs_basis
27 : use m_abicore
28 : use m_hide_blas
29 : use m_linalg_interfaces
30 : use m_xmpi
31 : use m_errors
32 : use m_copy
33 : use m_splines
34 : use m_lebedev
35 : use m_spectra
36 : use m_nctk
37 : use m_distribfft
38 : use netcdf
39 :
40 : use defs_abitypes, only : MPI_type
41 : use m_gwdefs, only : GW_TOLQ0, czero_gw, GW_Q0_DEFAULT
42 : use m_fstrings, only : toupper, endswith, sjoin, itoa, strcat
43 : use m_io_tools, only : open_file
44 : use m_time, only : cwtime, cwtime_report
45 : use m_numeric_tools, only : print_arr, hermitianize
46 : use m_special_funcs, only : k_fermi, k_thfermi
47 : use m_geometry, only : normv, vdotw, metric
48 : use m_hide_lapack, only : xginv
49 : use m_crystal, only : crystal_t
50 : use m_bz_mesh, only : kmesh_t, box_len
51 : use m_fft_mesh, only : g2ifft
52 : use m_fftcore, only : kgindex
53 : use m_fft, only : fourdp
54 : use m_gsphere, only : gsphere_t
55 : use m_vcoul, only : vcoul_t
56 : use m_io_screening, only : read_screening, write_screening, &
57 : HSCR_LATEST_HEADFORM, hscr_t, ncname_from_id, em1_ncname
58 : use m_paw_sphharm, only : ylmc
59 : use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq
60 : use m_pstat, only : pstat_proc
61 :
62 : implicit none
63 :
64 : private
65 : !!***
66 :
67 : !----------------------------------------------------------------------
68 :
69 : !!****t* m_screening/epsm1_t
70 : !! NAME
71 : !! epsm1_t
72 : !!
73 : !! FUNCTION
74 : !! For the GW part of ABINIT, the epsm1_t structured datatype
75 : !! gather the results of screening: the inverse dielectric matrix, and the omega matrices.
76 : !!
77 : !! SOURCE
78 :
79 : type,public :: epsm1_t
80 :
81 : integer :: id
82 : ! Matrix identifier: O if not yet defined, 1 for chi0,
83 : ! 2 for chi, 3 for epsilon, 4 for espilon^{-1}, 5 for W.
84 :
85 : integer :: ikxc
86 : ! Kxc kernel used, 0 for None (RPA), >0 for static TDDFT (=ixc), <0 for TDDFT
87 :
88 : integer :: fform
89 : ! File format: 1002 for SCR|SUSC files.
90 :
91 : integer :: mqmem
92 : ! =0 for out-of-core solution, =nqibz if entire matrix is stored in memory.
93 :
94 : integer :: nI,nJ
95 : ! Number of components (rows,columns) in chi|eps^-1. (1,1) if collinear.
96 :
97 : integer :: nqibz
98 : ! Number of q-points in the IBZ used.
99 :
100 : integer :: nqlwl
101 : ! Number of point used for the treatment of the long wave-length limit.
102 :
103 : integer :: nomega
104 : ! Total number of frequencies.
105 :
106 : integer :: nomega_i
107 : ! Number of purely imaginary frequencies used.
108 :
109 : integer :: nomega_i_conv
110 : ! Number of purely imaginary frequencies used for convolution.
111 :
112 : integer :: nomega_r
113 : ! Number of real frequencies used.
114 :
115 : integer :: npwe
116 : ! Number of G vectors.
117 :
118 : integer :: test_type
119 : ! 0 for None, 1 for TEST-PARTICLE, 2 for TEST-ELECTRON (only for TDDFT)
120 :
121 : integer :: tordering
122 : ! 0 if not defined, 1 for Time-Ordered, 2 for Advanced, 3 for Retarded.
123 :
124 : character(len=fnlen) :: fname
125 : ! Name of the file from which epsm1 is read.
126 :
127 : integer,allocatable :: gvec(:,:)
128 : ! gvec(3,npwe)
129 : ! G-vectors used to describe the two-point function (r.l.u.).
130 :
131 : real(dp),allocatable :: qibz(:,:)
132 : ! qibz(3,nqibz)
133 : ! q-points in reduced coordinates
134 :
135 : real(dp),allocatable :: qlwl(:,:)
136 : ! qlwl(3,nqlwl)
137 : ! q-points used for the long wave-length limit treatment.
138 :
139 : logical :: use_mpi_shared_win = .False.
140 : ! This flag allows us to understand if espm1 is associates to a MPI window or not.
141 : ! In the former case, we should free the window and not the pointer
142 :
143 : integer :: epsm1_win = xmpi_undefined
144 :
145 : type(xcomm_t) :: shared_comm
146 :
147 : complex(gwp), contiguous, pointer :: epsm1(:,:,:,:) => null()
148 : ! epsm1(npwe,npwe,nomega,nqibz)
149 : ! Contains the two-point function $\epsilon_{G,Gp}(q,omega)$ in frequency and reciprocal space.
150 : ! We use a pointer so that we can associated it to MPI shared memory window.
151 :
152 : complex(gwp), contiguous, pointer :: epsm1_qbz(:,:,:) => null()
153 : integer :: epsm1_qbz_win = xmpi_undefined
154 :
155 : complex(dp),allocatable :: omega(:)
156 : ! omega(nomega)
157 : ! Frequencies used both along the real and the imaginary axis.
158 :
159 : type(hscr_t) :: Hscr
160 : ! The header reported in the _SCR of _SUSC file.
161 : ! This object contains information on the susceptibility or the inverse dielectric matrix
162 : ! as stored in the external file. These quantities do *NOT* correspond to the quantities
163 : ! used during the GW calculation since some parameters might differ, actually they might be smaller.
164 : ! For example, the number of G-vectors used can be smaller than the number of G"s stored on file.
165 :
166 : contains
167 :
168 : procedure :: free => epsm1_free
169 : ! Free memory
170 :
171 : procedure :: print => epsm1_print
172 : ! Print basic info
173 :
174 : procedure :: rotate_iqbz => epsm1_rotate_iqbz
175 : ! Symmetrize two-point function at a q-point in the BZ.
176 :
177 : procedure :: rotate_iqbz_inplace => epsm1_rotate_iqbz_inplace
178 : ! In-place version of the above
179 :
180 : procedure :: from_file => epsm1_from_file
181 : ! Initialize the object from file
182 :
183 : procedure :: mkdump => epsm1_mkdump
184 : ! Dump the object to a file.
185 :
186 : procedure :: get_epsm1 => get_epsm1
187 :
188 : procedure :: decompose_epsm1 => decompose_epsm1
189 : ! Decompose the complex symmetrized dielectric
190 :
191 : procedure :: malloc_epsm1_qbz => epsm1_malloc_epsm1_qbz
192 :
193 : procedure :: free_epsm1_qbz => epsm1_free_epsm1_qbz
194 :
195 : end type epsm1_t
196 :
197 : public :: make_epsm1_driver ! Calculate the inverse symmetrical dielectric matrix starting from chi0
198 : public :: mkem1_q0 ! Construct the microscopic dielectric matrix for q-->0
199 :
200 : ! Routines for the model dielectric function
201 :
202 : public :: screen_mdielf ! Calculates W_{G,G'}(q,w) for a given q-point in the BZ using a model dielectric function.
203 : public :: rpa_symepsm1
204 : !!***
205 :
206 : !!****t* m_screening/chi_t
207 : !! NAME
208 : !! chi_t
209 : !!
210 : !! FUNCTION
211 : !! This object contains the head and the wings of the polarizability
212 : !! These quantities are used to treat the q-->0 limit
213 : !!
214 : !! SOURCE
215 :
216 : type,public :: chi_t
217 :
218 : integer :: npwe
219 : ! Number of G vectors.
220 :
221 : integer :: nomega
222 : ! Number of frequencies
223 :
224 : complex(gwp),allocatable :: mat(:,:,:)
225 : ! mat(npwe, npwe, nomega)
226 :
227 : complex(dp),allocatable :: head(:,:,:)
228 : ! head(3,3,nomega)
229 :
230 : complex(dp),allocatable :: lwing(:,:,:)
231 : ! lwing(3,npwe,nomega)
232 : ! Lower wings
233 :
234 : complex(dp),allocatable :: uwing(:,:,:)
235 : ! uwing(3,npwe,nomega)
236 : ! Upper wings.
237 :
238 : end type chi_t
239 :
240 : public :: chi_new ! Create new object (allocate memory)
241 : public :: chi_free ! Free memory.
242 : !!***
243 :
244 : !!****t* m_screening/lwl_t
245 : !! NAME
246 : !! lwl_t
247 : !!
248 : !! FUNCTION
249 : !!
250 : !! SOURCE
251 :
252 : type,public :: lwl_t
253 :
254 : integer :: npwe
255 : ! Number of G vectors.
256 :
257 : integer :: nomega
258 : ! Number of frequencies
259 :
260 : integer :: method
261 : ! 1 = Only head
262 : ! 2 = head + wings
263 : ! 3 = head + wings + body corrections.
264 :
265 : character(len=fnlen) :: fname
266 : ! Name of the file from which epsm1 is read.
267 :
268 : complex(dp),allocatable :: head(:,:,:)
269 : ! head(3,3,nomega)
270 :
271 : complex(dp),allocatable :: lwing(:,:,:)
272 : ! lwing(3,npwe,nomega)
273 : ! Lower wings
274 :
275 : complex(dp),allocatable :: uwing(:,:,:)
276 : ! uwing(3,npwe,nomega)
277 : ! Upper wings.
278 :
279 : complex(dp),allocatable :: body(:,:,:)
280 : ! uwing(npwe,npwe,nomega)
281 : ! Body terms
282 :
283 : end type lwl_t
284 :
285 : public :: lwl_write
286 : public :: lwl_init
287 : !public :: lwl_from_file
288 : public :: lwl_free
289 : !!***
290 :
291 : #define _MOK(integer) int(integer, kind=XMPI_OFFSET_KIND)
292 :
293 : CONTAINS !========================================================================================
294 : !!***
295 :
296 : !----------------------------------------------------------------------
297 :
298 : !!****f* m_screening/epsm1_free
299 : !! NAME
300 : !! epsm1_free
301 : !!
302 : !! FUNCTION
303 : !! Deallocate all the pointers in epsm1 that result to be associated.
304 : !! Perform also a cleaning of the Header.
305 : !!
306 : !! SOURCE
307 :
308 202 : subroutine epsm1_free(epsm1)
309 :
310 : !Arguments ------------------------------------
311 : class(epsm1_t),intent(inout) :: epsm1
312 :
313 : integer :: ierr
314 : ! *************************************************************************
315 :
316 : !integer
317 202 : ABI_SFREE(epsm1%gvec)
318 :
319 : !real
320 202 : ABI_SFREE(epsm1%qibz)
321 202 : ABI_SFREE(epsm1%qlwl)
322 :
323 : !complex
324 202 : if (epsm1%use_mpi_shared_win) then
325 30 : call xmpi_win_free(epsm1%epsm1_win, ierr)
326 30 : nullify(epsm1%epsm1)
327 30 : if (epsm1%epsm1_qbz_win /= xmpi_undefined) then
328 0 : call xmpi_win_free(epsm1%epsm1_qbz_win, ierr)
329 0 : nullify(epsm1%epsm1_qbz)
330 : end if
331 : else
332 172 : ABI_SFREE_PTR(epsm1%epsm1)
333 172 : ABI_SFREE_PTR(epsm1%epsm1_qbz)
334 : end if
335 202 : call epsm1%shared_comm%free()
336 :
337 202 : ABI_SFREE(epsm1%omega)
338 :
339 : !datatypes
340 202 : call epsm1%Hscr%free()
341 :
342 202 : end subroutine epsm1_free
343 : !!***
344 :
345 : !----------------------------------------------------------------------
346 :
347 : !!****f* m_screening/epsm1_print
348 : !! NAME
349 : !! epsm1_print
350 : !!
351 : !! FUNCTION
352 : !! Print the basic dimensions and the most important
353 : !! quantities reported in the epsm1_t data type.
354 : !!
355 : !! INPUTS
356 : !! epsm1<epsm1_t>=The data type.
357 : !! unit[optional]=the unit number for output.
358 : !! prtvol[optional]=verbosity level.
359 : !!
360 : !! OUTPUT
361 : !! Only printing.
362 : !!
363 : !! SOURCE
364 :
365 5 : subroutine epsm1_print(epsm1, units, prtvol)
366 :
367 : !Arguments ------------------------------------
368 : class(epsm1_t),intent(in) :: epsm1
369 : integer,intent(in) :: units(:)
370 : integer,intent(in) :: prtvol
371 :
372 : !Local variables-------------------------------
373 : integer :: iw,iq_ibz,iqlwl
374 : character(len=50) :: rfname,rforder,rfapprox,rftest,kxcname
375 : character(len=500) :: msg
376 : ! *************************************************************************
377 :
378 : ! chi0 or \epsilon^{-1} ?
379 5 : SELECT CASE (epsm1%ID)
380 : CASE (0)
381 0 : rfname = 'Undefined'
382 : CASE (1)
383 0 : rfname = 'Irreducible Polarizability'
384 : CASE (2)
385 0 : rfname = 'Polarizability'
386 : CASE (3)
387 0 : rfname = 'Symmetrical Dielectric Matrix'
388 : CASE (4)
389 5 : rfname = 'Symmetrical Inverse Dielectric Matrix'
390 : CASE DEFAULT
391 5 : ABI_BUG(sjoin('Wrong epsm1%ID:',itoa(epsm1%ID)))
392 : END SELECT
393 :
394 : ! For chi, \espilon or \epsilon^{-1}, define the approximation.
395 : rfapprox='None'
396 : if (epsm1%ID>=2.or.epsm1%ID<=4) then
397 5 : if (epsm1%ikxc==0) then
398 5 : rfapprox='RPA'
399 0 : else if (epsm1%ikxc>0) then
400 0 : rfapprox='Static TDDFT'
401 : else
402 0 : rfapprox='TDDFT'
403 : end if
404 : end if
405 :
406 : ! === If TDDFT and \epsilon^{-1}, define the type ===
407 5 : rftest='None'
408 : ! if (epsm1%ID==0) then
409 : ! if (epsm1%test_type==0) then
410 : ! rftest='TEST-PARTICLE'
411 : ! else if (epsm1%test_type==1) then
412 : ! rftest='TEST-ELECTRON'
413 : ! else
414 : ! write(msg,'(4a,i3)')ch10,&
415 : !& ' epsm1_print : BUG - ',ch10,&
416 : !& ' Wrong value of epsm1%test_type = ',epsm1%test_type
417 : ! ABI_ERROR(msg)
418 : ! end if
419 : ! end if
420 :
421 : ! === Define time-ordering ===
422 5 : rforder='Undefined'
423 5 : if (epsm1%Tordering==1) then
424 5 : rforder='Time-Ordered'
425 0 : else if (epsm1%Tordering==2) then
426 0 : rforder='Advanced'
427 0 : else if (epsm1%Tordering==3) then
428 0 : rforder='Retarded'
429 : else
430 0 : ABI_BUG(sjoin('Wrong epsm1%tordering= ',itoa(epsm1%Tordering)))
431 : end if
432 :
433 5 : kxcname='None'
434 5 : if (epsm1%ikxc/=0) then
435 : !TODO Add function to retrieve kxc name
436 0 : ABI_ERROR('Add function to retrieve kxc name')
437 0 : kxcname='XXXXX'
438 : end if
439 :
440 5 : write(msg,'(6a,5(3a))')ch10,&
441 5 : ' ==== Info on the Response Function ==== ',ch10,&
442 5 : ' Associated File ................ ',TRIM(epsm1%fname),ch10,&
443 5 : ' Response Function Type .......... ',TRIM(rfname),ch10,&
444 5 : ' Type of Approximation ........... ',TRIM(rfapprox),ch10,&
445 5 : ' XC kernel used .................. ',TRIM(kxcname),ch10,&
446 5 : ' Type of probing particle ........ ',TRIM(rftest),ch10,&
447 10 : ' Time-Ordering ................... ',TRIM(rforder),ch10
448 5 : call wrtout(units, msg)
449 : write(msg,'(a,2i4,a,3(a,i4,a),a,3i4,2a,i4,a)')&
450 5 : ' Number of components ............ ',epsm1%nI,epsm1%nJ,ch10,&
451 5 : ' Number of q-points in the IBZ ... ',epsm1%nqibz,ch10,&
452 5 : ' Number of q-points for q-->0 .... ',epsm1%nqlwl,ch10,&
453 5 : ' Number of G-vectors ............. ',epsm1%npwe,ch10,&
454 5 : ' Number of frequencies ........... ',epsm1%nomega,epsm1%nomega_r,epsm1%nomega_i,ch10,&
455 10 : ' Value of mqmem .................. ',epsm1%mqmem,ch10
456 5 : call wrtout(units,msg)
457 :
458 5 : if (epsm1%nqlwl/=0) then
459 5 : write(msg,'(a,i3)')' q-points for long wavelength limit: ',epsm1%nqlwl
460 5 : call wrtout(units, msg)
461 10 : do iqlwl=1,epsm1%nqlwl
462 5 : write(msg,'(1x,i5,a,3es16.8)')iqlwl,') ',epsm1%qlwl(:,iqlwl)
463 10 : call wrtout(units, msg)
464 : end do
465 : end if
466 :
467 5 : if (prtvol > 0) then
468 : ! Print out head and wings in the long-wavelength limit.
469 : ! TODO add additional stuff.
470 0 : write(msg,'(a,i4)')' Calculated Frequencies: ',epsm1%nomega
471 0 : call wrtout(units, msg)
472 0 : do iw=1,epsm1%nomega
473 0 : write(msg,'(i4,es14.6)')iw,epsm1%omega(iw)*Ha_eV
474 0 : call wrtout(units, msg)
475 : end do
476 :
477 0 : write(msg,'(a,i4)')' Calculated q-points: ',epsm1%nqibz
478 0 : call wrtout(units, msg)
479 0 : do iq_ibz=1,epsm1%nqibz
480 0 : write(msg,'(1x,i4,a,3es16.8)')iq_ibz,') ',epsm1%qibz(:,iq_ibz)
481 0 : call wrtout(units, msg)
482 : end do
483 : end if ! prtvol > 0
484 :
485 5 : end subroutine epsm1_print
486 : !!***
487 :
488 : !----------------------------------------------------------------------
489 :
490 : !!****f* m_screening/Epsm1_rotate_iqbz
491 : !! NAME
492 : !! Epsm1_rotate_iqbz
493 : !!
494 : !! FUNCTION
495 : !! Symmetrize the inverse dielectric matrix, namely calculate epsilon^{-1} at a generic
496 : !! q-point in the BZ starting from the knowledge of the matrix at a q-point in the IBZ.
497 : !! The procedure is quite generic and can be used for every two-point function which has
498 : !! the same symmetry as the crystal.
499 : !!
500 : !! INPUTS
501 : !! iq_bz=Index of the q-point in the BZ where epsilon^-1 is required.
502 : !! nomega=Number of frequencies required. All frequencies from 1 up to nomega are symmetrized.
503 : !! npwc=Number of G vectors in symmetrized matrix, has to be smaller than epsm1%npwe.
504 : !! Gsph<gsphere_t>=data related to the G-sphere
505 : !! Qmesh<kmesh_t>=Structure defining the q-mesh used for epsm1.
506 : !! remove_exchange=If .TRUE., return e^{-1}-1 namely remove the exchange part.
507 :
508 : !! OUTPUT
509 : !! epsm1_qbz(npwc,npwc,nomega)=The inverse dielectric matrix at the q-point defined by iq_bz.
510 : !! Exchange part can be subtracted out.
511 : !!
512 : !! NOTES
513 : !! In the present implementation we are not considering a possible umklapp vector G0 in the
514 : !! expression Sq = q+G0. Treating this case would require some changes in the G-sphere
515 : !! since we have to consider G-G0. The code however stops in sigma if a nonzero G0 is required
516 : !! to reconstruct the BZ.
517 : !!
518 : !! * Remember the symmetry properties of \tilde\espilon^{-1}
519 : !!
520 : !! If q_bz = S q_ibz + G0:
521 : !!
522 : !! $\epsilon^{-1}_{SG1-G0, SG2-G0}(q_bz) = e^{+iS(G2-G1).\tau} \epsilon^{-1}_{G1, G2)}(q)
523 : !!
524 : !! If time-reversal symmetry can be used then:
525 : !!
526 : !! $\epsilon^{-1}_{G1,G2}(-q_bz) = e^{+i(G1-G2).\tau} \epsilon^{-1}_{-S^{-1}(G1+Go), -S^{-1}(G2+G0)}^*(q)
527 : !!
528 : !! TODO
529 : !! Symmetrization can be skipped if iq_bz correspond to a point in the IBZ
530 : !!
531 : !! SOURCE
532 :
533 5981 : subroutine Epsm1_rotate_iqbz(epsm1, iq_bz, nomega, npwc, Gsph, Qmesh, remove_exchange)
534 :
535 : !Arguments ------------------------------------
536 : !scalars
537 : class(epsm1_t),intent(in) :: epsm1
538 : integer,intent(in) :: iq_bz, nomega, npwc
539 : logical,intent(in) :: remove_exchange
540 : type(gsphere_t),target,intent(in) :: Gsph
541 : type(kmesh_t),intent(in) :: Qmesh
542 :
543 : !Local variables-------------------------------
544 : !scalars
545 : integer :: iw,ii,jj,iq_ibz,itim_q,isym_q,iq_loc,sg1,sg2, ierr, g0(3)
546 : complex(gwp) :: phmsg1t,phmsg2t_star
547 : !arrays
548 : real(dp) :: qbz(3)
549 : ! *********************************************************************
550 :
551 5981 : ABI_CHECK(epsm1%nomega >= nomega, 'Too many frequencies required')
552 5981 : ABI_CHECK(epsm1%npwe >= npwc, 'Too many G-vectors required')
553 :
554 : ! Get iq_ibz, and symmetries from iq_ibz.
555 5981 : call Qmesh%get_BZ_item(iq_bz, qbz, iq_ibz, isym_q, itim_q, umklp=g0)
556 23924 : ABI_CHECK(all(g0 == 0), "non-zero g0 is not coded")
557 :
558 : ! If out-of-memory, only epsm1%espm1(:,:,:,1) has been allocated and filled.
559 5981 : iq_loc = iq_ibz; if (epsm1%mqmem == 0) iq_loc=1
560 :
561 5981 : if (epsm1%use_mpi_shared_win) call xmpi_win_fence(XMPI_MODE_NOPRECEDE, epsm1%epsm1_qbz_win, ierr) ! Start the RMA epoch.
562 :
563 : ! MG: rottb is a 1-1 mapping, hence we can collapse the loops (false sharing is not an issue here).
564 : !!$OMP PARALLEL DO COLLAPSE(2) PRIVATE(sg2,sg1,phmsg1t,phmsg2t_star) IF (NOT epsm1%use_mpi_shared_win)
565 118367 : do iw=1,nomega
566 112386 : if (epsm1%shared_comm%skip(iw)) cycle ! MPI parallelism with shared memory.
567 2831841 : do jj=1,npwc
568 2738050 : sg2 = Gsph%rottb(jj, itim_q, isym_q)
569 2738050 : phmsg2t_star = CONJG(Gsph%phmSGt(jj, isym_q))
570 119899638 : do ii=1,npwc
571 117073778 : sg1 = Gsph%rottb(ii,itim_q,isym_q)
572 117073778 : phmsg1t = Gsph%phmSGt(ii,isym_q)
573 119811828 : epsm1%epsm1_qbz(sg1,sg2,iw) = epsm1%epsm1(ii,jj,iw,iq_loc) * phmsg1t * phmsg2t_star
574 : end do
575 : end do
576 : end do
577 :
578 : ! Account for time-reversal
579 5981 : if (itim_q==2) then
580 : !!$OMP PARALLEL DO IF (nomega > 1)
581 0 : do iw=1,nomega
582 0 : if (epsm1%shared_comm%skip(iw)) cycle ! MPI parallelism with shared memory.
583 0 : call sqmat_itranspose(npwc, epsm1%epsm1_qbz(:,:,iw))
584 : end do
585 : end if
586 :
587 5981 : if (remove_exchange) then
588 : ! Subtract the exchange contribution
589 : ! If it's a pole screening, the exchange contribution is already removed
590 : !!$OMP PARALLEL DO IF (nomega > 1)
591 118367 : do iw=1,nomega
592 112386 : if (epsm1%shared_comm%skip(iw)) cycle ! MPI parallelism with shared memory.
593 2831841 : do ii=1,npwc
594 2825860 : epsm1%epsm1_qbz(ii,ii,iw) = epsm1%epsm1_qbz(ii,ii,iw) - cmplx(1.0_gwp,0.0_gwp)
595 : end do
596 : end do
597 : endif
598 :
599 5981 : if (epsm1%use_mpi_shared_win) call xmpi_win_fence(XMPI_MODE_NOSUCCEED, epsm1%epsm1_qbz_win, ierr) ! Close the RMA epoch.
600 :
601 5981 : end subroutine Epsm1_rotate_iqbz
602 : !!***
603 :
604 : !----------------------------------------------------------------------
605 :
606 : !!****f* m_screening/Epsm1_rotate_iqbz_inplace
607 : !! NAME
608 : !! Epsm1_rotate_iqbz_inplace
609 : !!
610 : !! FUNCTION
611 : !! Same function as Epsm1_rotate_iqbz, but now the array Ep%epsm1 is modified inplace
612 : !! via an auxiliary work array of shape (npwc,npwc)
613 : !!
614 : !! INPUTS
615 : !! iq_bz=Index of the q-point in the BZ where epsilon^-1 is required.
616 : !! nomega=Number of frequencies required. All frequencies from 1 up to nomega are symmetrized.
617 : !! npwc=Number of G vectors in symmetrized matrix, has to be smaller than epsm1%npwe.
618 : !! Gsph<gsphere_t>=data related to the G-sphere
619 : !! Qmesh<kmesh_t>=Structure defining the q-mesh used for epsm1.
620 : !! remove_exchange=If .TRUE., return e^{-1}-1 namely remove the exchange part.
621 : !!
622 : !! OUTPUT
623 : !! epsm1%epsm1(npwc,npwc,nomega,iq_loc) symmetrised
624 : !!
625 : !! NOTES
626 : !! In the present implementation we are not considering a possible umklapp vector G0 in the
627 : !! expression Sq = q+G0. Treating this case would require some changes in the G-sphere
628 : !! since we have to consider G-G0. The code however stops in sigma if a nonzero G0 is required
629 : !! to reconstruct the BZ.
630 : !!
631 : !! * Remember the symmetry properties of \tilde\espilon^{-1}
632 : !! If q_bz=Sq_ibz+G0:
633 : !!
634 : !! $\epsilon^{-1}_{SG1-G0,SG2-G0}(q_bz) = e^{+iS(G2-G1).\tau}\epsilon^{-1}_{G1,G2)}(q)
635 : !!
636 : !! If time-reversal symmetry can be used then :
637 : !! $\epsilon^{-1}_{G1,G2}(-q_bz) = e^{+i(G1-G2).\tau}\epsilon^{-1}_{-S^{-1}(G1+Go),-S^{-1}(G2+G0)}^*(q)
638 : !!
639 : !! TODO
640 : !! Symmetrization can be skipped if iq_bz correspond to a point in the IBZ
641 : !!
642 : !! SOURCE
643 :
644 0 : subroutine Epsm1_rotate_iqbz_inplace(epsm1, iq_bz, nomega, npwc, Gsph, Qmesh, remove_exchange)
645 :
646 : !Arguments ------------------------------------
647 : !scalars
648 : class(epsm1_t),intent(inout) :: epsm1
649 : integer,intent(in) :: iq_bz, nomega, npwc
650 : type(gsphere_t),target,intent(in) :: Gsph
651 : type(kmesh_t),intent(in) :: Qmesh
652 : logical,intent(in) :: remove_exchange
653 :
654 : !Local variables-------------------------------
655 : !scalars
656 : integer :: iw,ii,jj,iq_ibz,itim_q,isym_q,iq_loc,sg1,sg2, g0(3)
657 : !arrays
658 : real(dp) :: qbz(3)
659 : complex(gwp) :: phmsg1t,phmsg2t_star
660 0 : complex(gwp),allocatable :: work(:,:)
661 : ! *********************************************************************
662 :
663 0 : ABI_CHECK(epsm1%nomega>=nomega,'Too many frequencies required')
664 0 : ABI_CHECK(epsm1%npwe >=npwc , 'Too many G-vectors required')
665 :
666 0 : ABI_MALLOC(work, (npwc, npwc))
667 :
668 : ! Get iq_ibz, and symmetries from iq_ibz.
669 0 : call qmesh%get_BZ_item(iq_bz,qbz, iq_ibz, isym_q, itim_q, umklp=g0)
670 0 : ABI_CHECK(all(g0 == 0), "non-zero g0 is not coded")
671 :
672 : ! If out-of-memory, only epsm1%espm1(:,:,:,1) has been allocated and filled.
673 0 : iq_loc=iq_ibz; if (epsm1%mqmem==0) iq_loc=1
674 :
675 : !$OMP PARALLEL DO PRIVATE(sg2,sg1,phmsg1t,phmsg2t_star) IF (nomega > 1)
676 0 : do iw=1,nomega
677 0 : do jj=1,npwc
678 0 : sg2 = Gsph%rottb(jj,itim_q,isym_q)
679 0 : phmsg2t_star = CONJG(Gsph%phmSGt(jj,isym_q))
680 0 : do ii=1,npwc
681 0 : sg1 = Gsph%rottb(ii,itim_q,isym_q)
682 0 : phmsg1t = Gsph%phmSGt(ii,isym_q)
683 0 : work(sg1,sg2) = epsm1%epsm1(ii,jj,iw,iq_loc) * phmsg1t * phmsg2t_star
684 : end do
685 : end do
686 0 : epsm1%epsm1(:,:,iw,iq_loc) = work(:,:)
687 : end do
688 :
689 : ! Account for time-reversal
690 0 : if (itim_q==2) then
691 : !$OMP PARALLEL DO IF (nomega > 1)
692 0 : do iw=1,nomega
693 0 : call sqmat_itranspose(npwc,epsm1%epsm1(:,:,iw,iq_loc))
694 : end do
695 : end if
696 :
697 : ! Subtract the exchange contribution.
698 0 : if (remove_exchange) then
699 : !$OMP PARALLEL DO IF (nomega > 1)
700 0 : do iw=1,nomega
701 0 : do ii=1,npwc
702 0 : epsm1%epsm1(ii,ii,iw,iq_loc)=epsm1%epsm1(ii,ii,iw,iq_loc)-1.0_gwp
703 : end do
704 : end do
705 : endif
706 :
707 0 : ABI_FREE(work)
708 :
709 0 : end subroutine Epsm1_rotate_iqbz_inplace
710 : !!***
711 :
712 : !----------------------------------------------------------------------
713 :
714 : !!****f* m_screening/epsm1_from_file
715 : !! NAME
716 : !! epsm1_from_file
717 : !!
718 : !! FUNCTION
719 : !! Initialize basic dimensions and the important (small) arrays in an epsm1_t data type
720 : !! starting from a file containing either epsilon^{-1} (_SCR) or chi0 (_SUSC).
721 : !!
722 : !! INPUTS
723 : !! fname=The name of the external file used to read the matrix.
724 : !! mqmem=0 for out-of-core solution, /=0 if entire matrix has to be stored in memory.
725 : !! npwe_asked=Number of G-vector to be used in the calculation, if <=0 use Max allowed number.
726 : !! comm=MPI communicator.
727 : !!
728 : !! OUTPUT
729 : !! epsm1<epsm1_t>=The structure initialized with basic dimensions and arrays.
730 : !!
731 : !! SOURCE
732 :
733 167 : subroutine epsm1_from_file(epsm1, fname, mqmem, npwe_asked, comm)
734 :
735 : !Arguments ------------------------------------
736 : class(epsm1_t),intent(inout) :: epsm1
737 : character(len=*),intent(in) :: fname
738 : integer,intent(in) :: mqmem,npwe_asked,comm
739 :
740 : !Local variables-------------------------------
741 : !scalars
742 : integer,parameter :: master=0
743 : integer :: iw,fform,my_rank,unclassified
744 : real(dp) :: re, im, tol
745 : character(len=500) :: msg
746 : ! *********************************************************************
747 :
748 : !@epsm1_t
749 334 : my_rank = xmpi_comm_rank(comm)
750 :
751 : ! Read header from file.
752 167 : call wrtout(std_out, sjoin('epsm1_from_file- testing file: ', fname))
753 167 : call epsm1%hscr%from_file(fname, fform, comm)
754 :
755 : ! Master echoes the header.
756 306 : if (my_rank==master) call epsm1%hscr%print([std_out], 0)
757 :
758 : ! Generic Info
759 167 : epsm1%ID = 0 ! Not yet initialized as epsm1 is calculated in epsm1_mkdump.F90
760 167 : epsm1%fname = trim(fname)
761 167 : epsm1%fform = fform
762 167 : epsm1%Tordering = epsm1%Hscr%Tordering
763 :
764 : !TODO these quantitities should be checked and initialized in epsm1_mkdump
765 : !BEGIN HARDCODED
766 167 : epsm1%nI = 1
767 167 : epsm1%nJ = 1
768 167 : epsm1%ikxc = 0
769 167 : epsm1%test_type=-1
770 :
771 167 : epsm1%Hscr%headform = HSCR_LATEST_HEADFORM ! XG20090912
772 : !END HARDCODED
773 :
774 167 : epsm1%nqibz=epsm1%Hscr%nqibz
775 167 : epsm1%mqmem=mqmem ; if (mqmem/=0) epsm1%mqmem=epsm1%nqibz
776 501 : ABI_MALLOC(epsm1%qibz, (3,epsm1%nqibz))
777 4523 : epsm1%qibz(:,:) = epsm1%Hscr%qibz(:,:)
778 :
779 167 : epsm1%nqlwl = epsm1%Hscr%nqlwl
780 501 : ABI_MALLOC(epsm1%qlwl, (3, epsm1%nqlwl))
781 835 : epsm1%qlwl(:,:)=epsm1%Hscr%qlwl(:,:)
782 :
783 167 : epsm1%nomega=epsm1%Hscr%nomega
784 501 : ABI_MALLOC(epsm1%omega,(epsm1%nomega))
785 1327 : epsm1%omega(:)=epsm1%Hscr%omega(:)
786 :
787 : ! Count number of real, imaginary, and complex frequencies.
788 167 : epsm1%nomega_r = 1; epsm1%nomega_i = 0
789 167 : if (epsm1%nomega == 2) then
790 128 : epsm1%nomega_i = 1
791 : else
792 39 : unclassified = 0
793 39 : tol = tol6*Ha_eV
794 904 : do iw=2, epsm1%nomega
795 865 : re = REAL(epsm1%omega(iw))
796 865 : im = AIMAG(epsm1%omega(iw))
797 904 : if (re > tol .and. im < tol) then
798 207 : epsm1%nomega_r = iw ! Real freqs are packed in the first locations.
799 658 : else if (re < tol .and. im > tol) then
800 658 : epsm1%nomega_i = epsm1%nomega_i + 1
801 : else
802 0 : unclassified = unclassified + 1
803 : end if
804 : end do
805 39 : if (unclassified > 0) then
806 : write(msg,'(3a,i0)')&
807 0 : 'Some complex frequencies are too small to qualify as real or imaginary.',ch10,&
808 0 : 'Number of unidentified frequencies: ', unclassified
809 0 : ABI_WARNING(msg)
810 : end if
811 : end if
812 :
813 : ! Get G-vectors.
814 167 : epsm1%npwe = epsm1%Hscr%npwe
815 167 : if (npwe_asked>0) then
816 110 : if (npwe_asked > epsm1%Hscr%npwe) then
817 : write(msg,'(a,i8,2a,i8)')&
818 0 : 'Number of G-vectors saved on file is less than the value required = ',npwe_asked,ch10,&
819 0 : 'Calculation will proceed with Max available npwe = ',epsm1%Hscr%npwe
820 0 : ABI_WARNING(msg)
821 : else
822 : ! Redefine the no. of G"s for W.
823 110 : epsm1%npwe = npwe_asked
824 : end if
825 : end if
826 :
827 : ! pointer to epsm1%Hscr%gvec ?
828 501 : ABI_MALLOC(epsm1%gvec, (3,epsm1%npwe))
829 35378 : epsm1%gvec=epsm1%Hscr%gvec(:,1:epsm1%npwe)
830 :
831 167 : end subroutine epsm1_from_file
832 : !!***
833 :
834 : !----------------------------------------------------------------------
835 :
836 : !!****f* m_screening/epsm1_mkdump
837 : !! NAME
838 : !! epsm1_mkdump
839 : !!
840 : !! FUNCTION
841 : !! Dump the content of an epsm1_t data type on file.
842 : !!
843 : !! INPUTS
844 : !! id_required=Identifier of the matrix to be calculated
845 : !! Vcp<vcoul_t>=Structure gathering data on the Coulombian interaction
846 : !! ngfft(18)=Info on the FFT mesh.
847 : !! nfftot=Total number of point on the FFT mesh.
848 : !! gvec(3,npwe)=Reduced coordinates of plane waves for the response functions
849 : !! npwe=Number of plane waves.
850 : !! comm=MPI communicator.
851 : !!
852 : !! OUTPUT
853 : !!
854 : !! SOURCE
855 :
856 166 : subroutine epsm1_mkdump(epsm1, Vcp, npwe, gvec, nkxc, kxcg, id_required, approx_type, &
857 : ikxc_required, option_test, fname_dump, iomode, nfftot, ngfft, comm, &
858 0 : fxc_ADA) ! optional
859 :
860 : !Arguments ------------------------------------
861 : !scalars
862 : class(epsm1_t),intent(inout) :: epsm1
863 : integer,intent(in) :: id_required, approx_type, option_test, ikxc_required, nkxc
864 : integer,intent(in) :: iomode, nfftot, npwe, comm
865 : type(vcoul_t),intent(in) :: Vcp
866 : character(len=*),intent(in) :: fname_dump
867 : !arrays
868 : integer,intent(in) :: ngfft(18), gvec(3,npwe)
869 : complex(gwp),intent(in) :: kxcg(nfftot,nkxc)
870 : complex(gwp),intent(in), optional :: fxc_ADA(npwe*epsm1%nI,npwe*epsm1%nJ,epsm1%nqibz)
871 :
872 : !Local variables-------------------------------
873 : !scalars
874 : integer,parameter :: master=0
875 : integer :: dim_wing,iq_ibz,is_qeq0,mqmem_,npwe_asked,unt_dump,fform,rdwr,ierr,my_rank,iomode__, nprocs
876 : real(dp) :: ucvol
877 : character(len=500) :: msg
878 : character(len=fnlen) :: ofname
879 : character(len=nctk_slen) :: in_varname,out_varname
880 664 : type(hscr_t) :: Hscr_cp
881 166 : type(spectra_t) :: spectra
882 : integer :: comm__
883 : type(c_ptr) :: void_ptr
884 : integer(kind=XMPI_ADDRESS_KIND) :: count
885 : type(xcomm_t) :: xcomm
886 : !arrays
887 : integer :: units(2)
888 : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
889 166 : complex(gwp),allocatable :: tmp_epsm1(:,:,:)
890 166 : complex(dp),allocatable :: dummy_lwing(:,:,:),dummy_uwing(:,:,:),dummy_head(:,:,:)
891 : ! *********************************************************************
892 :
893 0 : ABI_CHECK_IEQ(id_required, 4, 'Value of id_required not coded')
894 166 : ABI_CHECK_IEQ(npwe, epsm1%npwe, "mismatch in npwe")
895 :
896 498 : units = [std_out, ab_out]
897 :
898 166 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
899 166 : call metric(gmet,gprimd,-1,rmet,Vcp%rprimd,ucvol)
900 :
901 : ! if (epsm1%ID/=0) call reset_Epsilonm1(epsm1)
902 166 : epsm1%ID=id_required
903 :
904 166 : ofname = fname_dump
905 166 : in_varname = ncname_from_id(epsm1%hscr%id)
906 166 : out_varname = ncname_from_id(id_required)
907 :
908 : !write(std_out,*)'epsm1%ID: ',epsm1%ID,', epsm1%Hscr%ID: ',epsm1%Hscr%ID
909 :
910 166 : if (epsm1%ID == epsm1%Hscr%ID) then
911 : ! The two-point function we are asking for is already stored on file.
912 : ! According to mqmem either read and store the entire matrix in memory or do nothing.
913 :
914 161 : if (epsm1%mqmem > 0) then
915 : ! In-core solution.
916 : epsm1%use_mpi_shared_win = .False.
917 : #ifdef HAVE_MPI_ALLOCATE_SHARED_CPTR
918 145 : epsm1%use_mpi_shared_win = nprocs > 1
919 : !epsm1%use_mpi_shared_win = .True. ! This to test shared memory with one proc.
920 : #endif
921 : !epsm1%use_mpi_shared_win = .False. ! This to go back to the old non-scalable version.
922 :
923 145 : iomode__ = iomode
924 145 : if (iomode__ == IO_MODE_MPI) then
925 0 : ABI_WARNING("SUSC files with IO_MODE_MPI is buggy. Using Fortran IO")
926 0 : iomode__ = IO_MODE_FORTRAN
927 : end if
928 :
929 145 : write(msg,'(a,f12.1,a)')' Memory for epsm1%epsm1: ',two*gwp*npwe**2*epsm1%nomega*epsm1%nqibz*b2Mb,' [Mb] <<< MEM'
930 145 : call wrtout(std_out, msg)
931 :
932 145 : if (.not. epsm1%use_mpi_shared_win) then
933 115 : call epsm1%shared_comm%set_to_self()
934 :
935 115 : if (nprocs > 1) then
936 : msg = strcat("- WARNING: Cannot use MPI shared memory as MPI library does not support MPI_WIN_ALLOCATE_SHARED_CPTR", ch10, &
937 0 : "- Memory for epsm1 will increase with nprocs per node!")
938 0 : ABI_WARNING(msg)
939 0 : call wrtout(ab_out, msg)
940 : end if
941 :
942 690 : ABI_MALLOC_OR_DIE(epsm1%epsm1, (npwe, npwe, epsm1%nomega, epsm1%nqibz), ierr)
943 115 : call read_screening(in_varname, epsm1%fname, epsm1%npwe, epsm1%nqibz, epsm1%nomega, epsm1%epsm1, iomode__, comm)
944 :
945 : else
946 30 : call wrtout(std_out, "- HAPPY: Using MPI shared memory, memory for epsm1 won't increase with nprocs per node!")
947 30 : comm__ = comm
948 30 : xcomm = xcomm_from_mpi_int(comm__)
949 30 : epsm1%shared_comm = xcomm%split_type()
950 :
951 30 : count = _MOK(2 * npwe) * _MOK(npwe) * _MOK(epsm1%nomega * epsm1%nqibz)
952 30 : call epsm1%shared_comm%allocate_shared_master(count, gwp, xmpi_info_null, void_ptr, epsm1%epsm1_win)
953 150 : call c_f_pointer(void_ptr, epsm1%epsm1, shape=[npwe, npwe, epsm1%nomega, epsm1%nqibz])
954 :
955 : ! Only one proc in shared_comm reads from file.
956 30 : call xmpi_win_fence(XMPI_MODE_NOPRECEDE, epsm1%epsm1_win, ierr) ! Start the RMA epoch.
957 30 : ABI_CHECK_MPI(ierr, "")
958 30 : if (epsm1%shared_comm%me == 0) then
959 10 : call read_screening(in_varname, epsm1%fname, epsm1%npwe, epsm1%nqibz, epsm1%nomega, epsm1%epsm1, iomode__, xmpi_comm_self)
960 : end if
961 30 : call xmpi_win_fence(XMPI_MODE_NOSUCCEED, epsm1%epsm1_win, ierr) ! Close the RMA epoch.
962 30 : ABI_CHECK_MPI(ierr, "")
963 30 : call xcomm%free()
964 : end if
965 :
966 : else
967 : ! Out-of-core solution
968 16 : ABI_COMMENT("mqmem==0 => allocating a single q-slice of (W|chi0) (slower but less memory).")
969 : continue
970 : end if
971 :
972 161 : return
973 :
974 : else
975 : ! === The matrix stored on file do not correspond to the quantity required ===
976 : ! * Presently only the transformation chi0 => e^-1 is coded
977 : ! * According to epsm1%mqmem either calculate e^-1 dumping the result to a file
978 : ! for a subsequent use or calculate e^-1 keeping everything in memory.
979 :
980 5 : if (epsm1%mqmem == 0) then
981 : ! Open file and write the header for the SCR file. For the moment only master works.
982 :
983 1 : if (my_rank==master) then
984 1 : if (iomode == IO_MODE_ETSF) then
985 1 : ofname = nctk_ncify(ofname)
986 1 : NCF_CHECK(nctk_open_create(unt_dump, ofname, xmpi_comm_self))
987 : else
988 0 : if (open_file(ofname,msg,newunit=unt_dump,form="unformatted",status="unknown",action="write") /= 0) then
989 0 : ABI_ERROR(msg)
990 : end if
991 : end if
992 1 : call wrtout(std_out,sjoin('epsm1_mkdump: calculating and writing epsilon^-1 matrix on file: ',ofname))
993 :
994 : ! Update the entries in the header that have been modified.
995 : ! TODO, write function to return title, just for info
996 1 : call epsm1%Hscr%copy(Hscr_cp)
997 1 : Hscr_cp%ID = id_required
998 1 : Hscr_cp%ikxc = ikxc_required
999 1 : Hscr_cp%test_type = option_test
1000 1 : Hscr_cp%titles(1) = 'SCR file: epsilon^-1'
1001 1 : Hscr_cp%titles(2) = 'TESTPARTICLE'
1002 : ! Treat the case in which a smaller matrix is used.
1003 1 : Hscr_cp%npwe = npwe
1004 :
1005 1 : rdwr=2; fform=Hscr_cp%fform
1006 1 : call hscr_cp%io(fform,rdwr,unt_dump,xmpi_comm_self,master,iomode)
1007 1 : call Hscr_cp%free()
1008 :
1009 5 : ABI_MALLOC_OR_DIE(tmp_epsm1, (npwe, npwe, epsm1%nomega), ierr)
1010 :
1011 7 : do iq_ibz=1,epsm1%nqibz
1012 6 : is_qeq0=0
1013 6 : if (normv(epsm1%qibz(:,iq_ibz),gmet,'G')<GW_TOLQ0) is_qeq0=1
1014 :
1015 6 : call read_screening(in_varname,epsm1%fname,npwe,1,epsm1%nomega,tmp_epsm1,iomode,xmpi_comm_self,iqiA=iq_ibz)
1016 :
1017 6 : dim_wing=0; if (is_qeq0==1) dim_wing=3
1018 30 : ABI_MALLOC(dummy_lwing,(npwe*epsm1%nI,epsm1%nomega,dim_wing))
1019 30 : ABI_MALLOC(dummy_uwing,(npwe*epsm1%nJ,epsm1%nomega,dim_wing))
1020 30 : ABI_MALLOC(dummy_head,(dim_wing,dim_wing,epsm1%nomega))
1021 :
1022 6 : if (approx_type<2 .or. approx_type>3) then
1023 6 : ABI_WARNING('Entering out-of core RPA or Kxc branch')
1024 : call make_epsm1_driver(iq_ibz,dim_wing,npwe,epsm1%nI,epsm1%nJ,epsm1%nomega,epsm1%omega,&
1025 : approx_type,option_test,Vcp,nfftot,ngfft,nkxc,kxcg,gvec,dummy_head,&
1026 6 : dummy_lwing,dummy_uwing,tmp_epsm1,spectra,xmpi_comm_self,xmpi_undefined)
1027 : else
1028 0 : ABI_WARNING('Entering out-of core fxc_ADA branch')
1029 : call make_epsm1_driver(iq_ibz,dim_wing,npwe,epsm1%nI,epsm1%nJ,epsm1%nomega,epsm1%omega,&
1030 : approx_type,option_test,Vcp,nfftot,ngfft,nkxc,kxcg,gvec,dummy_head,&
1031 : dummy_lwing,dummy_uwing,tmp_epsm1,spectra,xmpi_comm_self, xmpi_undefined, &
1032 0 : fxc_ADA=fxc_ADA(:,:,iq_ibz))
1033 : end if
1034 :
1035 6 : ABI_FREE(dummy_head)
1036 6 : ABI_FREE(dummy_uwing)
1037 6 : ABI_FREE(dummy_lwing)
1038 :
1039 6 : if (is_qeq0==1) then
1040 1 : call spectra%repr(msg)
1041 1 : call wrtout(units, msg)
1042 : end if
1043 6 : call spectra%free()
1044 :
1045 7 : call write_screening(out_varname,unt_dump,iomode,npwe,epsm1%nomega,iq_ibz,tmp_epsm1)
1046 : end do
1047 :
1048 1 : if (iomode == IO_MODE_ETSF) then
1049 1 : NCF_CHECK(nf90_close(unt_dump))
1050 : else
1051 0 : close(unt_dump)
1052 : endif
1053 :
1054 1 : ABI_FREE(tmp_epsm1)
1055 : end if !master
1056 :
1057 : ! Master broadcasts ofname.
1058 : ! NOTE: A synchronization is required here, else the other procs start to read the
1059 : ! SCR file before it is written by the master. xmpi_bcast will synch the procs.
1060 1 : call xmpi_bcast(ofname, master, comm, ierr)
1061 :
1062 : ! Now epsm1% "belongs" to the file "ofname", thus
1063 : ! each proc has to destroy and re-initialize the object.
1064 1 : call epsm1%free()
1065 :
1066 1 : mqmem_=epsm1%mqmem; npwe_asked=npwe
1067 1 : call epsm1%from_file(ofname, mqmem_, npwe_asked, comm)
1068 :
1069 : ! Now epsm1% has been reinitialized and ready-to-use.
1070 1 : epsm1%id = id_required
1071 2 : call epsm1%print([std_out], 0)
1072 : else
1073 : ! ========================
1074 : ! === In-core solution ===
1075 : ! ========================
1076 : ! In-core solution.
1077 : epsm1%use_mpi_shared_win = .False.
1078 : #ifdef HAVE_MPI_ALLOCATE_SHARED_CPTR
1079 4 : epsm1%use_mpi_shared_win = nprocs > 1
1080 : !epsm1%use_mpi_shared_win = .True. ! This to test shared memory with one proc.
1081 : #endif
1082 : !epsm1%use_mpi_shared_win = .False. ! This to go back to the old non-scalable version.
1083 :
1084 4 : if (.not. epsm1%use_mpi_shared_win) then
1085 4 : call epsm1%shared_comm%set_to_self()
1086 :
1087 : if (nprocs > 1) then
1088 : msg = strcat("- WARNING: Cannot use MPI shared memory as MPI library does not support MPI_WIN_ALLOCATE_SHARED_CPTR", ch10, &
1089 : "- Memory for epsm1 will increase with nprocs per node!")
1090 : ABI_WARNING(msg)
1091 : call wrtout(ab_out, msg)
1092 : end if
1093 :
1094 24 : ABI_MALLOC_OR_DIE(epsm1%epsm1, (npwe,npwe,epsm1%nomega,epsm1%nqibz), ierr)
1095 4 : call read_screening(in_varname, epsm1%fname, npwe, epsm1%nqibz, epsm1%nomega, epsm1%epsm1, iomode, comm)
1096 :
1097 : else
1098 0 : call wrtout(std_out, "- HAPPY: Using MPI shared memory, memory for epsm1 won't increase with nprocs per node!")
1099 0 : comm__ = comm
1100 0 : xcomm = xcomm_from_mpi_int(comm__)
1101 0 : epsm1%shared_comm = xcomm%split_type()
1102 :
1103 0 : count = _MOK(2 * npwe) * _MOK(npwe) * _MOK(epsm1%nomega * epsm1%nqibz)
1104 0 : call epsm1%shared_comm%allocate_shared_master(count, gwp, xmpi_info_null, void_ptr, epsm1%epsm1_win)
1105 0 : call c_f_pointer(void_ptr, epsm1%epsm1, shape=[npwe, npwe, epsm1%nomega, epsm1%nqibz])
1106 :
1107 : ! Only one proc in shared_comm reads from file.
1108 0 : call xmpi_win_fence(XMPI_MODE_NOPRECEDE, epsm1%epsm1_win, ierr) ! Start the RMA epoch.
1109 0 : ABI_CHECK_MPI(ierr, "")
1110 0 : if (epsm1%shared_comm%me == 0) then
1111 0 : call read_screening(in_varname, epsm1%fname, epsm1%npwe, epsm1%nqibz, epsm1%nomega, epsm1%epsm1, iomode__, xmpi_comm_self)
1112 : end if
1113 0 : call xmpi_win_fence(XMPI_MODE_NOSUCCEED, epsm1%epsm1_win, ierr) ! Close the RMA epoch.
1114 0 : ABI_CHECK_MPI(ierr, "")
1115 0 : call xcomm%free()
1116 : end if
1117 :
1118 26 : do iq_ibz=1,epsm1%nqibz
1119 22 : is_qeq0=0; if (normv(epsm1%qibz(:,iq_ibz),gmet,'G')<GW_TOLQ0) is_qeq0=1
1120 :
1121 22 : dim_wing=0; if (is_qeq0==1) dim_wing=3 ! FIXME
1122 110 : ABI_MALLOC(dummy_lwing,(npwe*epsm1%nI,epsm1%nomega,dim_wing))
1123 110 : ABI_MALLOC(dummy_uwing,(npwe*epsm1%nJ,epsm1%nomega,dim_wing))
1124 110 : ABI_MALLOC(dummy_head,(dim_wing,dim_wing,epsm1%nomega))
1125 :
1126 22 : comm__ = comm
1127 22 : if (epsm1%use_mpi_shared_win) comm__ = epsm1%shared_comm%value
1128 :
1129 22 : if (approx_type<2 .or. approx_type>3) then
1130 22 : ABI_COMMENT('Entering in-core RPA and Kxc branch')
1131 : call make_epsm1_driver(iq_ibz,dim_wing,npwe,epsm1%nI,epsm1%nJ,epsm1%nomega,epsm1%omega,&
1132 : approx_type,option_test,Vcp,nfftot,ngfft,nkxc,kxcg,gvec,dummy_head,&
1133 22 : dummy_lwing,dummy_uwing,epsm1%epsm1(:,:,:,iq_ibz),spectra,comm__, epsm1%epsm1_win)
1134 : else
1135 0 : ABI_COMMENT('Entering in-core fxc_ADA branch')
1136 : call make_epsm1_driver(iq_ibz,dim_wing,npwe,epsm1%nI,epsm1%nJ,epsm1%nomega,epsm1%omega,&
1137 : approx_type,option_test,Vcp,nfftot,ngfft,nkxc,kxcg,gvec,dummy_head,&
1138 : dummy_lwing,dummy_uwing,epsm1%epsm1(:,:,:,iq_ibz),spectra,comm__, epsm1%epsm1_win, &
1139 0 : fxc_ADA=fxc_ADA(:,:,iq_ibz))
1140 : end if
1141 :
1142 22 : ABI_FREE(dummy_lwing)
1143 22 : ABI_FREE(dummy_uwing)
1144 22 : ABI_FREE(dummy_head)
1145 :
1146 22 : if (is_qeq0==1) then
1147 4 : call spectra%repr(msg)
1148 4 : call wrtout(units, msg)
1149 : end if
1150 :
1151 26 : call spectra%free()
1152 : end do
1153 :
1154 4 : epsm1%id = id_required
1155 8 : call epsm1%print([std_out], 0)
1156 : end if
1157 : end if
1158 :
1159 498 : end subroutine epsm1_mkdump
1160 : !!***
1161 :
1162 132 : subroutine epsm1_malloc_epsm1_qbz(epsm1, npwc, nomega)
1163 :
1164 : !Arguments ------------------------------------
1165 : class(epsm1_t),intent(inout) :: epsm1
1166 : integer,intent(in) :: npwc, nomega
1167 :
1168 : !Local variables-------------------------------
1169 : integer :: ierr
1170 : integer(kind=XMPI_ADDRESS_KIND) :: count
1171 : type(c_ptr) :: void_ptr
1172 : ! *********************************************************************
1173 :
1174 132 : ABI_CHECK(epsm1%npwe >= npwc, 'Too many G-vectors required')
1175 :
1176 132 : if (.not. epsm1%use_mpi_shared_win) then
1177 420 : ABI_MALLOC_OR_DIE(epsm1%epsm1_qbz, (npwc, npwc, nomega), ierr)
1178 : else
1179 48 : count = _MOK(2 * npwc) * _MOK(npwc) * _MOK(nomega)
1180 48 : call epsm1%shared_comm%allocate_shared_master(count, gwp, xmpi_info_null, void_ptr, epsm1%epsm1_qbz_win)
1181 192 : call c_f_pointer(void_ptr, epsm1%epsm1_qbz, shape=[npwc, npwc, nomega])
1182 : end if
1183 :
1184 132 : end subroutine epsm1_malloc_epsm1_qbz
1185 : !!***
1186 :
1187 : !----------------------------------------------------------------------
1188 :
1189 : !!****f* m_screening/epsm1_free_epsm1_qbz
1190 : !! NAME
1191 : !! epsm1_free_epsm1_qbz
1192 : !!
1193 : !! FUNCTION
1194 : !! Free the internal buffer %epsm1_qbz
1195 :
1196 460 : subroutine epsm1_free_epsm1_qbz(epsm1)
1197 :
1198 : !Arguments ------------------------------------
1199 : class(epsm1_t),intent(inout) :: epsm1
1200 :
1201 : !Local variables-------------------------------
1202 : integer :: ierr
1203 : ! *********************************************************************
1204 :
1205 460 : if (.not. epsm1%use_mpi_shared_win) then
1206 304 : ABI_SFREE_PTR(epsm1%epsm1_qbz)
1207 : else
1208 156 : if (epsm1%epsm1_qbz_win /= xmpi_undefined) then
1209 48 : call xmpi_win_free(epsm1%epsm1_qbz_win, ierr)
1210 48 : nullify(epsm1%epsm1_qbz)
1211 : end if
1212 : end if
1213 :
1214 460 : end subroutine epsm1_free_epsm1_qbz
1215 : !!***
1216 :
1217 : !----------------------------------------------------------------------
1218 :
1219 : !!****f* m_screening/get_epsm1
1220 : !! NAME
1221 : !! get_epsm1
1222 : !!
1223 : !! FUNCTION
1224 : !! Work in progress but the main is idea is as follows:
1225 : !!
1226 : !! Return the symmetrized inverse dielectric matrix.
1227 : !! This method implements both in-core and the out-of-core solution
1228 : !! In the later, epsilon^-1 or chi0 are read from file.
1229 : !! It is possible to specify options to retrieve (RPA |TDDDT, [TESTCHARGE|TESTPARTICLE]).
1230 : !! All dimensions are already initialized in the epsm1% object, this method
1231 : !! should act as a wrapper around rdscr and make_epsm1_driver. A better
1232 : !! implementation will be done in the following once the coding of file handlers is completed.
1233 : !!
1234 : !! INPUTS
1235 : !! Vcp<vcoul_t>=Structure gathering data on the Coulombian interaction
1236 : !! iqibzA[optional]=Index of the q-point to be read from file (only for out-of-memory solutions)
1237 : !! iomode=option definig the file format.
1238 : !! option_test
1239 : !! comm=MPI communicator.
1240 : !!
1241 : !! OUTPUT
1242 : !! epsm1%epsm1
1243 : !!
1244 : !! TODO
1245 : !! Remove this routine. Now everything should be done with epsm1_mkdump
1246 : !!
1247 : !! SOURCE
1248 :
1249 290 : subroutine get_epsm1(epsm1,Vcp,approx_type,option_test,iomode,comm,iqibzA)
1250 :
1251 : !Arguments ------------------------------------
1252 : !scalars
1253 : class(epsm1_t),intent(inout) :: epsm1
1254 : integer,intent(in) :: iomode,option_test,approx_type,comm
1255 : integer,optional,intent(in) :: iqibzA
1256 : type(vcoul_t),intent(in) :: Vcp
1257 :
1258 : !Local variables-------------------------------
1259 : !scalars
1260 : integer :: my_approx_type,my_option_test,ng,ierr
1261 : ! *********************************************************************
1262 :
1263 290 : my_approx_type = approx_type; my_option_test = option_test
1264 :
1265 : ! Vcp not yet used.
1266 290 : ng = Vcp%ng
1267 :
1268 580 : select case (epsm1%mqmem)
1269 : case (0)
1270 : ! Out-of-core solution
1271 290 : ABI_SFREE_PTR(epsm1%epsm1)
1272 1450 : ABI_MALLOC_OR_DIE(epsm1%epsm1,(epsm1%npwe,epsm1%npwe,epsm1%nomega,1), ierr)
1273 :
1274 290 : call read_screening(em1_ncname,epsm1%fname,epsm1%npwe,epsm1%nqibz,epsm1%nomega,epsm1%epsm1,iomode,comm,iqiA=iqibzA)
1275 :
1276 290 : if (epsm1%id == 4) then
1277 : ! If q-slice of epsilon^-1 has been read then return
1278 : !call epsm1%print()
1279 : return
1280 : else
1281 0 : ABI_ERROR(sjoin('Wrong epsm1%ID', itoa(epsm1%id)))
1282 : end if
1283 :
1284 : case default
1285 : ! In-core solution.
1286 290 : ABI_ERROR("you should not be here")
1287 : end select
1288 :
1289 : end subroutine get_epsm1
1290 : !!***
1291 :
1292 : !----------------------------------------------------------------------
1293 :
1294 : !!****f* m_screening/decompose_epsm1
1295 : !! NAME
1296 : !! decompose_epsm1
1297 : !!
1298 : !! FUNCTION
1299 : !! Decompose the complex symmetrized dielectric
1300 : !!
1301 : !! INPUTS
1302 : !!
1303 : !! OUTPUT
1304 : !!
1305 : !! SOURCE
1306 :
1307 0 : subroutine decompose_epsm1(epsm1, iq_ibz, eigs)
1308 :
1309 : !Arguments ------------------------------------
1310 : !scalars
1311 : class(epsm1_t),intent(in) :: epsm1
1312 : integer,intent(in) :: iq_ibz
1313 : !arrays
1314 : complex(dp),intent(out) :: eigs(epsm1%npwe,epsm1%nomega)
1315 :
1316 : !Local variables-------------------------------
1317 : !scalars
1318 : integer :: info,lwork,iw,negw,ig1,ig2,idx,sdim,npwe,ierr
1319 : character(len=500) :: msg
1320 : !arrays
1321 0 : real(dp),allocatable :: ww(:),rwork(:)
1322 0 : complex(dp),allocatable :: work(:),Adpp(:),eigvec(:,:),Afull(:,:),vs(:,:),wwc(:)
1323 0 : logical,allocatable :: bwork(:)
1324 : logical :: sortcplx !BUG in abilint
1325 : ! *********************************************************************
1326 :
1327 0 : ABI_CHECK(epsm1%mqmem /= 0, 'mqmem==0 not implemented')
1328 :
1329 0 : npwe = epsm1%npwe
1330 :
1331 0 : do iw=1,epsm1%nomega
1332 :
1333 0 : if (ABS(REAL(epsm1%omega(iw)))>0.00001) then
1334 : ! Eigenvalues for a generic complex matrix
1335 0 : lwork=4*2*npwe
1336 0 : ABI_MALLOC(wwc,(npwe))
1337 0 : ABI_MALLOC(work,(lwork))
1338 0 : ABI_MALLOC(rwork,(npwe))
1339 0 : ABI_MALLOC(bwork,(npwe))
1340 0 : ABI_MALLOC(vs,(npwe,npwe))
1341 0 : ABI_MALLOC(Afull,(npwe,npwe))
1342 :
1343 0 : Afull=epsm1%epsm1(:,:,iw,iq_ibz)
1344 :
1345 : !for the moment no sort, maybe here I should sort using the real part?
1346 0 : call ZGEES('V','N',sortcplx,npwe,Afull,npwe,sdim,wwc,vs,npwe,work,lwork,rwork,bwork,info)
1347 0 : if (info/=0) then
1348 0 : ABI_ERROR(sjoin("ZGEES returned info:",itoa(info)))
1349 : end if
1350 :
1351 0 : eigs(:,iw)=wwc(:)
1352 :
1353 0 : ABI_FREE(wwc)
1354 0 : ABI_FREE(work)
1355 0 : ABI_FREE(rwork)
1356 0 : ABI_FREE(bwork)
1357 0 : ABI_FREE(vs)
1358 0 : ABI_FREE(Afull)
1359 :
1360 : else
1361 : ! Hermitian version.
1362 0 : lwork=2*npwe-1
1363 0 : ABI_MALLOC(ww,(npwe))
1364 0 : ABI_MALLOC(work,(lwork))
1365 0 : ABI_MALLOC(rwork,(3*npwe-2))
1366 0 : ABI_MALLOC(eigvec,(npwe,npwe))
1367 0 : ABI_MALLOC_OR_DIE(Adpp,(npwe*(npwe+1)/2), ierr)
1368 :
1369 : idx=0 ! Pack the matrix
1370 0 : do ig2=1,npwe
1371 0 : do ig1=1,ig2
1372 0 : idx=idx+1
1373 0 : Adpp(idx)=epsm1%epsm1(ig1,ig2,iw,iq_ibz)
1374 : end do
1375 : end do
1376 :
1377 : ! For the moment we require also the eigenvectors.
1378 0 : call ZHPEV('V','U',npwe,Adpp,ww,eigvec,npwe,work,rwork,info)
1379 0 : if (info/=0) then
1380 0 : ABI_ERROR(sjoin('ZHPEV returned info=', itoa(info)))
1381 : end if
1382 :
1383 0 : negw=(COUNT((REAL(ww)<tol6)))
1384 0 : if (negw/=0) then
1385 : write(msg,'(a,i5,a,i3,a,f8.4)')&
1386 0 : 'Found negative eigenvalues. No. ',negw,' at iq_ibz= ',iq_ibz,' minval= ',MINVAL(REAL(ww))
1387 0 : ABI_WARNING(msg)
1388 : end if
1389 :
1390 0 : eigs(:,iw)=ww(:)
1391 :
1392 0 : ABI_FREE(ww)
1393 0 : ABI_FREE(work)
1394 0 : ABI_FREE(rwork)
1395 0 : ABI_FREE(eigvec)
1396 0 : ABI_FREE(Adpp)
1397 : end if
1398 : end do !iw
1399 :
1400 : ! contains
1401 : ! function sortcplx(carg) result(res)
1402 : ! implicit none
1403 : ! complex(dp),intent(in) :: carg
1404 : ! logical :: res
1405 : ! res=.TRUE.
1406 : ! end function sortcplx
1407 :
1408 0 : end subroutine decompose_epsm1
1409 : !!***
1410 :
1411 : !----------------------------------------------------------------------
1412 :
1413 : !!****f* m_screening/make_epsm1_driver
1414 : !! NAME
1415 : !! make_epsm1_driver
1416 : !!
1417 : !! FUNCTION
1418 : !! Driver routine to calculate the inverse symmetrical dielectric matrix starting
1419 : !! from the irreducible polarizability. The routine considers a single q-point, and
1420 : !! performs the following tasks:
1421 : !!
1422 : !! 1) Calculate $\tilde\epsilon^{-1}$ using different approximations:
1423 : !! * RPA
1424 : !! * ALDA within TDDFT
1425 : !!
1426 : !! 2) Use a special treatment of non-Analytic behavior of heads and wings in reciprocal space
1427 : !! calculating these quantities for different small q-directions specified by the user
1428 : !! (Not yet operative)
1429 : !!
1430 : !! 3) Output the electron energy loss function and the macroscopic dielectric function with and
1431 : !! without local field effects (only if non-zero real frequencies are available)
1432 : !!
1433 : !! INPUTS
1434 : !! iq_ibz=index of the q-point in the array Vcp%qibz where epsilon^-1 has to be calculated
1435 : !! dim_wing=Dimension of the wings (0 or 3 if q-->0)
1436 : !! npwe=Number of G-vectors in chi0.
1437 : !! nI,nJ=Number of rows/columns in chi0_ij (1,1 if collinear case)
1438 : !! nomega=Number of frequencies.
1439 : !! omega(nomega)=Frequencines in Hartree
1440 : !! approx_type=Integer flag defining the type of approximation
1441 : !! == 0 for RPA ==
1442 : !! == 1 for TDDFT ==
1443 : !! option_test=Only for TDDFT:
1444 : !! == 0 for TESTPARTICLE ==
1445 : !! == 1 for TESTELECTRON ==
1446 : !! Vcp<vcoul_t>=Structure gathering data on the Coulombian interaction
1447 : !! nfftot=Total number of points in the FFT mesh.
1448 : !! ngfft(18)=Info on the FFT mesh.
1449 : !! nkxc=Dimension of the kernel in reciprocal space. 0 if kernel is not needed
1450 : !! kxcg(nfftot,nkxc)=TDDFT kernel in reciprocal space on the FFT mesh. Used only if approx_type==1
1451 : !! gvec(3,npwe)=G-vectors
1452 : !! comm=MPI communicator.
1453 : !! chi0_lwing(npwe*nI,nomega,dim_wing)=Lower wings of chi0 (only for q-->0)
1454 : !! chi0_uwing(npwe*nJ,nomega,dim_wing)=Upper wings of chi0 (only for q-->0)
1455 : !! chi0_head(dim_wing,dim_wing,nomega)=Head of of chi0 (only for q-->0)
1456 : !!
1457 : !! OUTPUT
1458 : !! spectra<spectra_t>Object containing e_macro(w) and EELS(w)
1459 : !!
1460 : !! SIDE EFFECTS
1461 : !! chi0(npwe*nI,npwe*nJ,nomega): in input the irreducible polarizability, in output
1462 : !! the symmetrized inverse dielectric matrix.
1463 : !!
1464 : !! SOURCE
1465 :
1466 950 : subroutine make_epsm1_driver(iq_ibz, dim_wing, npwe, nI, nJ, nomega, omega,&
1467 950 : approx_type, option_test, Vcp, nfftot, ngfft, nkxc, kxcg, gvec, &
1468 950 : chi0_head, chi0_lwing, chi0_uwing, chi0, spectra, comm, epsm1_win, &
1469 950 : fxc_ADA, rhor) ! optional argument
1470 :
1471 : !Arguments ------------------------------------
1472 : !scalars
1473 : integer,intent(in) :: iq_ibz,nI,nJ,npwe,nomega,dim_wing,approx_type,option_test,nkxc,nfftot,comm, epsm1_win
1474 : real(dp),intent(in),optional :: rhor
1475 : type(vcoul_t),target,intent(in) :: Vcp
1476 : type(spectra_t),intent(out) :: Spectra
1477 : !arrays
1478 : integer,intent(in) :: ngfft(18),gvec(3,npwe)
1479 : complex(gwp),intent(in) :: kxcg(nfftot,nkxc)
1480 : complex(dp),intent(in) :: omega(nomega)
1481 : complex(dp),intent(inout) :: chi0_lwing(:,:,:) !(npwe*nI,nomega,dim_wing)
1482 : complex(dp),intent(inout) :: chi0_uwing(:,:,:) !(npwe*nJ,nomega,dim_wing)
1483 : complex(dp),intent(inout) :: chi0_head(:,:,:) !(dim_wing,dim_wing,nomega)
1484 : complex(gwp),intent(inout) :: chi0(npwe*nI,npwe*nJ,nomega)
1485 : complex(gwp),intent(in),optional :: fxc_ADA(npwe*nI,npwe*nJ)
1486 :
1487 : !Local variables-------------------------------
1488 : !scalars
1489 : integer,parameter :: master=0
1490 : integer :: i1,i2,ig1,ig2,io,ierr,irank,my_nqlwl !iqlwl
1491 : integer :: nor,my_rank,nprocs,g1mg2_idx
1492 : real(dp) :: ucvol
1493 : logical :: is_qeq0, use_MPI, use_mpi_shared_win
1494 : character(len=500) :: msg
1495 : !arrays
1496 1900 : integer :: omega_distrb(nomega)
1497 950 : integer,allocatable :: istart(:),istop(:)
1498 : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
1499 950 : real(dp),allocatable :: eelf(:,:),tmp_eelf(:)
1500 950 : complex(dp),allocatable :: epsm_lf(:,:),epsm_nlf(:,:),tmp_lf(:),tmp_nlf(:)
1501 950 : complex(dp),allocatable :: buffer_lwing(:,:),buffer_uwing(:,:)
1502 950 : complex(gwp),allocatable :: kxcg_mat(:,:)
1503 :
1504 : !bootstrap and LR
1505 : integer :: istep,nstep
1506 : real(dp) :: conv_err, alpha, Zr, qpg2(3), qpg2_nrm, cpu, wall, gflops
1507 : real(gwp) :: chi00_head, fxc_head
1508 950 : complex(gwp),allocatable :: vfxc_boot(:,:), vfxc_boot0(:,:), vfxc_lr(:,:), vfxc_tmp(:,:), chi0_tmp(:,:), chi0_save(:,:,:)
1509 950 : complex(gwp), contiguous, pointer :: vc_sqrt(:)
1510 : ! *************************************************************************
1511 :
1512 950 : if (nI/=1.or.nJ/=1) then
1513 0 : ABI_ERROR("nI or nJ=/1 not yet implemented")
1514 : end if
1515 :
1516 950 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
1517 950 : call cwtime(cpu, wall, gflops, "start")
1518 :
1519 950 : call metric(gmet,gprimd,-1,rmet,Vcp%rprimd,ucvol)
1520 :
1521 950 : is_qeq0 = normv(Vcp%qibz(:,iq_ibz),gmet,'G') < GW_TOLQ0
1522 :
1523 6713 : omega_distrb = my_rank
1524 950 : use_mpi_shared_win = (epsm1_win /= xmpi_undefined)
1525 950 : use_MPI = .FALSE.
1526 950 : use_MPI = nprocs >= nomega ! Parallelism is not used
1527 :
1528 : ! FIXME: MPI mode is temporarily disabled here because we need to know if
1529 : ! screening is allocated in shared memory or not.
1530 : ! Perhaps now it makes mores sense to use Scalapack/ELPA instead of parallelizing the loop over frequencies
1531 : !use_MPI = .FALSE.
1532 :
1533 950 : if (use_MPI) then
1534 : ! Initialize distribution table for frequencies.
1535 891 : ABI_MALLOC(istart,(nprocs))
1536 594 : ABI_MALLOC(istop,(nprocs))
1537 297 : call xmpi_split_work2_i4b(nomega,nprocs,istart,istop)
1538 852 : omega_distrb(:)=xmpi_undefined_rank
1539 1196 : do irank=0,nprocs-1
1540 899 : i1 = istart(irank+1)
1541 899 : i2 = istop (irank+1)
1542 1751 : if (i1<=i2) omega_distrb(i1:i2) = irank
1543 : end do
1544 297 : ABI_FREE(istart)
1545 297 : ABI_FREE(istop)
1546 : end if
1547 :
1548 : ! Initialize container for spectral results
1549 3677 : do nor=1,nomega
1550 3677 : if (ABS(AIMAG(omega(nor)))>1.e-3) EXIT
1551 : end do
1552 950 : nor=nor-1; if (nor==0) nor = 1 ! only imag !?
1553 :
1554 950 : if (dim_wing==3) then
1555 151 : call wrtout(std_out,' Analyzing long wavelength limit for several q')
1556 564 : call spectra_init(Spectra,nor,REAL(omega(1:nor)),Vcp%nqlwl,Vcp%qlwl)
1557 151 : my_nqlwl = 1
1558 : !my_nqlwl = dim_wing ! TODO
1559 : !ABI_CHECK(dim_wing==SIZE(Vcp%vcqlwl_sqrt,DIM=2),"WRONG DIMS")
1560 : else
1561 3119 : call spectra_init(Spectra,nor,REAL(omega(1:nor)),1,Vcp%qibz(:,iq_ibz))
1562 799 : my_nqlwl = 1
1563 : end if
1564 : !
1565 : ! NOTE: all processors have to perform this operation in order to have the
1566 : ! epsm1 matrix when performing a sigma calculation starting with the file _SUS
1567 : !
1568 : ! Temporary arrays to store spectra.
1569 3800 : ABI_MALLOC(epsm_lf,(nomega,my_nqlwl))
1570 2850 : ABI_MALLOC(epsm_nlf,(nomega,my_nqlwl))
1571 3800 : ABI_MALLOC(eelf,(nomega,my_nqlwl))
1572 21089 : epsm_lf=czero; epsm_nlf=czero; eelf=zero
1573 :
1574 : ! Temporary arrays used to store output results.
1575 2850 : ABI_MALLOC(tmp_lf, (my_nqlwl))
1576 1900 : ABI_MALLOC(tmp_nlf, (my_nqlwl))
1577 2850 : ABI_MALLOC(tmp_eelf, (my_nqlwl))
1578 :
1579 950 : if (use_mpi_shared_win) call xmpi_win_fence(XMPI_MODE_NOPRECEDE, epsm1_win, ierr) ! Start the RMA epoch.
1580 :
1581 1872 : select case (approx_type)
1582 :
1583 : case (0)
1584 : ! RPA: \tepsilon = 1 - Vc^{1/2} chi0 Vc^{1/2}
1585 : ! vc_sqrt contains vc^{1/2}(q,G), complex-valued to allow for a possible cutoff.
1586 6273 : do io=1,nomega
1587 6273 : if (omega_distrb(io) == my_rank) then
1588 : !write(std_out,*)"dim_wing",dim_wing
1589 : call rpa_symepsm1(iq_ibz,Vcp,npwe,nI,nJ,chi0(:,:,io),my_nqlwl,dim_wing, &
1590 : chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:), &
1591 5007 : tmp_lf,tmp_nlf,tmp_eelf,xmpi_comm_self)
1592 :
1593 : ! Store results.
1594 10014 : epsm_lf(io,:) = tmp_lf
1595 10014 : epsm_nlf(io,:) = tmp_nlf
1596 10014 : eelf(io,:) = tmp_eelf
1597 : end if
1598 : end do ! nomega
1599 :
1600 : case (1)
1601 : ! Vertex correction from Adiabatic TDDFT. chi_{G1,G2} = [\delta -\chi0 (vc+kxc)]^{-1}_{G1,G3} \chi0_{G3,G2}
1602 12 : ABI_CHECK(Vcp%nqlwl==1,"nqlwl/=1 not coded")
1603 12 : ABI_CHECK(nkxc==1,"nkxc/=1 not coded")
1604 :
1605 : ! Make kxcg_mat(G1,G2) = kxcg(G1-G2) from kxcg defined on the FFT mesh.
1606 48 : ABI_MALLOC_OR_DIE(kxcg_mat,(npwe,npwe), ierr)
1607 :
1608 12 : ierr=0
1609 264 : do ig2=1,npwe
1610 5988 : do ig1=1,npwe
1611 22896 : g1mg2_idx = g2ifft(gvec(:,ig1)-gvec(:,ig2),ngfft)
1612 5976 : if (g1mg2_idx>0) then
1613 5724 : kxcg_mat(ig1,ig2) = kxcg(g1mg2_idx,1)
1614 : else
1615 0 : ierr=ierr+1
1616 0 : kxcg_mat(ig1,ig2) = czero
1617 : end if
1618 : end do
1619 : end do
1620 :
1621 12 : if (ierr/=0) then
1622 : write(msg,'(a,i4,3a)')&
1623 0 : 'Found ',ierr,' G1-G2 vectors falling outside the FFT box. ',ch10,&
1624 0 : 'Enlarge the FFT mesh to get rid of this problem. '
1625 0 : ABI_WARNING(msg)
1626 : end if
1627 :
1628 : !FIXME "recheck TDDFT code and parallel"
1629 12 : ABI_CHECK(nkxc==1,"nkxc/=1 not coded")
1630 264 : do io=1,nomega
1631 264 : if (omega_distrb(io) == my_rank) then
1632 : call atddft_symepsm1(iq_ibz,Vcp,npwe,nI,nJ,chi0(:,:,io),kxcg_mat,option_test,my_nqlwl,dim_wing,omega(io),&
1633 21404 : chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf, xmpi_comm_self)
1634 :
1635 : ! Store results.
1636 504 : epsm_lf(io,:) = tmp_lf
1637 504 : epsm_nlf(io,:) = tmp_nlf
1638 504 : eelf(io,:) = tmp_eelf
1639 : end if
1640 : end do
1641 :
1642 12 : ABI_FREE(kxcg_mat)
1643 264 : do io=1,nomega
1644 252 : write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]'
1645 252 : call wrtout(std_out, msg)
1646 516 : call print_arr([std_out], chi0(:,:,io))
1647 : end do
1648 :
1649 : case (2)
1650 : ! ADA nonlocal vertex correction contained in fxc_ADA
1651 0 : ABI_WARNING('Entered fxc_ADA branch: EXPERIMENTAL!')
1652 : ! Test that argument was passed
1653 0 : if (.not.present(fxc_ADA)) then
1654 0 : ABI_ERROR('make_epsm1_driver was not called with optional argument fxc_ADA')
1655 : end if
1656 0 : ABI_CHECK(Vcp%nqlwl==1,"nqlwl/=1 not coded")
1657 :
1658 0 : do io=1,nomega
1659 0 : if (omega_distrb(io) == my_rank) then
1660 : call atddft_symepsm1(iq_ibz,Vcp,npwe,nI,nJ,chi0(:,:,io),fxc_ADA,option_test,my_nqlwl,dim_wing,omega(io),&
1661 0 : chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,xmpi_comm_self)
1662 :
1663 : ! Store results.
1664 0 : epsm_lf(io,:) = tmp_lf
1665 0 : epsm_nlf(io,:) = tmp_nlf
1666 0 : eelf(io,:) = tmp_eelf
1667 : end if
1668 : end do
1669 :
1670 0 : do io=1,nomega
1671 0 : write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]'
1672 0 : call wrtout(std_out,msg)
1673 0 : call print_arr([std_out], chi0(:,:,io))
1674 : end do
1675 :
1676 : case (4)
1677 : ! Bootstrap vertex kernel by Sharma [[cite:Sharma2011]]
1678 32 : ABI_MALLOC_OR_DIE(vfxc_boot,(npwe*nI,npwe*nJ), ierr)
1679 24 : ABI_MALLOC_OR_DIE(vfxc_boot0,(npwe*nI,npwe*nJ), ierr)
1680 24 : ABI_MALLOC_OR_DIE(chi0_tmp,(npwe*nI,npwe*nJ), ierr)
1681 40 : ABI_MALLOC_OR_DIE(chi0_save,(npwe*nI,npwe*nJ,nomega), ierr)
1682 :
1683 8 : if (iq_ibz==1) then
1684 2 : vc_sqrt => Vcp%vcqlwl_sqrt(:,1) ! Use Coulomb term for q-->0
1685 : else
1686 6 : vc_sqrt => Vcp%vc_sqrt(:,iq_ibz)
1687 : end if
1688 :
1689 19296 : chi0_save = chi0 ! a copy of chi0 (ks)
1690 8 : nstep = 50 ! max iteration steps
1691 8 : alpha = 0.6 ! mixing
1692 8 : chi00_head = chi0(1,1,1)*vc_sqrt(1)**2
1693 3848 : fxc_head = czero; vfxc_boot = czero; chi0_tmp = czero
1694 272 : epsm_lf = czero; epsm_nlf = czero; eelf = zero
1695 8 : write(msg,'(a,2f10.6)') ' -> chi0_dft(head): ', chi00_head
1696 8 : call wrtout(std_out,msg)
1697 : ! loop
1698 8 : conv_err = 0.1
1699 71 : do istep=1, nstep
1700 171181 : chi0 = chi0_save
1701 71 : io=1 ! for now only at omega=0
1702 : call atddft_symepsm1(iq_ibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_boot,0,my_nqlwl,dim_wing,omega(io),&
1703 6151 : chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,xmpi_comm_self)
1704 71 : conv_err = smallest_real
1705 1136 : do ig2=1,npwe*nJ
1706 17111 : do ig1=1,npwe*nI
1707 17040 : conv_err= MAX(conv_err, ABS(chi0(ig1,ig2,1) - chi0_tmp(ig1,ig2)))
1708 : end do
1709 : end do
1710 71 : write(msg,'(a,i4,a,f10.6)') ' => bootstrap itr# ', istep, ', eps^-1 max error: ', conv_err
1711 71 : call wrtout(std_out,msg)
1712 71 : write(msg,'(a,2f10.6)') ' eps^-1(head): ', chi0(1,1,1)
1713 71 : call wrtout(std_out,msg)
1714 71 : write(msg,'(a,2f10.6)') ' v^-1*fxc(head): ', fxc_head
1715 71 : call wrtout(std_out,msg)
1716 71 : if (conv_err <= tol4) exit
1717 : !
1718 15246 : chi0_tmp = chi0(:,:,1)
1719 15246 : vfxc_boot = chi0(:,:,1)/chi00_head
1720 63 : if (istep > 1) then
1721 13310 : vfxc_boot = alpha*vfxc_boot0 + (one-alpha)*vfxc_boot
1722 : end if
1723 15246 : vfxc_boot0 = vfxc_boot
1724 63 : fxc_head = vfxc_boot(1,1)
1725 1016 : do ig1=1,npwe
1726 15183 : vfxc_boot(ig1,:) = vc_sqrt(ig1)*vc_sqrt(:)*vfxc_boot(ig1,:)
1727 : end do
1728 : end do
1729 : ! end loop
1730 8 : if (istep <= nstep) then
1731 8 : write(msg,'(a,i4,a)') ' => bootstrap fxc converged after ', istep, ' iterations'
1732 8 : call wrtout(std_out,msg)
1733 : else
1734 0 : write(msg,'(a,i4,a)') ' -> bootstrap fxc not converged after ', nstep, ' iterations'
1735 0 : ABI_WARNING(msg)
1736 : end if
1737 : !
1738 19288 : chi0 = chi0_save
1739 88 : do io=1,nomega
1740 88 : if (omega_distrb(io) == my_rank) then
1741 : call atddft_symepsm1(iq_ibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_boot,option_test,my_nqlwl,dim_wing,omega(io),&
1742 7760 : chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,xmpi_comm_self)
1743 160 : epsm_lf(io,:) = tmp_lf
1744 160 : epsm_nlf(io,:) = tmp_nlf
1745 160 : eelf(io,:) = tmp_eelf
1746 : end if
1747 : end do
1748 :
1749 8 : ABI_FREE(chi0_tmp)
1750 8 : ABI_FREE(chi0_save)
1751 8 : ABI_FREE(vfxc_boot)
1752 8 : ABI_FREE(vfxc_boot0)
1753 :
1754 88 : do io=1,nomega
1755 80 : write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]'
1756 80 : call wrtout(std_out,msg)
1757 168 : call print_arr([std_out], chi0(:,:,io))
1758 : end do
1759 :
1760 : case (5)
1761 : ! One-shot scalar bootstrap approximation
1762 16 : ABI_MALLOC_OR_DIE(vfxc_boot,(npwe*nI,npwe*nJ), ierr)
1763 20 : ABI_MALLOC_OR_DIE(chi0_save,(npwe*nI,npwe*nJ,nomega), ierr)
1764 :
1765 4 : if (iq_ibz==1) then
1766 1 : vc_sqrt => Vcp%vcqlwl_sqrt(:,1) ! Use Coulomb term for q-->0
1767 : else
1768 3 : vc_sqrt => Vcp%vc_sqrt(:,iq_ibz)
1769 : end if
1770 :
1771 9648 : chi0_save = chi0 ! a copy of chi0
1772 964 : fxc_head = czero; vfxc_boot = czero;
1773 136 : epsm_lf = czero; epsm_nlf = czero; eelf = zero
1774 4 : chi00_head = chi0(1,1,1)*vc_sqrt(1)**2
1775 4 : write(msg,'(a,2f10.6)') ' -> chi0_dft(head): ',chi00_head
1776 4 : call wrtout(std_out,msg)
1777 :
1778 4 : fxc_head = vc_sqrt(1)**2/chi00_head + vc_sqrt(1)**2/chi00_head - vc_sqrt(1)**2
1779 4 : fxc_head = 0.5*fxc_head + 0.5*sqrt(fxc_head**2 - 4.0*vc_sqrt(1)**4/(chi00_head*chi00_head))
1780 4 : vfxc_boot(1,1) = fxc_head
1781 4 : write(msg,'(a,2f10.6)') ' -> v^-1*fxc(head): ',fxc_head/vc_sqrt(1)**2
1782 4 : call wrtout(std_out,msg)
1783 :
1784 9644 : chi0 = chi0_save
1785 44 : do io=1,nomega
1786 44 : if (omega_distrb(io) == my_rank) then
1787 : call atddft_symepsm1(iq_ibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_boot,option_test,my_nqlwl,dim_wing,omega(io),&
1788 3880 : chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,xmpi_comm_self)
1789 80 : epsm_lf(io,:) = tmp_lf
1790 80 : epsm_nlf(io,:) = tmp_nlf
1791 80 : eelf(io,:) = tmp_eelf
1792 : end if
1793 : end do
1794 4 : write(msg,'(a,2f10.6)') ' eps^-1(head): ',chi0(1,1,1)
1795 4 : call wrtout(std_out,msg)
1796 :
1797 4 : ABI_FREE(chi0_save)
1798 4 : ABI_FREE(vfxc_boot)
1799 :
1800 44 : do io=1,nomega
1801 40 : write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]'
1802 40 : call wrtout(std_out,msg)
1803 84 : call print_arr([std_out], chi0(:,:,io))
1804 : end do
1805 :
1806 : case (6)
1807 : ! RPA bootstrap by Rigamonti [[cite:Rigamonti2015]] and Berger [[cite:Berger2015]]
1808 0 : ABI_MALLOC_OR_DIE(vfxc_boot,(npwe*nI,npwe*nJ), ierr)
1809 0 : ABI_MALLOC_OR_DIE(chi0_save,(npwe*nI,npwe*nJ,nomega), ierr)
1810 0 : ABI_MALLOC_OR_DIE(chi0_tmp,(npwe*nI,npwe*nJ), ierr)
1811 :
1812 0 : if (iq_ibz==1) then
1813 0 : vc_sqrt => Vcp%vcqlwl_sqrt(:,1) ! Use Coulomb term for q-->0
1814 : else
1815 0 : vc_sqrt => Vcp%vc_sqrt(:,iq_ibz)
1816 : end if
1817 :
1818 0 : chi0_save = chi0 ! a copy of chi0
1819 0 : fxc_head = czero; vfxc_boot = czero;
1820 0 : epsm_lf = czero; epsm_nlf = czero; eelf = zero
1821 : !chi00_head = chi0(1,1,1)*vc_sqrt(1)**2
1822 0 : write(msg,'(a,2f10.6)') ' -> chi0_dft(head): ',chi00_head
1823 0 : call wrtout(std_out,msg)
1824 :
1825 0 : io = 1 ! static
1826 : call atddft_symepsm1(iq_ibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_boot,0,my_nqlwl,dim_wing,omega(io),&
1827 0 : chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,xmpi_comm_self)
1828 0 : epsm_lf(1,:) = tmp_lf
1829 :
1830 : ! chi(RPA) = chi0 * (1 - chi0 * v_c)^-1
1831 0 : chi0 = chi0_save
1832 0 : do ig2=2,npwe
1833 0 : do ig1=2,npwe
1834 0 : chi0(ig1,ig2,io)=-vc_sqrt(ig1)*chi0(ig1,ig2,io)*vc_sqrt(ig2)
1835 : end do
1836 0 : chi0(ig2,ig2,io)=one+chi0(ig2,ig2,io)
1837 : end do
1838 0 : chi0(1,:,io) = czero; chi0(:,1,io) = czero; chi0(1,1,io) = one
1839 0 : chi0_tmp = chi0(:,:,io)
1840 0 : call xginv(chi0_tmp,npwe,comm=comm)
1841 0 : chi0 = chi0_save
1842 0 : chi0_tmp = MATMUL(chi0(:,:,io), chi0_tmp(:,:)) ! chi(RPA)
1843 0 : do ig1=1,npwe
1844 0 : chi0_tmp(ig1,:) = vc_sqrt(ig1)*vc_sqrt(:)*chi0_tmp(ig1,:)
1845 : end do
1846 : !call xginv(chi0_tmp,npwe,comm=comm) ! chi(RPA)^-1
1847 : !vfxc_boot = chi0_tmp/epsm_lf(1,1)
1848 : !
1849 : !vfxc_boot(1,1) = chi0_tmp(1,1)/epsm_lf(1,1)
1850 0 : vfxc_boot(1,1) = one/chi0_tmp(1,1)/epsm_lf(1,1)
1851 : !@WC: alternatively:
1852 : !chi00_head = chi0(1,1,io)*vc_sqrt(1)**2
1853 : !vfxc_boot(1,1) = one/chi00_head/epsm_lf(1,1)
1854 0 : fxc_head = vfxc_boot(1,1)
1855 0 : do ig1=1,npwe
1856 0 : vfxc_boot(ig1,:) = vc_sqrt(ig1)*vc_sqrt(:)*vfxc_boot(ig1,:)
1857 : end do
1858 0 : write(msg,'(a,2f10.6)') ' -> v^-1*fxc(head): ',fxc_head
1859 0 : call wrtout(std_out,msg)
1860 :
1861 0 : chi0 = chi0_save
1862 0 : do io=1,nomega
1863 0 : if (omega_distrb(io) == my_rank) then
1864 : call atddft_symepsm1(iq_ibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_boot,option_test,my_nqlwl,dim_wing,omega(io),&
1865 0 : chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,xmpi_comm_self)
1866 0 : epsm_lf(io,:) = tmp_lf
1867 0 : epsm_nlf(io,:) = tmp_nlf
1868 0 : eelf(io,:) = tmp_eelf
1869 : end if
1870 : end do
1871 0 : write(msg,'(a,2f10.6)') ' eps^-1(head): ',chi0(1,1,1)
1872 0 : call wrtout(std_out,msg)
1873 :
1874 0 : ABI_FREE(chi0_save)
1875 0 : ABI_FREE(vfxc_boot)
1876 0 : ABI_FREE(chi0_tmp)
1877 :
1878 0 : do io=1,nomega
1879 0 : write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]'
1880 0 : call wrtout(std_out,msg)
1881 0 : call print_arr([std_out], chi0(:,:,io))
1882 : end do
1883 :
1884 : case (7)
1885 : ! LR+ALDA hybrid vertex kernel by Tal
1886 : ! First ALDA
1887 4 : ABI_CHECK(Vcp%nqlwl==1,"nqlwl/=1 not coded")
1888 4 : ABI_CHECK(nkxc==1,"nkxc/=1 not coded")
1889 : ! Make kxcg_mat(G1,G2) = kxcg(G1-G2) from kxcg defined on the FFT mesh.
1890 16 : ABI_MALLOC_OR_DIE(kxcg_mat,(npwe,npwe), ierr)
1891 4 : ierr=0
1892 64 : do ig2=1,npwe
1893 964 : do ig1=1,npwe
1894 3600 : g1mg2_idx = g2ifft(gvec(:,ig1)-gvec(:,ig2),ngfft)
1895 960 : if (g1mg2_idx>0) then
1896 900 : kxcg_mat(ig1,ig2) = kxcg(g1mg2_idx,1)
1897 : else
1898 0 : ierr=ierr+1
1899 0 : kxcg_mat(ig1,ig2) = czero
1900 : end if
1901 : end do
1902 : end do
1903 4 : if (ierr/=0) then
1904 : write(msg,'(a,i4,3a)')&
1905 0 : 'Found ',ierr,' G1-G2 vectors falling outside the FFT box. ',ch10,&
1906 0 : 'Enlarge the FFT mesh to get rid of this problem. '
1907 0 : ABI_WARNING(msg)
1908 : end if
1909 : !FIXME "recheck TDDFT code and parallel"
1910 4 : ABI_CHECK(nkxc==1,"nkxc/=1 not coded")
1911 :
1912 : ! Now LR: (1-Z)*chi0^-1
1913 16 : ABI_MALLOC_OR_DIE(vfxc_lr,(npwe*nI,npwe*nJ), ierr)
1914 12 : ABI_MALLOC_OR_DIE(vfxc_tmp,(npwe*nI,npwe*nJ), ierr)
1915 12 : ABI_MALLOC_OR_DIE(chi0_tmp,(npwe*nI,npwe*nJ), ierr)
1916 :
1917 4 : if (iq_ibz==1) then
1918 1 : vc_sqrt => Vcp%vcqlwl_sqrt(:,1) ! Use Coulomb term for q-->0
1919 : else
1920 3 : vc_sqrt => Vcp%vc_sqrt(:,iq_ibz)
1921 : end if
1922 :
1923 4 : Zr = 0.78
1924 4 : chi00_head = chi0(1,1,1)*vc_sqrt(1)**2
1925 1924 : fxc_head = czero; vfxc_lr = czero; vfxc_tmp = czero
1926 136 : epsm_lf = czero; epsm_nlf = czero; eelf = zero
1927 4 : write(msg,'(a,2f10.6)') ' -> chi0_dft(head): ', chi00_head
1928 4 : call wrtout(std_out,msg)
1929 : !
1930 968 : chi0_tmp = chi0(:,:,1)
1931 4 : call xginv(chi0_tmp,npwe,comm=comm)
1932 968 : vfxc_lr = (one-Zr)*chi0_tmp(:,:)
1933 4 : write(msg,'(a)') ' Constructing LR+ALDA fxc kernel'
1934 4 : call wrtout(std_out,msg)
1935 : !
1936 64 : do ig1=1,npwe
1937 964 : do ig2=1,npwe
1938 3600 : qpg2 = Vcp%qibz(:,iq_ibz) + gvec(:,ig1)
1939 900 : qpg2_nrm = normv(qpg2,gmet,"G")
1940 3600 : qpg2 = Vcp%qibz(:,iq_ibz) + gvec(:,ig2)
1941 900 : qpg2_nrm = SQRT(qpg2_nrm * normv(qpg2,gmet,"G"))
1942 : vfxc_tmp(ig1,ig2) = vfxc_lr(ig1,ig2)*exp(-(qpg2_nrm/k_thfermi(rhor))**2) + &
1943 960 : kxcg_mat(ig1,ig2)*(one - exp(-(qpg2_nrm/k_thfermi(rhor))**2))
1944 : !write(std_out,*) ig1, qpg2_nrm, k_thfermi(rhor), vfxc_lr(ig1,ig1), kxcg_mat(ig1,ig1), vfxc_tmp(ig1,ig1)
1945 : end do
1946 : end do
1947 : !
1948 968 : vfxc_lr = vfxc_tmp
1949 :
1950 44 : do io=1,nomega
1951 44 : if (omega_distrb(io) == my_rank) then
1952 : call atddft_hyb_symepsm1(iq_ibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_lr,kxcg_mat,option_test,my_nqlwl,dim_wing,omega(io),&
1953 3880 : chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,xmpi_comm_self)
1954 80 : epsm_lf(io,:) = tmp_lf
1955 80 : epsm_nlf(io,:) = tmp_nlf
1956 80 : eelf(io,:) = tmp_eelf
1957 : end if
1958 : end do
1959 :
1960 4 : ABI_FREE(kxcg_mat)
1961 4 : ABI_FREE(chi0_tmp)
1962 4 : ABI_FREE(vfxc_lr)
1963 4 : ABI_FREE(vfxc_tmp)
1964 :
1965 44 : do io=1,nomega
1966 40 : write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]'
1967 40 : call wrtout(std_out,msg)
1968 84 : call print_arr([std_out], chi0(:,:,io))
1969 : end do
1970 :
1971 : case default
1972 950 : ABI_BUG(sjoin('Wrong approx_type:',itoa(approx_type)))
1973 : end select
1974 :
1975 950 : if (use_mpi_shared_win) call xmpi_win_fence(XMPI_MODE_NOSUCCEED, epsm1_win, ierr) ! Close the RMA epoch.
1976 :
1977 950 : if (use_MPI .and. .not. use_mpi_shared_win) then
1978 : ! Collect results on each node.
1979 1188 : ABI_MALLOC(buffer_lwing, (size(chi0_lwing,dim=1), size(chi0_lwing, dim=3)))
1980 1188 : ABI_MALLOC(buffer_uwing, (size(chi0_uwing,dim=1), size(chi0_uwing, dim=3)))
1981 :
1982 852 : do io=1,nomega
1983 555 : if (omega_distrb(io)/=my_rank) then
1984 : ! Zero arrays.
1985 967208 : chi0(:,:,io) = czero_gw
1986 344 : if(dim_wing>0) then
1987 5904 : chi0_lwing(:,io,:) = zero
1988 5904 : chi0_uwing(:,io,:) = zero
1989 624 : chi0_head(:,:,io) = czero
1990 : endif
1991 688 : epsm_lf(io,:) = czero
1992 688 : epsm_nlf(io,:) = czero
1993 688 : eelf(io,:) = zero
1994 : end if
1995 :
1996 555 : call xmpi_sum(chi0(:,:,io), comm,ierr)
1997 :
1998 852 : if (dim_wing>0) then
1999 : ! Build contiguous arrays
2000 10562 : buffer_lwing = chi0_lwing(:,io,:)
2001 10562 : buffer_uwing = chi0_uwing(:,io,:)
2002 79 : call xmpi_sum(buffer_lwing,comm,ierr)
2003 79 : call xmpi_sum(buffer_uwing,comm,ierr)
2004 10483 : chi0_lwing(:,io,:) = buffer_lwing
2005 10483 : chi0_uwing(:,io,:) = buffer_uwing
2006 237 : if (size(chi0_head(:,:,io))/= zero) then
2007 79 : call xmpi_sum(chi0_head(:,:,io),comm,ierr)
2008 : end if
2009 : end if
2010 :
2011 : end do ! iw
2012 :
2013 297 : call xmpi_sum(epsm_lf, comm, ierr)
2014 297 : call xmpi_sum(epsm_nlf,comm, ierr)
2015 297 : call xmpi_sum(eelf, comm, ierr)
2016 297 : ABI_FREE(buffer_lwing)
2017 594 : ABI_FREE(buffer_uwing)
2018 : end if
2019 :
2020 : ! Save results in Spectra%, mind the slicing.
2021 4633 : Spectra%emacro_nlf(:,:) = epsm_nlf(1:nor,:)
2022 4633 : Spectra%emacro_lf (:,:) = epsm_lf (1:nor,:)
2023 4633 : Spectra%eelf (:,:) = eelf (1:nor,:)
2024 :
2025 950 : ABI_FREE(epsm_lf)
2026 950 : ABI_FREE(epsm_nlf)
2027 950 : ABI_FREE(eelf)
2028 950 : ABI_FREE(tmp_lf)
2029 950 : ABI_FREE(tmp_nlf)
2030 950 : ABI_FREE(tmp_eelf)
2031 :
2032 950 : call cwtime_report(" make_epsm1_driver", cpu, wall, gflops)
2033 :
2034 2850 : end subroutine make_epsm1_driver
2035 : !!***
2036 :
2037 : !----------------------------------------------------------------------
2038 :
2039 : !!****f* m_screening/rpa_symepsm1
2040 : !! NAME
2041 : !! rpa_symepsm1
2042 : !!
2043 : !! FUNCTION
2044 : !! Calculate RPA $\tilde\epsilon^{-1}$
2045 : !!
2046 : !! Use a special treatment of non-Analytic behavior of heads and wings in reciprocal space
2047 : !! calculating these quantities for different small q-directions specified by the user
2048 : !! (Not yet operative)
2049 : !!
2050 : !! INPUTS
2051 : !! iq_ibz=index of the q-point in the array Vcp%qibz where epsilon^-1 has to be calculated
2052 : !! Vcp<vcoul_t>=Structure gathering data on the Coulombian interaction
2053 : !! npwe=Number of G-vectors in chi0.
2054 : !! nI,nJ=Number of rows/columns in chi0_ij (1,1 in collinear case)
2055 : !! dim_wing=Dimension of the wings (0 or 3 if q-->0)
2056 : !! chi0_head(dim_wing,dim_wing)=Head of of chi0 (only for q-->0)
2057 : !! chi0_lwing(npwe*nI,dim_wing)=Lower wings of chi0 (only for q-->0)
2058 : !! chi0_uwing(npwe*nJ,dim_wing)=Upper wings of chi0 (only for q-->0)
2059 : !! comm=MPI communicator.
2060 : !!
2061 : !! SIDE EFFECTS
2062 : !! chi0(npwe*nI,npwe*nJ): in input the irreducible polarizability,
2063 : !! in output the symmetrized inverse dielectric matrix.
2064 : !!
2065 : !! SOURCE
2066 :
2067 5007 : subroutine rpa_symepsm1(iq_ibz,Vcp,npwe,nI,nJ,chi0,my_nqlwl,dim_wing,chi0_head,chi0_lwing,chi0_uwing,epsm_lf,epsm_nlf,eelf,comm)
2068 :
2069 : !Arguments ------------------------------------
2070 : !scalars
2071 : integer,intent(in) :: iq_ibz,nI,nJ,npwe,dim_wing,my_nqlwl,comm
2072 : type(vcoul_t),target,intent(in) :: Vcp
2073 : !arrays
2074 : complex(gwp),intent(inout) :: chi0(npwe*nI,npwe*nJ)
2075 : complex(dp),intent(inout) :: chi0_lwing(:,:) !(npwe*nI,dim_wing)
2076 : complex(dp),intent(inout) :: chi0_uwing(:,:) !(npwe*nJ,dim_wing)
2077 : complex(dp),intent(inout) :: chi0_head(:,:) !(dim_wing,dim_wing)
2078 : real(dp),intent(out) :: eelf(my_nqlwl)
2079 : complex(dp),intent(out) :: epsm_lf(my_nqlwl),epsm_nlf(my_nqlwl)
2080 :
2081 : !Local variables-------------------------------
2082 : !scalars
2083 : integer,parameter :: master=0,prtvol=0
2084 : integer :: ig1,ig2,iqlwl,my_rank,nprocs
2085 : real(dp) :: ucvol
2086 : logical :: is_qeq0
2087 : !character(len=500) :: msg
2088 : !arrays
2089 : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
2090 5007 : complex(gwp), contiguous, pointer :: vc_sqrt(:)
2091 5007 : complex(gwp),allocatable :: chi0_save(:,:)
2092 : ! *************************************************************************
2093 :
2094 : ABI_UNUSED(chi0_head(1,1))
2095 :
2096 5007 : if (nI/=1.or.nJ/=1) then
2097 0 : ABI_ERROR("nI or nJ=/1 not yet implemented")
2098 : end if
2099 :
2100 5007 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
2101 :
2102 5007 : call metric(gmet,gprimd,-1,rmet,Vcp%rprimd,ucvol)
2103 :
2104 5007 : is_qeq0 = (normv(Vcp%qibz(:,iq_ibz),gmet,'G')<GW_TOLQ0)
2105 5007 : if (is_qeq0) then
2106 790 : ABI_CHECK(iq_ibz==1, "q is 0 but iq_ibz /= 1")
2107 : end if
2108 :
2109 5007 : if (my_nqlwl>1) then
2110 0 : ABI_MALLOC(chi0_save,(npwe*nI,npwe*nJ))
2111 0 : chi0_save = chi0
2112 : end if
2113 : !
2114 : ! Symmetrized RPA epsilon: 1 - Vc^{1/2} chi0 Vc^{1/2}
2115 : ! vc_sqrt contains vc^{1/2}(q, G)
2116 : !
2117 : ! Loop over small q"s (if any) to treat the nonanalytical behavior.
2118 10014 : do iqlwl=my_nqlwl,1,-1
2119 :
2120 5007 : if (my_nqlwl>1) then
2121 0 : chi0(:,:) = chi0_save ! restore pristine polarizability
2122 0 : chi0(:,1) = chi0_lwing(:,iqlwl) ! change the wings
2123 0 : chi0(1,:) = chi0_uwing(:,iqlwl)
2124 : end if
2125 :
2126 5007 : if (iq_ibz==1) then
2127 790 : vc_sqrt => Vcp%vcqlwl_sqrt(:,iqlwl) ! Use Coulomb term for q-->0
2128 : else
2129 4217 : vc_sqrt => Vcp%vc_sqrt(:,iq_ibz)
2130 : end if
2131 :
2132 227716 : do ig2=1,npwe*nJ
2133 17303524 : do ig1=1,npwe*nI
2134 17303524 : chi0(ig1,ig2) = -vc_sqrt(ig1) * chi0(ig1,ig2) * vc_sqrt(ig2)
2135 : end do
2136 227716 : chi0(ig2,ig2) = one + chi0(ig2,ig2)
2137 : end do
2138 : ! chi0, now contains \tepsilon.
2139 :
2140 5007 : epsm_nlf(iqlwl)=chi0(1,1)
2141 :
2142 : if (prtvol > 0) then
2143 : call wrtout(std_out,' Symmetrical epsilon(G,G'') ')
2144 : call print_arr([std_out], chi0)
2145 : end if
2146 : !
2147 : ! === Invert tepsilon and calculate macroscopic dielectric constant ===
2148 : ! * epsm_lf(w) = 1 / epsm1(G=0,Gp=0,w).
2149 : ! * Since G=Gp=0, there is no difference btw symmetrical and not symmetrical.
2150 :
2151 5007 : call xginv(chi0,npwe,comm=comm)
2152 :
2153 5007 : epsm_lf(iqlwl) = one/chi0(1,1)
2154 5007 : eelf(iqlwl) = -AIMAG(chi0(1,1))
2155 :
2156 : if (prtvol > 0) then
2157 : call wrtout(std_out," Symmetrical epsilon^-1(G,G'')")
2158 : call print_arr([std_out], chi0)
2159 : end if
2160 : !
2161 : ! Save wings of e^-1 overwriting input values.
2162 5007 : if (dim_wing>0.and..FALSE.) then
2163 : chi0_lwing(:,iqlwl) = chi0(:,1)
2164 : chi0_uwing(:,iqlwl) = chi0(1,:)
2165 : end if
2166 :
2167 : end do !iqlwl
2168 :
2169 5007 : ABI_SFREE(chi0_save)
2170 :
2171 5007 : end subroutine rpa_symepsm1
2172 : !!***
2173 :
2174 : !----------------------------------------------------------------------
2175 :
2176 : !!****f* m_screening/atddft_symepsm1
2177 : !! NAME
2178 : !! atddft_symepsm1
2179 : !!
2180 : !! FUNCTION
2181 : !! Calculate $\tilde\epsilon^{-1}$ using ALDA within TDDFT
2182 : !!
2183 : !! 2) Use a special treatment of non-Analytic behavior of heads and wings in reciprocal space
2184 : !! calculating these quantities for different small q-directions specified by the user
2185 : !! (Not yet operative)
2186 : !!
2187 : !! Output the electron energy loss function and the macroscopic dielectric function with and
2188 : !! without local field effects (only if non-zero real frequencies are available)
2189 : !!
2190 : !! INPUTS
2191 : !! iq_ibz=index of the q-point in the array Vcp%qibz where epsilon^-1 has to be calculated
2192 : !! npwe=Number of G-vectors in chi0.
2193 : !! nI,nJ=Number of rows/columns in chi0_ij (1,1 in collinear case)
2194 : !! dim_wing=Dimension of the wings (0 or 3 if q-->0)
2195 : !! option_test=Only for TDDFT:
2196 : !! == 0 for TESTPARTICLE ==
2197 : !! == 1 for TESTELECTRON ==
2198 : !! Vcp<vcoul_t>=Structure gathering data on the Coulombian interaction
2199 : !! %nqibz=Number of q-points.
2200 : !! %qibz(3,nqibz)=q-points in the IBZ.
2201 : !! comm=MPI communicator.
2202 : !! chi0_lwing(npwe*nI,dim_wing)=Lower wings of chi0 (only for q-->0)
2203 : !! chi0_uwing(npwe*nJ,dim_wing)=Upper wings of chi0 (only for q-->0)
2204 : !! chi0_head(dim_wing,dim_wing)=Head of of chi0 (only for q-->0)
2205 : !!
2206 : !! OUTPUT
2207 : !!
2208 : !! SIDE EFFECTS
2209 : !! chi0(npwe*nI,npwe*nJ): in input the irreducible polarizability, in output
2210 : !! the symmetrized inverse dielectric matrix.
2211 : !!
2212 : !! SOURCE
2213 :
2214 443 : subroutine atddft_symepsm1(iq_ibz,Vcp,npwe,nI,nJ,chi0,kxcg_mat,option_test,my_nqlwl,dim_wing,omega,&
2215 443 : & chi0_head,chi0_lwing,chi0_uwing,epsm_lf,epsm_nlf,eelf,comm)
2216 :
2217 : !Arguments ------------------------------------
2218 : !scalars
2219 : integer,intent(in) :: iq_ibz,nI,nJ,npwe,dim_wing,my_nqlwl
2220 : integer,intent(in) :: option_test,comm
2221 : type(vcoul_t),target,intent(in) :: Vcp
2222 : !arrays
2223 : complex(gwp),intent(in) :: kxcg_mat(npwe,npwe)
2224 : complex(dp),intent(in) :: omega
2225 : complex(dp),intent(inout) :: chi0_lwing(npwe*nI,dim_wing)
2226 : complex(dp),intent(inout) :: chi0_uwing(npwe*nJ,dim_wing)
2227 : complex(dp),intent(inout) :: chi0_head(dim_wing,dim_wing)
2228 : complex(gwp),intent(inout) :: chi0(npwe*nI,npwe*nJ)
2229 : real(dp),intent(out) :: eelf(my_nqlwl)
2230 : complex(dp),intent(out) :: epsm_lf(my_nqlwl),epsm_nlf(my_nqlwl)
2231 :
2232 : !Local variables-------------------------------
2233 : !scalars
2234 : integer,parameter :: master=0,prtvol=0
2235 : integer :: ig1,ig2,my_rank,nprocs,ierr
2236 : real(dp) :: ucvol
2237 : logical :: is_qeq0
2238 : character(len=500) :: msg
2239 : !arrays
2240 : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
2241 443 : complex(gwp),allocatable :: chitmp(:,:)
2242 443 : complex(gwp), contiguous, pointer :: vc_sqrt(:)
2243 : ! *************************************************************************
2244 :
2245 : ABI_UNUSED(chi0_head(1,1))
2246 : ABI_UNUSED(chi0_lwing(1,1))
2247 : ABI_UNUSED(chi0_uwing(1,1))
2248 :
2249 443 : if (nI/=1.or.nJ/=1) then
2250 0 : ABI_ERROR("nI or nJ=/1 not yet implemented")
2251 : end if
2252 :
2253 443 : ABI_CHECK(Vcp%nqlwl==1,"nqlwl/=1 not coded")
2254 443 : ABI_CHECK(my_nqlwl==1,"my_nqlwl/=1 not coded")
2255 :
2256 443 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
2257 :
2258 443 : call metric(gmet,gprimd,-1,rmet,Vcp%rprimd,ucvol)
2259 :
2260 443 : is_qeq0 = (normv(Vcp%qibz(:,iq_ibz),gmet,'G')<GW_TOLQ0)
2261 :
2262 443 : if (iq_ibz==1) then
2263 : !%vc_sqrt => Vcp%vcqlwl_sqrt(:,iqlwl) ! Use Coulomb term for q-->0
2264 84 : vc_sqrt => Vcp%vcqlwl_sqrt(:,1) ! TODO add treatment of non-Analytic behavior
2265 : else
2266 359 : vc_sqrt => Vcp%vc_sqrt(:,iq_ibz)
2267 : end if
2268 :
2269 443 : write(msg,'(a,f8.2,a)')" chitmp requires: ",npwe**2*gwp*b2Mb," Mb"
2270 1772 : ABI_MALLOC_OR_DIE(chitmp, (npwe,npwe), ierr)
2271 : !
2272 : ! Calculate chi0*fxc.
2273 3607552 : chitmp = MATMUL(chi0,kxcg_mat)
2274 : ! First calculate the NLF contribution
2275 7232 : do ig1=1,npwe
2276 112512 : do ig2=1,npwe
2277 112512 : chitmp(ig1,ig2)=-chitmp(ig1,ig2)
2278 : end do
2279 7232 : chitmp(ig1,ig1)=chitmp(ig1,ig1)+one
2280 : end do
2281 :
2282 443 : call xginv(chitmp,npwe,comm=comm)
2283 :
2284 7102592 : chitmp = MATMUL(chitmp,chi0)
2285 : !if (.not. ABS(REAL(omega))> tol3) call hermitianize(chitmp,"All")
2286 443 : chitmp(1,1)=-vc_sqrt(1)*chitmp(1,1)*vc_sqrt(1)
2287 443 : chitmp(1,1)=chitmp(1,1)+one
2288 :
2289 443 : epsm_nlf(1)=chitmp(1,1)
2290 :
2291 3607109 : chitmp = MATMUL(chi0,kxcg_mat)
2292 : ! * Calculate (1-chi0*Vc-chi0*Kxc) and put it in chitmp.
2293 7232 : do ig1=1,npwe
2294 112512 : do ig2=1,npwe
2295 112512 : chitmp(ig1,ig2)=-chitmp(ig1,ig2)-chi0(ig1,ig2)*vc_sqrt(ig2)**2
2296 : end do
2297 7232 : chitmp(ig1,ig1)=chitmp(ig1,ig1)+one
2298 : end do
2299 :
2300 : ! * Invert (1-chi0*Vc-chi0*Kxc) and Multiply by chi0.
2301 443 : call xginv(chitmp,npwe,comm=comm)
2302 7102592 : chitmp=MATMUL(chitmp,chi0)
2303 :
2304 : ! * Save result, now chi0 contains chi.
2305 112955 : chi0=chitmp
2306 :
2307 874 : select case (option_test)
2308 : case (0)
2309 : ! Symmetrized TESTPARTICLE epsilon^-1
2310 431 : call wrtout(std_out,' Calculating TESTPARTICLE epsilon^-1(G,G") = 1 + Vc*chi')
2311 6896 : do ig1=1,npwe
2312 103440 : chi0(ig1,:)=(vc_sqrt(ig1)*vc_sqrt(:))*chi0(ig1,:)
2313 6896 : chi0(ig1,ig1)=one+chi0(ig1,ig1)
2314 : end do
2315 :
2316 : case (1)
2317 : ! Symmetrized TESTELECTRON epsilon^-1
2318 12 : call wrtout(std_out,' Calculating TESTELECTRON epsilon^-1(G,G") = 1 + (Vc + fxc)*chi')
2319 490560 : chitmp=MATMUL(kxcg_mat,chi0)
2320 :
2321 : ! Perform hermitianization, only valid along the imaginary axis.
2322 12 : if (.not. ABS(REAL(omega))> tol3) call hermitianize(chitmp,"All")
2323 :
2324 336 : do ig1=1,npwe
2325 9072 : chi0(ig1,:)=(vc_sqrt(ig1)*vc_sqrt(:))*chi0(ig1,:)+chitmp(ig1,:)
2326 336 : chi0(ig1,ig1)=one+chi0(ig1,ig1)
2327 : end do
2328 :
2329 : case default
2330 443 : ABI_BUG(sjoin('Wrong option_test:',itoa(option_test)))
2331 : end select
2332 :
2333 443 : ABI_FREE(chitmp)
2334 : !
2335 : ! === chi0 now contains symmetrical epsm1 ===
2336 : ! * Calculate macroscopic dielectric constant epsm_lf(w)=1/epsm1(G=0,Gp=0,w).
2337 443 : epsm_lf(1) = one/chi0(1,1)
2338 443 : eelf (1) = -AIMAG(chi0(1,1))
2339 :
2340 : if (prtvol > 0) then
2341 : write(msg,'(a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at omega',omega*Ha_eV,' [eV]'
2342 : call wrtout(std_out,msg)
2343 : call print_arr([std_out], chi0)
2344 : end if
2345 :
2346 443 : end subroutine atddft_symepsm1
2347 : !!***
2348 :
2349 : !----------------------------------------------------------------------
2350 :
2351 : !!****f* m_screening/atddft_hyb_symepsm1
2352 : !! NAME
2353 : !! atddft_hyb_symepsm1
2354 : !!
2355 : !! FUNCTION
2356 : !! Calculate $\tilde\epsilon^{-1}$ using ALDA within TDDFT
2357 : !!
2358 : !! Based on atddft_symepsm1, modified for the LR+ALDA hybrid kernel
2359 : !!
2360 : !! Output the electron energy loss function and the macroscopic dielectric function with and
2361 : !! without local field effects (only if non-zero real frequencies are available)
2362 : !!
2363 : !! INPUTS
2364 : !! iq_ibz=index of the q-point in the array Vcp%qibz where epsilon^-1 has to be calculated
2365 : !! npwe=Number of G-vectors in chi0.
2366 : !! nI,nJ=Number of rows/columns in chi0_ij (1,1 in collinear case)
2367 : !! dim_wing=Dimension of the wings (0 or 3 if q-->0)
2368 : !! kxcg_mat_sr=Short-range fxc kernel used in the test-electron epsilon^-1
2369 : !! option_test=Only for TDDFT:
2370 : !! == 0 for TESTPARTICLE ==
2371 : !! == 1 for TESTELECTRON ==
2372 : !! Vcp<vcoul_t>=Structure gathering data on the Coulombian interaction
2373 : !! %nqibz=Number of q-points.
2374 : !! %qibz(3,nqibz)=q-points in the IBZ.
2375 : !! comm=MPI communicator.
2376 : !! chi0_lwing(npwe*nI,dim_wing)=Lower wings of chi0 (only for q-->0)
2377 : !! chi0_uwing(npwe*nJ,dim_wing)=Upper wings of chi0 (only for q-->0)
2378 : !! chi0_head(dim_wing,dim_wing)=Head of of chi0 (only for q-->0)
2379 : !!
2380 : !! OUTPUT
2381 : !!
2382 : !! SIDE EFFECTS
2383 : !! chi0(npwe*nI,npwe*nJ): in input the irreducible polarizability, in output
2384 : !! the symmetrized inverse dielectric matrix.
2385 : !!
2386 : !! SOURCE
2387 :
2388 40 : subroutine atddft_hyb_symepsm1(iq_ibz,Vcp,npwe,nI,nJ,chi0,kxcg_mat,kxcg_mat_sr,option_test,my_nqlwl,dim_wing,omega,&
2389 40 : & chi0_head,chi0_lwing,chi0_uwing,epsm_lf,epsm_nlf,eelf,comm)
2390 :
2391 : !Arguments ------------------------------------
2392 : !scalars
2393 : integer,intent(in) :: iq_ibz,nI,nJ,npwe,dim_wing,my_nqlwl
2394 : integer,intent(in) :: option_test,comm
2395 : type(vcoul_t),target,intent(in) :: Vcp
2396 : !arrays
2397 : complex(gwp),intent(in) :: kxcg_mat(npwe,npwe), kxcg_mat_sr(npwe,npwe)
2398 : complex(dp),intent(in) :: omega
2399 : complex(dp),intent(inout) :: chi0_lwing(npwe*nI,dim_wing)
2400 : complex(dp),intent(inout) :: chi0_uwing(npwe*nJ,dim_wing)
2401 : complex(dp),intent(inout) :: chi0_head(dim_wing,dim_wing)
2402 : complex(gwp),intent(inout) :: chi0(npwe*nI,npwe*nJ)
2403 : real(dp),intent(out) :: eelf(my_nqlwl)
2404 : complex(dp),intent(out) :: epsm_lf(my_nqlwl),epsm_nlf(my_nqlwl)
2405 :
2406 : !Local variables-------------------------------
2407 : !scalars
2408 : integer,parameter :: master=0,prtvol=0
2409 : integer :: ig1,ig2,my_rank,nprocs,ierr
2410 : real(dp) :: ucvol, Zr
2411 : logical :: is_qeq0
2412 : character(len=500) :: msg
2413 : !arrays
2414 : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
2415 40 : complex(gwp),allocatable :: chitmp(:,:)
2416 40 : complex(gwp), contiguous, pointer :: vc_sqrt(:)
2417 : ! *************************************************************************
2418 :
2419 : ABI_UNUSED(chi0_head(1,1))
2420 : ABI_UNUSED(chi0_lwing(1,1))
2421 : ABI_UNUSED(chi0_uwing(1,1))
2422 :
2423 40 : if (nI/=1.or.nJ/=1) then
2424 0 : ABI_ERROR("nI or nJ=/1 not yet implemented")
2425 : end if
2426 :
2427 40 : ABI_CHECK(Vcp%nqlwl==1,"nqlwl/=1 not coded")
2428 40 : ABI_CHECK(my_nqlwl==1,"my_nqlwl/=1 not coded")
2429 :
2430 40 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
2431 :
2432 40 : call metric(gmet,gprimd,-1,rmet,Vcp%rprimd,ucvol)
2433 :
2434 40 : is_qeq0 = (normv(Vcp%qibz(:,iq_ibz),gmet,'G')<GW_TOLQ0)
2435 :
2436 40 : if (iq_ibz==1) then
2437 : !%vc_sqrt => Vcp%vcqlwl_sqrt(:,iqlwl) ! Use Coulomb term for q-->0
2438 10 : vc_sqrt => Vcp%vcqlwl_sqrt(:,1) ! TODO add treatment of non-Analytic behavior
2439 : else
2440 30 : vc_sqrt => Vcp%vc_sqrt(:,iq_ibz)
2441 : end if
2442 :
2443 40 : write(msg,'(a,f8.2,a)')" chitmp requires: ",npwe**2*gwp*b2Mb," Mb"
2444 160 : ABI_MALLOC_OR_DIE(chitmp,(npwe,npwe), ierr)
2445 : !
2446 : ! * Calculate chi0*fxc.
2447 289280 : chitmp = MATMUL(chi0,kxcg_mat)
2448 : ! * First calculate the NLF contribution
2449 640 : do ig1=1,npwe
2450 9600 : do ig2=1,npwe
2451 9600 : chitmp(ig1,ig2)=-chitmp(ig1,ig2)
2452 : end do
2453 640 : chitmp(ig1,ig1)=chitmp(ig1,ig1)+one
2454 : end do
2455 :
2456 40 : call xginv(chitmp,npwe,comm=comm)
2457 :
2458 568960 : chitmp = MATMUL(chitmp,chi0)
2459 : !if (.not. ABS(REAL(omega))> tol3) call hermitianize(chitmp,"All")
2460 40 : chitmp(1,1)=-vc_sqrt(1)*chitmp(1,1)*vc_sqrt(1)
2461 40 : chitmp(1,1)=chitmp(1,1)+one
2462 :
2463 40 : epsm_nlf(1)=chitmp(1,1)
2464 :
2465 289240 : chitmp = MATMUL(chi0,kxcg_mat)
2466 : ! * Calculate (1-chi0*Vc-chi0*Kxc) and put it in chitmp.
2467 640 : do ig1=1,npwe
2468 9600 : do ig2=1,npwe
2469 9600 : chitmp(ig1,ig2)=-chitmp(ig1,ig2)-chi0(ig1,ig2)*vc_sqrt(ig2)**2
2470 : end do
2471 640 : chitmp(ig1,ig1)=chitmp(ig1,ig1)+one
2472 : end do
2473 :
2474 : ! * Invert (1-chi0*Vc-chi0*Kxc) and Multiply by chi0.
2475 40 : call xginv(chitmp,npwe,comm=comm)
2476 568960 : chitmp=MATMUL(chitmp,chi0)
2477 :
2478 : ! * Save result, now chi0 contains chi.
2479 9640 : chi0=chitmp
2480 :
2481 40 : select case (option_test)
2482 : case (0)
2483 : ! Symmetrized TESTPARTICLE epsilon^-1
2484 0 : call wrtout(std_out,' Calculating TESTPARTICLE epsilon^-1(G,G") = 1 + Vc*chi')
2485 0 : do ig1=1,npwe
2486 0 : chi0(ig1,:)=(vc_sqrt(ig1)*vc_sqrt(:))*chi0(ig1,:)
2487 0 : chi0(ig1,ig1)=one+chi0(ig1,ig1)
2488 : end do
2489 :
2490 : case (1)
2491 : ! Symmetrized TESTELECTRON epsilon^-1
2492 40 : call wrtout(std_out,' Calculating TESTELECTRON epsilon^-1(G,G") = 1 + Vc*chi + Zr*Kxc_sr*chi')
2493 289280 : chitmp=MATMUL(kxcg_mat_sr,chi0)
2494 40 : Zr = 0.78
2495 :
2496 : ! Perform hermitianization, only valid along the imaginary axis.
2497 40 : if (.not. ABS(REAL(omega))> tol3) call hermitianize(chitmp,"All")
2498 :
2499 640 : do ig1=1,npwe
2500 9600 : chi0(ig1,:)=(vc_sqrt(ig1)*vc_sqrt(:))*chi0(ig1,:)+chitmp(ig1,:)*Zr
2501 640 : chi0(ig1,ig1)=one+chi0(ig1,ig1)
2502 : end do
2503 :
2504 : case default
2505 40 : ABI_BUG(sjoin('Wrong option_test:',itoa(option_test)))
2506 : end select
2507 :
2508 40 : ABI_FREE(chitmp)
2509 : !
2510 : ! === chi0 now contains symmetrical epsm1 ===
2511 : ! * Calculate macroscopic dielectric constant epsm_lf(w)=1/epsm1(G=0,Gp=0,w).
2512 40 : epsm_lf(1) = one/chi0(1,1)
2513 40 : eelf (1) = -AIMAG(chi0(1,1))
2514 :
2515 : if (prtvol > 0) then
2516 : write(msg,'(a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at omega',omega*Ha_eV,' [eV]'
2517 : call wrtout(std_out,msg)
2518 : call print_arr([std_out], chi0)
2519 : end if
2520 :
2521 40 : end subroutine atddft_hyb_symepsm1
2522 : !!***
2523 :
2524 : !----------------------------------------------------------------------
2525 :
2526 : !!****f* m_screening/mkem1_q0
2527 : !! NAME
2528 : !! mkem1_q0
2529 : !!
2530 : !! FUNCTION
2531 : !! This routine construct the microscopic dieletric matrix for q-->0 starting from the heads, wings and the body
2532 : !! of the irreducible polarizability. Afterwards it calculates the symmetrized inverse dieletric matrix
2533 : !! via a block wise inversion thus obtaining the heads and the wings of e^{-1} that can be
2534 : !! used to describe the non-analytic behavior for q-->0.
2535 : !!
2536 : !! INPUTS
2537 : !! npwe=Number of Gs used to describe chi0
2538 : !! nomega=Number of frequencies in chi0.
2539 : !! n1,n2=Factors used to define the same of the chi0 matrix (1,1 if collinear, the typical case)
2540 : !! Cryst<crystal_t>=Structure describing the crystal structure.
2541 : !! Vcp<vcoul_t>=datatypes gathering info on the Coulomb term
2542 : !! gvec(3,npwe)=G-vector for chi0 in reduced coordinates.
2543 : !! comm=MPI communicator
2544 : !!
2545 : !! OUTPUT
2546 : !! eps_head(3,3,nomega)=The macroscopic dieletric tensor in reduced coordinates.
2547 : !! The dieletric matrix along versor \hat q can be obtained with
2548 : !! e(\hat q) = \hat q.eps_head \hat q if all quantities are given in Cartesian coordinates.
2549 : !!
2550 : !! SIDE EFFECTS
2551 : !! chi0(npwe*n1,npwe*n2,nomega)= Input: polarizability. output: inverse dieletric matrix (only the body is used)
2552 : !! chi0_lwing(npwe*n1,nomega,3)
2553 : !! chi0_uwing(npwe*n2,nomega,3) Input: the lower and upper wings of the polarizability
2554 : !! Output: the "lower" and "upper" wings of the inverse dieletric matrix. See notes below.
2555 : !! chi0_head(3,3,nomega)= Input: the polarizability tensor in Cartesian coordinates.
2556 : !! Output: The "head" of the inverse dieletric matrix. See notes below.
2557 : !!
2558 : !! NOTES
2559 : !! Matrix inversion in block form.
2560 : !!
2561 : !! 1 n-1
2562 : !! M = | c u^t| 1 ==> M^{-1} = | 1/k -u^t A^{-1}/k |
2563 : !! | v A | n-1 | -A^{-1} v/k A^{-1} + (A^{-1}v u^t A^{-1})/k |
2564 : !!
2565 : !! where k = c - u^t A^{-1} v
2566 : !!
2567 : !! Let q be a versor in reciprocal space, the symmetrized dielectric matrix with bare coulomb interaction
2568 : !! can be written as
2569 : !!
2570 : !! \tilde\epsilon = | q.Eq q.Wl(G2) | where E_ij = \delta_ij -4\pi chi0_head_ij
2571 : !! | q.Wl(G1) B(G1,G2 | Wl(G1) = -4\pi chi0_lwing(G1)
2572 : !! Wu(G2) = -4\pi chi0_uwing(G1)
2573 : !! therefore, in Cartesian coordinates, we have:
2574 : !!
2575 : !! 1) e^{-1}_{00}(q) = [ q_i q_j (E_{ij} - \sum_{GG'} Wu_i(G)a B_{GG'}^{-1} Wl_j(G')) ]^{-1} = 1/(q.Lq)
2576 : !!
2577 : !! 2) e^{-1}_{0G'}(q) = -e^{-1}_{00}(q) [ \sum_{iG} q_i Wu_i(G)a B_{GG'}^{-1} ] = (q.Su) /(q.Lq)
2578 : !!
2579 : !! 3) e^{-1}_{G0}(q) = -e^{-1}_{00}(q) [ \sum_{iG'} q_i B_{GG'}^{-1} Wl_i(G') ] = (q.Sl) /(q.Lq)
2580 : !!
2581 : !! 4) e^{-1}_{GG'}(q) = B_{GG'}^{-1} +
2582 : !! [ \sum_{ij} q_i q_j ( \sum_T B^{-1}_{GT}^{-1} Wl_i(T)) (\sum_T' Wu_j(T') B^{-1}_{T'G'}^{-1} ] / (q.Lq)
2583 : !!
2584 : !! where Su(G,3) and Sl(G,3) are the "upper" and "lower" wings of the inverse dielectric matrix and
2585 : !! L is the inverse dielectric tensor. Similar equations hold even if vectors and tensors are given in terms
2586 : !! of the reciprocal lattice vectors provided that the metric tensor is taken into account.
2587 : !! The main difference is in the expression for the tensor as only one metric tensor can be
2588 : !! absorbed in the scalar product, the second metric multiplies one of the wings.
2589 : !!
2590 : !! *) The present implementation assumes that no cutoff technique is used in the Coulomb term.
2591 : !!
2592 : !! *) Once the tensor in know it is possible to average the quadratic form on the sphere exactly.
2593 : !! In Cartesian coordinates one obtains.
2594 : !!
2595 : !! \dfrac{1}{4\pi} \int v.Tv d\Omega = Trace(T)/3
2596 : !!
2597 : !! For the inverse dielectric matrix we have to resort to a numerical integration
2598 : !!
2599 : !! SOURCE
2600 :
2601 0 : subroutine mkem1_q0(npwe,n1,n2,nomega,Cryst,Vcp,gvec,chi0_head,chi0_lwing,chi0_uwing,chi0,eps_head,comm)
2602 :
2603 : !Arguments ------------------------------------
2604 : !scalars
2605 : integer,intent(in) :: npwe,nomega,n1,n2,comm
2606 : type(crystal_t),intent(in) :: Cryst
2607 : type(vcoul_t),intent(in) :: Vcp
2608 : !arrays
2609 : integer,intent(in) :: gvec(3,npwe)
2610 : complex(gwp),intent(in) :: chi0(npwe*n1,npwe*n2,nomega)
2611 : complex(dp),intent(inout) :: chi0_lwing(npwe*n1,nomega,3)
2612 : complex(dp),intent(inout) :: chi0_uwing(npwe*n2,nomega,3)
2613 : complex(dp),intent(inout) :: chi0_head(3,3,nomega)
2614 : complex(dp),intent(out) :: eps_head(3,3,nomega)
2615 :
2616 : !Local variables ------------------------------
2617 : !scalars
2618 : integer :: iw,ig,ig1,ig2,idir,jdir
2619 : !arrays
2620 0 : real(dp),allocatable :: modg_inv(:)
2621 0 : complex(dp),allocatable :: eps_lwing(:,:),eps_uwing(:,:),eps_body(:,:),cvec(:)
2622 : !************************************************************************
2623 :
2624 0 : ABI_CHECK(npwe /= 1, "npwe must be >1")
2625 : ABI_UNUSED(comm)
2626 :
2627 : ! Precompute 1/|G|.
2628 0 : ABI_MALLOC(modg_inv,(npwe-1))
2629 0 : do ig=1,npwe-1
2630 0 : modg_inv(ig) = one/normv(gvec(:,ig+1),Cryst%gmet,'G')
2631 : end do
2632 :
2633 0 : ABI_MALLOC(eps_uwing,((npwe-1)*n1,3))
2634 0 : ABI_MALLOC(eps_lwing,((npwe-1)*n2,3))
2635 0 : ABI_MALLOC(eps_body,(npwe-1,npwe-1))
2636 0 : ABI_MALLOC(cvec,(npwe-1))
2637 :
2638 0 : do iw=1,nomega
2639 : !
2640 : ! Head and wings of the symmetrized epsilon.
2641 0 : eps_head(:,:,iw) = -four_pi*chi0_head(:,:,iw)
2642 0 : do idir=1,3
2643 0 : eps_head(idir,idir,iw) = one + eps_head(idir,idir,iw)
2644 0 : eps_lwing(:,idir) = -four_pi * modg_inv * chi0_lwing(2:,iw,idir)
2645 0 : eps_uwing(:,idir) = -four_pi * modg_inv * chi0_uwing(2:,iw,idir)
2646 : !eps_lwing(:,idir) = -chi0_lwing(2:,iw,idir) * SQRT(four_pi) * Vcp%vcqlwl_sqrt(2:npwe,1)
2647 : !eps_uwing(:,idir) = -chi0_uwing(2:,iw,idir) * SQRT(four_pi) * Vcp%vcqlwl_sqrt(2:npwe,1)
2648 : end do
2649 :
2650 0 : write(std_out,*)" espilon head"
2651 0 : call print_arr([std_out], eps_head(:,:,iw))
2652 : !
2653 : ! Construct the body of the symmetrized epsilon then invert it.
2654 0 : do ig2=1,npwe-1
2655 0 : do ig1=1,npwe-1
2656 0 : eps_body(ig1,ig2) = -four_pi * modg_inv(ig1)*chi0(ig1+1,ig2+1,iw )*modg_inv(ig2)
2657 : !eps_body(ig1,ig2) = -Vcp%vcqlwl_sqrt(ig1+1,1)*chi0(ig1+1,ig2+1,iw)* Vcp%vcqlwl_sqrt(ig2+1,1)
2658 : end do
2659 0 : eps_body(ig2,ig2) = one + eps_body(ig2,ig2)
2660 : end do
2661 :
2662 0 : call xginv(eps_body,npwe-1)
2663 : !
2664 : ! Overwrite chi0_head and chi0_wings with the head and the wings of the inverse dielectric matrix.
2665 0 : do jdir=1,3
2666 : !
2667 : ! Head.
2668 0 : cvec=czero
2669 0 : do idir=1,3
2670 0 : cvec = cvec + two_pi*Cryst%gmet(jdir,idir)*MATMUL(eps_body,eps_lwing(:,idir)) ! as we work in reciprocal coords.
2671 : end do
2672 : !cvec = MATMUL(eps_body,eps_lwing(:,jdir))
2673 0 : do idir=1,3
2674 0 : chi0_head(idir,jdir,iw) = eps_head(idir,jdir,iw) - xdotu(npwe-1,eps_uwing(:,idir),1,cvec,1)
2675 : end do
2676 : !
2677 : ! Now the wings.
2678 0 : chi0_uwing(2:,iw,jdir) = -MATMUL(eps_uwing(:,jdir),eps_body)
2679 0 : chi0_lwing(2:,iw,jdir) = -MATMUL(eps_body,eps_lwing(:,jdir))
2680 : !
2681 : end do !jdir
2682 :
2683 0 : call wrtout(std_out, "espilon^1 head after block inversion")
2684 0 : call print_arr([std_out], chi0_head(:,:,iw))
2685 : !
2686 : ! Change the body but do not add the corrections due to the head and the wings.
2687 : ! since they can be obtained on the fly from eps_body and the wings of eps^{-1}.
2688 : !%chi0(2:,2:,iw) = eps_body
2689 : end do !iw
2690 :
2691 0 : ABI_FREE(modg_inv)
2692 0 : ABI_FREE(cvec)
2693 0 : ABI_FREE(eps_lwing)
2694 0 : ABI_FREE(eps_uwing)
2695 0 : ABI_FREE(eps_body)
2696 :
2697 : RETURN
2698 : ABI_UNUSED(Vcp%ng)
2699 :
2700 : end subroutine mkem1_q0
2701 : !!***
2702 :
2703 : !----------------------------------------------------------------------
2704 :
2705 : !!****f* m_screening/lebedev_laikov_int
2706 : !! NAME
2707 : !! lebedev_laikov_int
2708 : !!
2709 : !! FUNCTION
2710 : !!
2711 : !! INPUTS
2712 : !!
2713 : !! OUTPUT
2714 : !!
2715 : !! SOURCE
2716 :
2717 : subroutine lebedev_laikov_int()
2718 :
2719 : !Arguments ------------------------------------
2720 :
2721 : !Local variables-------------------------------
2722 : !scalars
2723 : integer :: on,npts,ii,ll,mm,lmax,leb_idx !ierr,
2724 : real(dp) :: accuracy
2725 : complex(dp) :: ang_int
2726 : !arrays
2727 : real(dp) :: cart_vpt(3) !,real_pars(0)
2728 : real(dp),allocatable :: vx(:),vy(:),vz(:),ww(:)
2729 : complex(dp) :: tensor(3,3),cplx_pars(9)
2730 : complex(dp),allocatable :: ref_func(:),expd_func(:) !tmp_momenta(:)
2731 : ! *************************************************************************
2732 :
2733 : ABI_ERROR("lebedev_laikov_int is still under development")
2734 :
2735 : !tensor=RESHAPE((/4.0,2.0,4.0,0.5,2.1,0.0,5.4,2.1,5.0/),(/3,3/))
2736 : tensor=RESHAPE((/4.0,0.0,0.0,0.0,4.0,0.0,0.0,0.0,5.0/),(/3,3/))
2737 : !tensor=RESHAPE((/1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0/),(/3,3/))
2738 :
2739 : npts=26
2740 : ABI_MALLOC(vx,(npts))
2741 : ABI_MALLOC(vy,(npts))
2742 : ABI_MALLOC(vz,(npts))
2743 : ABI_MALLOC(ww,(npts))
2744 :
2745 : !call LD0026(vx,vy,vz,ww,on)
2746 :
2747 : ang_int=czero
2748 : do ii=1,npts
2749 : cart_vpt = [vx(ii),vy(ii),vz(ii)]
2750 : ang_int = ang_int + ww(ii)*DOT_PRODUCT(cart_vpt,MATMUL(tensor,cart_vpt))
2751 : end do
2752 :
2753 : !write(std_out,*)"quadratic form associated to tensor=",tensor
2754 : write(std_out,*)"on ang_int",on,ang_int
2755 :
2756 : ABI_FREE(vx)
2757 : ABI_FREE(vy)
2758 : ABI_FREE(vz)
2759 : ABI_FREE(ww)
2760 :
2761 : !call init_lebedev_gridset()
2762 : cplx_pars = RESHAPE(tensor,(/9/)); accuracy=tol10
2763 :
2764 : ! This is the function to be expanded evaluated on the lebedev_laikov grid of index leb_idx
2765 : leb_idx=3; npts=lebedev_npts(leb_idx)
2766 : ABI_MALLOC(ref_func,(npts))
2767 : do ii=1,npts
2768 : !cart_vpt = Lgridset(leb_idx)%versor(:,ii)
2769 : ref_func(ii) = one/DOT_PRODUCT(cart_vpt,MATMUL(tensor,cart_vpt))
2770 : end do
2771 :
2772 : ! Calculate the expansion in angular momenta of 1/{q.Tq}.
2773 : ! Only even l-components contribute thanks to the parity of the integrand.
2774 : ! tol6 seems to be an acceptable error, convergence wrt lmax is very slow even for simple tensors.
2775 : ABI_MALLOC(expd_func,(npts))
2776 : expd_func=czero
2777 : lmax=10
2778 : do ll=0,lmax,2
2779 : !allocate(tmp_momenta(-ll:ll))
2780 : do mm=-ll,ll
2781 : ! MG: Commented because it causes problems with the new version of abilint
2782 : !call lebedev_quadrature(ylmstar_over_qTq,(/ll,mm/),real_pars,cplx_pars,ang_int,ierr,accuracy)
2783 : write(std_out,*)ll,mm,ang_int
2784 : !tmp_momenta(mm) = ang_int
2785 : do ii=1,npts
2786 : !cart_vpt = Lgridset(leb_idx)%versor(:,ii)
2787 : expd_func(ii) = expd_func(ii) + four_pi*ang_int*ylmc(ll,mm,cart_vpt)
2788 : end do
2789 : end do
2790 : !deallocate(tmp_momenta)
2791 : write(std_out,*)"Error in angular expansion at l=",ll," is ",MAXVAL(ABS(expd_func-ref_func))
2792 : end do
2793 :
2794 : !BEGINDEBUG
2795 : ! do ii=1,npts
2796 : ! write(777,*)ref_func(ii)
2797 : ! write(778,*)expd_func(ii)
2798 : ! end do
2799 : !ENDDEBUG
2800 :
2801 : ABI_FREE(expd_func)
2802 : ABI_FREE(ref_func)
2803 :
2804 : ABI_ERROR("Exiting from lebedev_laikov_int")
2805 :
2806 : end subroutine lebedev_laikov_int
2807 : !!***
2808 :
2809 : !----------------------------------------------------------------------
2810 :
2811 : !!****f* m_screening/ylmstar_over_qTq
2812 : !! NAME
2813 : !! ylmstar_over_qTq
2814 : !!
2815 : !! FUNCTION
2816 : !! Return Ylm(q)^*/(q.Tq) where q is a versor in Cartesian coordinates.
2817 : !! and Ylm is a complex spherical Harmonics whose index (l,m) are
2818 : !! passed via int_pars(1:2). T is a tensore in Cartesian coordinates
2819 : !! passed via cplx_pars(1:9).
2820 : !!
2821 : !! INPUTS
2822 : !! cart_vers(3)=Cartesian components of the versor
2823 : !! int_pars(1:2)=(l,m) indices in Ylm. l>=0 and m \in [-l,l]
2824 : !! cplx_pars(1:9)=Tensor T in Cartesian coordinates.
2825 : !! real_pars=Not used.
2826 : !!
2827 : !! OUTPUT
2828 : !! Value of Ylm(q)^*/(q.Tq)
2829 : !!
2830 : !! SOURCE
2831 :
2832 : function ylmstar_over_qTq(cart_vers,int_pars,real_pars,cplx_pars)
2833 :
2834 : !Arguments ------------------------------------
2835 : !scalars
2836 : real(dp),intent(in) :: cart_vers(3)
2837 : integer,intent(in) :: int_pars(:)
2838 : real(dp),intent(in) :: real_pars(:)
2839 : complex(dp),intent(in) :: cplx_pars(:)
2840 : complex(dp) :: ylmstar_over_qTq
2841 : !arrays
2842 :
2843 : !Local variables-------------------------------
2844 : !scalars
2845 : integer :: ll,mm
2846 : !arrays
2847 : complex(dp) :: tensor(3,3)
2848 : ! *************************************************************************
2849 :
2850 : tensor = RESHAPE(cplx_pars(1:9),(/3,3/))
2851 : ll = int_pars(1) ! ll starts from zero.
2852 : mm = int_pars(2) ! m \in [-l,l]
2853 :
2854 : ylmstar_over_qTq = CONJG(ylmc(ll,mm,cart_vers))/DOT_PRODUCT(cart_vers,MATMUL(tensor,cart_vers))
2855 :
2856 : RETURN
2857 : ABI_UNUSED(real_pars(1))
2858 :
2859 : end function ylmstar_over_qTq
2860 : !!***
2861 :
2862 : !----------------------------------------------------------------------
2863 :
2864 : !!****f* m_screening/ylmstar_wtq_over_qTq
2865 : !! NAME
2866 : !! ylmstar_wtq_over_qTq
2867 : !!
2868 : !! FUNCTION
2869 : !! Return Ylm(q)^* weight(q)/(q.Tq) where q is a versor in Cartesian coordinates.
2870 : !! Ylm is a complex spherical Harmonics whose index (l,m) are
2871 : !! passed via int_pars(1:2). T is a tensor in Cartesian coordinates
2872 : !! passed via cplx_pars(1:9). weight(q) is the weighting function giving
2873 : !! the length of the vector parallel to versor q that connects the origin
2874 : !! of the lattice to one of the boundaries of the small cell centered at Gamma
2875 : !!
2876 : !! INPUTS
2877 : !! cart_vers(3)=Cartesian components of the versor
2878 : !! int_pars(1:2)=(l,m) indices in Ylm. l>=0 and m \in [-l,l]
2879 : !! cplx_pars(1:9)=Tensor T in Cartesian coordinates.
2880 : !! real_pars(1:9)=The Cartesian vectors defining the small box centered around gamma point
2881 : !! when referred to this vectors the points in the box are given by {(x,y,z) | x,y,z \in [-1,1]}.
2882 : !!
2883 : !! OUTPUT
2884 : !! Value of Ylm(q)^* weigh(q)/(q.Tq)
2885 : !!
2886 : !! SOURCE
2887 :
2888 : function ylmstar_wtq_over_qTq(cart_vers,int_pars,real_pars,cplx_pars)
2889 :
2890 : !Arguments ------------------------------------
2891 : !scalars
2892 : real(dp),intent(in) :: cart_vers(3)
2893 : integer,intent(in) :: int_pars(:)
2894 : real(dp),intent(in) :: real_pars(:)
2895 : complex(dp),intent(in) :: cplx_pars(:)
2896 : complex(dp) :: ylmstar_wtq_over_qTq
2897 : !arrays
2898 :
2899 : !Local variables-------------------------------
2900 : !scalars
2901 : integer :: ll,mm
2902 : real(dp) :: wtq
2903 : !arrays
2904 : real(dp) :: gprimd(3,3),rprimd(3,3),red_vers(3)
2905 : complex(dp) :: tensor(3,3)
2906 : ! *************************************************************************
2907 :
2908 : ABI_ERROR("Work in progress")
2909 : ! box_len has to be tested
2910 :
2911 : gprimd = RESHAPE(real_pars(1:9),(/3,3/))
2912 : red_vers = MATMUL(rprimd,cart_vers)
2913 : wtq = box_len(red_vers,gprimd)
2914 :
2915 : tensor = RESHAPE(cplx_pars(1:9),(/3,3/))
2916 : ll = int_pars(1) ! true ll i.e. not shifted
2917 : mm = int_pars(2)
2918 :
2919 : ylmstar_wtq_over_qTq = CONJG(ylmc(ll,mm,cart_vers))*wtq/DOT_PRODUCT(cart_vers,MATMUL(tensor,cart_vers))
2920 :
2921 : end function ylmstar_wtq_over_qTq
2922 : !!***
2923 :
2924 : !----------------------------------------------------------------------
2925 :
2926 : !!****f* m_screening/mdielf_bechstedt
2927 : !! NAME
2928 : !! mdielf_bechstedt
2929 : !!
2930 : !! FUNCTION
2931 : !! Calculates the model dielectric function for the homogeneous system
2932 : !! as proposed by F. Bechstedt, in Solid State Commun. 84, 765 1992.
2933 : !!
2934 : !! INPUTS
2935 : !! eps_inf=Dielectric constant of the material
2936 : !! qnrm=The modulus of the q-point.
2937 : !! rhor=The local value of the density
2938 : !!
2939 : !! SOURCE
2940 :
2941 140216832 : elemental function mdielf_bechstedt(eps_inf, qnrm, rhor) result(mdielf)
2942 :
2943 : !Arguments ------------------------------------
2944 : !scalars
2945 : real(dp),intent(in) :: eps_inf,qnrm,rhor
2946 : real(dp) :: mdielf
2947 : ! *************************************************************************
2948 :
2949 : mdielf = one + &
2950 140216832 : one / ( one/(eps_inf-one) + (qnrm/k_thfermi(rhor))**2 + (three*qnrm**4)/(four*k_fermi(rhor)**2 * k_thfermi(rhor)**2) )
2951 :
2952 140216832 : end function mdielf_bechstedt
2953 : !!***
2954 :
2955 : !----------------------------------------------------------------------
2956 :
2957 : !!****f* m_screening/screen_mdielf
2958 : !! NAME
2959 : !! screen_mdielf
2960 : !!
2961 : !! FUNCTION
2962 : !! Calculates W_{G,G'}(q,w) for a given q-point in the BZ using a model dielectric function.
2963 : !!
2964 : !! INPUTS
2965 : !! iq_bz=The index of the q-point in the BZ where W(q) is calculated.
2966 : !! npw=Number of plane waves for W
2967 : !! nomega=Number of frequency points.
2968 : !! model_type=Flag defining the model.
2969 : !! eps_inf=Dielectric constant of the material.
2970 : !! Cryst<crystal_t>=Info on the unit cell
2971 : !! Qmesh<kmesh_t>=Info on the set of q-points.
2972 : !! Vcp<vcoul_t datatype>= containing information on the cutoff technique
2973 : !! Gsph<Gsphere>=The G-sphere for W.
2974 : !! nspden=Number of spin density components of the density.
2975 : !! nfft=Number of FFT points on the dense FFT mesh
2976 : !! ngfft(18)=contain all needed information about 3D FFT.
2977 : !! rhor(nfft,nspden)=Electron density in real space (The PAW AE term is included)
2978 : !! which= Set it to "EM1" if the symmetrized inverse dielectric matrix is wanted.
2979 : !! By default the routines returns W.
2980 : !! comm=MPI communicator.
2981 : !!
2982 : !! OUTPUT
2983 : !! w_qbz(npw,npw,nomega)
2984 : !!
2985 : !! NOTES
2986 : !! W_{G1,G2} = 1/2 {
2987 : !! v(q+G1) \int em1(|q+G1|,r) e^{-i(G1-G2).r} dr +
2988 : !! v(q+G2) \int em1(|q+G2|,r) e^{-i(G1-G2).r} dr } / \Omega
2989 : !!
2990 : !! SOURCE
2991 :
2992 631 : subroutine screen_mdielf(iq_bz,npw,nomega,model_type,eps_inf,Cryst,Qmesh,Vcp,Gsph,nspden,nfft,ngfft,rhor,which,w_qbz,comm)
2993 :
2994 : !Arguments ------------------------------------
2995 : !scalars
2996 : integer,intent(in) :: npw,nomega,nfft,nspden,iq_bz,comm,model_type
2997 : real(dp),intent(in) :: eps_inf
2998 : character(len=*),intent(in) :: which
2999 : type(kmesh_t),intent(in) :: Qmesh
3000 : type(crystal_t),intent(in) :: Cryst
3001 : type(vcoul_t),target,intent(in) :: Vcp
3002 : type(gsphere_t),intent(in) :: Gsph
3003 : !arrays
3004 : integer,intent(in) :: ngfft(18)
3005 : real(dp),intent(in) :: rhor(nfft,nspden)
3006 : complex(gwp),intent(out) :: w_qbz(npw,npw,nomega)
3007 :
3008 : !Local variables-------------------------------
3009 : !scalars
3010 : integer,parameter :: tim_fourdp0=0,paral_kgb0=0,cplex1=1
3011 : integer :: my_gstart,my_gstop,iq_ibz,ig,itim_q,isym_q
3012 : integer :: ig1,ig2,g1mg2_fft,iw,ii,ierr,nprocs,isg,ifft !,row ,col
3013 : real(dp) :: qpg2_nrm
3014 : complex(dp) :: ph_mqbzt
3015 : logical :: is_qeq0,isirred
3016 : !character(len=500) :: msg
3017 631 : type(MPI_type) :: MPI_enreg_seq
3018 : !arrays
3019 : integer :: umklp(3)
3020 631 : integer,allocatable :: igfft(:),g1mg2(:,:)
3021 : real(dp) :: qpg2(3),qpt_bz(3)
3022 631 : real(dp),allocatable :: em1_qpg2r(:),fofg(:,:)
3023 631 : complex(gwp),contiguous, pointer :: vc_sqrt_ibz(:)
3024 631 : complex(gwp),allocatable :: vc_qbz(:),ctmp(:,:)
3025 631 : logical,allocatable :: mask(:)
3026 : ! *************************************************************************
3027 :
3028 0 : ABI_CHECK(nomega==1,"screen_mdielf does not support nomega>1")
3029 :
3030 : ! Fake MPI_type for the sequential part.
3031 631 : call initmpi_seq(MPI_enreg_seq)
3032 631 : call MPI_enreg_seq%distribfft%init_seq('c',ngfft(2),ngfft(3),'all')
3033 :
3034 631 : nprocs = xmpi_comm_size(comm)
3035 631 : call xmpi_split_work(npw,comm,my_gstart,my_gstop)
3036 :
3037 631 : call qmesh%get_bz_item(iq_bz,qpt_bz,iq_ibz,isym_q,itim_q,ph_mqbzt,umklp,isirred)
3038 :
3039 : !if (itim_q/=1.or.isym_q/=1.or.ANY(umklp/=0) ) then
3040 : ! ABI_ERROR("Bug in mdielf_bechstedt")
3041 : !end if
3042 : !
3043 : ! Symmetrize Vc in the full BZ.
3044 631 : is_qeq0 = (normv(qpt_bz,Cryst%gmet,'G')<GW_TOLQ0) ! Check if q==0
3045 631 : if (is_qeq0) then
3046 49 : vc_sqrt_ibz => Vcp%vcqlwl_sqrt(:,1) ! Use Coulomb term for q-->0, only first Q is used, shall we average if nqlwl>1?
3047 : else
3048 582 : vc_sqrt_ibz => Vcp%vc_sqrt(:,iq_ibz)
3049 : end if
3050 :
3051 1893 : ABI_MALLOC(vc_qbz,(npw))
3052 33260 : do ig=1,npw
3053 32629 : isg = Gsph%rottb(ig,itim_q,isym_q)
3054 33260 : vc_qbz(isg) = vc_sqrt_ibz(ig)**2
3055 : end do
3056 :
3057 1893 : ABI_MALLOC(igfft,(npw))
3058 1893 : ABI_MALLOC(g1mg2,(3,npw))
3059 1893 : ABI_MALLOC(fofg,(2,nfft))
3060 1893 : ABI_MALLOC(em1_qpg2r,(nfft))
3061 1262 : ABI_MALLOC(mask,(npw))
3062 :
3063 1724402 : w_qbz=czero
3064 31372 : do ig2=my_gstart,my_gstop
3065 : !
3066 : ! Compute the index of G-G2 wave in the FFT grid.
3067 1609860 : do ii=1,npw
3068 6347217 : g1mg2(:,ii) = Gsph%gvec(:,ii) - Gsph%gvec(:,ig2)
3069 : end do
3070 30741 : call kgindex(igfft,g1mg2,mask,MPI_enreg_seq,ngfft,npw)
3071 :
3072 : ! TODO can use zero-padding FFT to speed up the transform.
3073 : !call sphereboundary(gbound,istwfk1,g1mg2,mgfft,npw)
3074 :
3075 : ! Evaluate em1_qpg2r = \int em1(|q+G2|,r) e^{-i(G1-G2).r} dr }.
3076 122964 : qpg2 = qpt_bz + Gsph%gvec(:,ig2)
3077 30741 : qpg2_nrm = normv(qpg2,Cryst%gmet,"G")
3078 :
3079 62113 : do iw=1,nomega
3080 : !
3081 30741 : select case (model_type)
3082 : case (1)
3083 140247573 : do ifft=1,nfft
3084 140247573 : em1_qpg2r(ifft) = one / mdielf_bechstedt(eps_inf,qpg2_nrm,rhor(ifft,1))
3085 : end do
3086 : case default
3087 30741 : ABI_ERROR(sjoin("Unknown model_type:",itoa(model_type)))
3088 : end select
3089 :
3090 30741 : call fourdp(cplex1,fofg,em1_qpg2r,-1,MPI_enreg_seq,nfft,1,ngfft,tim_fourdp0)
3091 : !
3092 : ! Here, unlike the other parts of the code, the unsymmetrized e^{-1} is used.
3093 1640601 : do ig1=1,npw
3094 1579119 : g1mg2_fft = igfft(ig1)
3095 1609860 : w_qbz(ig1,ig2,iw) = DCMPLX(fofg(1,g1mg2_fft), fofg(2,g1mg2_fft)) * vc_qbz(ig2) !/ Cryst%ucvol
3096 : end do
3097 : end do ! iw
3098 : end do ! ig2
3099 :
3100 631 : ABI_FREE(em1_qpg2r)
3101 631 : ABI_FREE(fofg)
3102 631 : ABI_FREE(igfft)
3103 631 : ABI_FREE(g1mg2)
3104 631 : ABI_FREE(mask)
3105 : !
3106 : ! W = 1/2 * (A + A^H)
3107 : ! The MPI sum is done inside the loop to avoid problems with the size of the packet.
3108 2524 : ABI_MALLOC_OR_DIE(ctmp,(npw,npw), ierr)
3109 :
3110 1262 : do iw=1,nomega
3111 : !ctmp = TRANSPOSE(CONJG(w_qbz(:,:,iw)))
3112 1724402 : ctmp = GWPC_CONJG(w_qbz(:,:,iw))
3113 631 : call sqmat_itranspose(npw,ctmp)
3114 1723771 : w_qbz(:,:,iw) = half * (ctmp + w_qbz(:,:,iw))
3115 1262 : call xmpi_sum(w_qbz(:,:,iw),comm,ierr)
3116 : end do
3117 : !
3118 : ! Calculate the symmetrized Em1. W = vc(G1)^{1/2} \tilde Em1 vc(G2)^{1/2} -------------------------
3119 631 : if (toupper(which)=="EM1") then
3120 33260 : do ig=1,npw
3121 32629 : isg = Gsph%rottb(ig,itim_q,isym_q)
3122 33260 : vc_qbz(isg) = vc_sqrt_ibz(ig) ! Workspace storing vc*{1/2}(q_BZ,G).
3123 : end do
3124 :
3125 33260 : do ig2=1,npw
3126 1723771 : do ig1=1,npw
3127 1723140 : ctmp(ig1,ig2) = one / (vc_qbz(ig1) * vc_qbz(ig2))
3128 : end do
3129 : end do
3130 :
3131 1893 : do iw=1,nomega
3132 1724402 : w_qbz(:,:,iw) = w_qbz(:,:,iw) * ctmp(:,:)
3133 : end do
3134 : end if
3135 :
3136 631 : call destroy_mpi_enreg(MPI_enreg_seq)
3137 :
3138 631 : ABI_FREE(vc_qbz)
3139 631 : ABI_SFREE(ctmp)
3140 :
3141 1893 : end subroutine screen_mdielf
3142 : !!***
3143 :
3144 : !----------------------------------------------------------------------
3145 :
3146 : !!****f* m_screening/chi_new
3147 : !! NAME
3148 : !! chi_new
3149 : !!
3150 : !! FUNCTION
3151 : !!
3152 : !! INPUTS
3153 : !!
3154 : !! OUTPUT
3155 : !!
3156 : !! SOURCE
3157 :
3158 0 : type(chi_t) function chi_new(npwe, nomega) result(chi)
3159 :
3160 : !Arguments ------------------------------------
3161 : integer,intent(in) :: npwe,nomega
3162 : ! *************************************************************************
3163 :
3164 0 : chi%nomega = nomega; chi%npwe = npwe
3165 :
3166 0 : ABI_MALLOC(chi%mat, (npwe,npwe,nomega))
3167 :
3168 0 : ABI_MALLOC(chi%head, (3,3,nomega))
3169 0 : ABI_MALLOC(chi%lwing, (npwe, nomega,3))
3170 0 : ABI_MALLOC(chi%uwing, (npwe, nomega,3))
3171 :
3172 0 : end function chi_new
3173 : !!***
3174 :
3175 : !----------------------------------------------------------------------
3176 :
3177 : !!****f* m_screening/chi_free
3178 : !! NAME
3179 : !! chi_free
3180 : !!
3181 : !! FUNCTION
3182 : !!
3183 : !! INPUTS
3184 : !!
3185 : !! OUTPUT
3186 : !!
3187 : !! SOURCE
3188 :
3189 0 : subroutine chi_free(chi)
3190 :
3191 : !Arguments ------------------------------------
3192 : !scalars
3193 : class(chi_t),intent(inout) :: chi
3194 : ! *************************************************************************
3195 :
3196 0 : ABI_SFREE(chi%mat)
3197 0 : ABI_SFREE(chi%head)
3198 0 : ABI_SFREE(chi%lwing)
3199 0 : ABI_SFREE(chi%uwing)
3200 :
3201 0 : end subroutine chi_free
3202 : !!***
3203 :
3204 : !----------------------------------------------------------------------
3205 :
3206 : !!****f* m_screening/lwl_write
3207 : !! NAME
3208 : !! lwl_write
3209 : !!
3210 : !! FUNCTION
3211 : !!
3212 : !! INPUTS
3213 : !!
3214 : !! OUTPUT
3215 : !!
3216 : !! SOURCE
3217 :
3218 0 : subroutine lwl_write(path, cryst, vcp, npwe, nomega, gvec, chi0, chi0_head, chi0_lwing, chi0_uwing, comm)
3219 :
3220 : !Arguments ------------------------------------
3221 : !scalars
3222 : integer,intent(in) :: npwe,nomega,comm
3223 : character(len=*),intent(in) :: path
3224 : type(crystal_t),intent(in) :: cryst
3225 : type(vcoul_t),intent(in) :: Vcp
3226 : !arrays
3227 : integer,intent(in) :: gvec(3,npwe)
3228 : complex(gwp),intent(in) :: chi0(npwe,npwe,nomega)
3229 : complex(dp),intent(inout) :: chi0_head(3,3,nomega),chi0_lwing(npwe,nomega,3),chi0_uwing(npwe,nomega,3)
3230 :
3231 : !Local variables-------------------------------
3232 : !scalars
3233 : integer,parameter :: master=0,prtvol=0
3234 : integer :: iw,ii,iomode,unt,my_rank
3235 : character(len=500) :: msg
3236 : real(dp) :: length
3237 : complex(dp) :: wng(3),em1_00
3238 : ! type(hscr_t),intent(out) :: hscr
3239 : !arrays
3240 0 : complex(dp),allocatable :: wtest(:),eps_head(:,:,:)
3241 : ! *************************************************************************
3242 :
3243 : !if (xmpi_comm_rank(comm) /= master) goto 100
3244 0 : my_rank = xmpi_comm_rank(comm)
3245 :
3246 0 : ABI_MALLOC(wtest,(npwe))
3247 :
3248 : if (prtvol > 0 .and. my_rank == master) then
3249 : call wrtout(std_out, "head of chi0")
3250 : do iw=1,nomega
3251 : call print_arr([std_out], chi0_head(:,:,iw), max_r=3, max_c=3)
3252 : end do
3253 :
3254 : do iw=1,nomega
3255 : call wrtout(std_out, "symmetrized e_00 via tensor")
3256 : wng = MATMUL(chi0_head(:,:,iw),GW_Q0_DEFAULT)
3257 : write(std_out,*) one - vdotw(GW_Q0_DEFAULT,wng,cryst%gmet,"G") * Vcp%vcqlwl_sqrt(1,1)*Vcp%vcqlwl_sqrt(1,1)
3258 :
3259 : call wrtout(std_out, "symmetrized e_0G via tensor")
3260 : do ii=1,npwe
3261 : wng = chi0_uwing(ii,iw,:)
3262 : wtest(ii) = - vdotw(GW_Q0_DEFAULT,wng,cryst%gmet,"G") * Vcp%vcqlwl_sqrt(1,1) * Vcp%vcqlwl_sqrt(ii,1)
3263 : end do
3264 : call print_arr([std_out], wtest,max_r=9)
3265 :
3266 : call wrtout(std_out, "symmetrized e_G0 via tensor")
3267 : do ii=1,npwe
3268 : wng = chi0_lwing(ii,iw,:)
3269 : wtest(ii) = - vdotw(GW_Q0_DEFAULT,wng,cryst%gmet,"G") * Vcp%vcqlwl_sqrt(1,1) * Vcp%vcqlwl_sqrt(ii,1)
3270 : end do
3271 : call print_arr([std_out], wtest,max_r=9)
3272 : end do
3273 : end if
3274 :
3275 : ! Write chi0 data
3276 0 : iomode = IO_MODE_FORTRAN; if (endswith(path, ".nc")) iomode = IO_MODE_ETSF
3277 :
3278 0 : if (my_rank == master) then
3279 0 : if (iomode == IO_MODE_FORTRAN) then
3280 0 : if (open_file(path,msg,newunit=unt,form="unformatted", action="write") /= 0) then
3281 0 : ABI_ERROR(msg)
3282 : end if
3283 : !call epsm1%hscr%io(fform,rdwr,unt,comm,master,iomode)
3284 0 : do iw=1,nomega
3285 0 : write(unt)chi0_head(:,:,iw)
3286 : end do
3287 : !do iw=1,nomega
3288 : ! write(unt)chi0_lwing(:,iw,:)
3289 : !end do
3290 : !do iw=1,nomega
3291 : ! write(unt)chi0_uwing(:,iw,:)
3292 : !end do
3293 :
3294 : else
3295 0 : ABI_ERROR(sjoin("iomode", itoa(iomode), "is not supported"))
3296 : end if
3297 : end if
3298 :
3299 0 : ABI_MALLOC(eps_head,(3,3,nomega))
3300 0 : call mkem1_q0(npwe,1,1,nomega,cryst,Vcp,gvec,chi0_head,chi0_lwing,chi0_uwing,chi0,eps_head,comm)
3301 :
3302 0 : if (my_rank == master) then
3303 0 : if (iomode == IO_MODE_FORTRAN) then
3304 0 : do iw=1,nomega
3305 0 : write(unt)eps_head(:,:,iw)
3306 : end do
3307 : !do iw=1,nomega
3308 : ! write(unt)chi0_lwing(:,iw,:)
3309 : !end do
3310 : !do iw=1,nomega
3311 : ! write(unt)chi0_uwing(:,iw,:)
3312 : !end do
3313 : else
3314 0 : ABI_ERROR(sjoin("iomode:", itoa(iomode), "is not supported"))
3315 : end if
3316 : end if
3317 :
3318 0 : ABI_FREE(eps_head)
3319 :
3320 : if (prtvol > 0 .and. my_rank == master) then
3321 : length = normv(GW_Q0_DEFAULT,cryst%gmet,"G")
3322 :
3323 : do iw=1,nomega
3324 : em1_00 = one / vdotw(GW_Q0_DEFAULT/length, MATMUL(chi0_head(:,:,iw),GW_Q0_DEFAULT/length),cryst%gmet,"G")
3325 : call wrtout(std_out, "e^1_{00} from tensor")
3326 : write(std_out,*) em1_00
3327 :
3328 : call wrtout(std_out, "symmetrized e^-1_0G via tensor")
3329 : do ii=1,npwe
3330 : wng = chi0_uwing(ii,iw,:)
3331 : wtest(ii) = em1_00*vdotw(GW_Q0_DEFAULT/length,wng,cryst%gmet,"G")
3332 : end do
3333 : wtest(1) = em1_00
3334 : call print_arr([std_out], wtest,max_r=9)
3335 :
3336 : call wrtout(std_out, "symmetrized e^-1_G0 via tensor")
3337 : do ii=1,npwe
3338 : wng = chi0_lwing(ii,iw,:)
3339 : wtest(ii) = em1_00*vdotw(GW_Q0_DEFAULT/length,wng,cryst%gmet,"G")
3340 : end do
3341 : wtest(1) = em1_00
3342 : call print_arr([std_out], wtest,max_r=9)
3343 : end do !iw
3344 : end if
3345 :
3346 0 : ABI_FREE(wtest)
3347 :
3348 0 : if (my_rank == master) then
3349 0 : if (iomode == IO_MODE_FORTRAN) then
3350 0 : close(unt)
3351 : else
3352 0 : NCF_CHECK(nf90_close(unt))
3353 : end if
3354 : end if
3355 :
3356 : !100 call xmpi_barrier(comm)
3357 :
3358 0 : end subroutine lwl_write
3359 : !!***
3360 :
3361 : !----------------------------------------------------------------------
3362 :
3363 : !!****f* m_screening/lwl_init
3364 : !! NAME
3365 : !! lwl_init
3366 : !!
3367 : !! FUNCTION
3368 : !!
3369 : !! INPUTS
3370 : !!
3371 : !! OUTPUT
3372 : !!
3373 : !! SOURCE
3374 :
3375 0 : subroutine lwl_init(lwl, path, method, cryst, vcp, npwe, gvec, comm)
3376 :
3377 : !Arguments ------------------------------------
3378 : !scalars
3379 : integer,intent(in) :: comm,method,npwe
3380 : character(len=*),intent(in) :: path
3381 : type(crystal_t),intent(in) :: cryst
3382 : type(vcoul_t),intent(in) :: vcp
3383 : type(lwl_t),intent(out) :: lwl
3384 : !arrays
3385 : integer,intent(in) :: gvec(3,npwe)
3386 :
3387 : !Local variables-------------------------------
3388 : !scalars
3389 : integer,parameter :: master=0
3390 : integer :: iomode,my_rank,nproc,unt
3391 : character(len=500) :: msg
3392 : ! *************************************************************************
3393 :
3394 : ABI_UNUSED((/cryst%natom, gvec(1,1), vcp%ng/))
3395 :
3396 0 : my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
3397 :
3398 0 : lwl%fname = path
3399 0 : lwl%method = method
3400 0 : ABI_CHECK(any(method == [1,2,3]), sjoin("Wrong method:", itoa(method)))
3401 :
3402 0 : iomode = IO_MODE_FORTRAN; if (endswith(path, ".nc")) iomode = IO_MODE_ETSF
3403 :
3404 : ! Only master reads.
3405 0 : if (my_rank == master) then
3406 :
3407 0 : select case (iomode)
3408 : case (IO_MODE_FORTRAN)
3409 0 : if (open_file(path, msg, newunit=unt, action="read", form="unformatted", status="old") /= 0) then
3410 0 : ABI_ERROR(msg)
3411 : end if
3412 :
3413 0 : close(unt)
3414 :
3415 : case default
3416 0 : ABI_ERROR(sjoin("iomode:", itoa(iomode), "is not coded"))
3417 : end select
3418 : end if
3419 :
3420 : ! Broad cast data
3421 : if (nproc > 1) then
3422 : end if
3423 :
3424 0 : end subroutine lwl_init
3425 : !!***
3426 :
3427 : !----------------------------------------------------------------------
3428 :
3429 : !!****f* m_screening/lwl_free
3430 : !! NAME
3431 : !! lwl_free
3432 : !!
3433 : !! FUNCTION
3434 : !!
3435 : !! INPUTS
3436 : !!
3437 : !! OUTPUT
3438 : !!
3439 : !! SOURCE
3440 :
3441 0 : subroutine lwl_free(lwl)
3442 :
3443 : !Arguments ------------------------------------
3444 : !scalars
3445 : type(lwl_t),intent(inout) :: lwl
3446 : ! *************************************************************************
3447 :
3448 0 : ABI_SFREE(lwl%head)
3449 0 : ABI_SFREE(lwl%lwing)
3450 0 : ABI_SFREE(lwl%uwing)
3451 0 : ABI_SFREE(lwl%body)
3452 :
3453 0 : end subroutine lwl_free
3454 : !!***
3455 :
3456 0 : end module m_screening
3457 : !!***
|
