Branch data 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 : iomode__ = iomode
1085 [ - + ]: 4 : if (iomode__ == IO_MODE_MPI) then
1086 : 0 : ABI_WARNING("SUSC files with IO_MODE_MPI is buggy. Using Fortran IO")
1087 : 0 : iomode__ = IO_MODE_FORTRAN
1088 : : end if
1089 : :
1090 [ + - ]: 4 : if (.not. epsm1%use_mpi_shared_win) then
1091 : 4 : call epsm1%shared_comm%set_to_self()
1092 : :
1093 [ - + ]: 4 : if (nprocs > 1) then
1094 : : msg = strcat("- WARNING: Cannot use MPI shared memory as MPI library does not support MPI_WIN_ALLOCATE_SHARED_CPTR", ch10, &
1095 [ # # ]: 0 : "- Memory for epsm1 will increase with nprocs per node!")
1096 : 0 : ABI_WARNING(msg)
1097 : 0 : call wrtout(ab_out, msg)
1098 : : end if
1099 : :
1100 [ - + + - : 24 : ABI_MALLOC_OR_DIE(epsm1%epsm1, (npwe,npwe,epsm1%nomega,epsm1%nqibz), ierr)
- + + - -
+ + - - +
+ - + - +
- + - + -
- + ]
1101 [ + - - - : 4 : call read_screening(in_varname, epsm1%fname, npwe, epsm1%nqibz, epsm1%nomega, epsm1%epsm1, iomode, comm)
- - - - -
- - - - +
- - - - -
- - - ]
1102 : :
1103 : : else
1104 : 0 : call wrtout(std_out, "- HAPPY: Using MPI shared memory, memory for epsm1 won't increase with nprocs per node!")
1105 : 0 : comm__ = comm
1106 : 0 : xcomm = xcomm_from_mpi_int(comm__)
1107 : 0 : epsm1%shared_comm = xcomm%split_type()
1108 : :
1109 : 0 : count = _MOK(2 * npwe) * _MOK(npwe) * _MOK(epsm1%nomega * epsm1%nqibz)
1110 : 0 : call epsm1%shared_comm%allocate_shared_master(count, gwp, xmpi_info_null, void_ptr, epsm1%epsm1_win)
1111 [ # # ]: 0 : call c_f_pointer(void_ptr, epsm1%epsm1, shape=[npwe, npwe, epsm1%nomega, epsm1%nqibz])
1112 : :
1113 : : ! Only one proc in shared_comm reads from file.
1114 : 0 : call xmpi_win_fence(XMPI_MODE_NOPRECEDE, epsm1%epsm1_win, ierr) ! Start the RMA epoch.
1115 : 0 : ABI_CHECK_MPI(ierr, "")
1116 [ # # ]: 0 : if (epsm1%shared_comm%me == 0) then
1117 [ # # # # : 0 : call read_screening(in_varname, epsm1%fname, epsm1%npwe, epsm1%nqibz, epsm1%nomega, epsm1%epsm1, iomode__, xmpi_comm_self)
# # # # #
# # # # #
# # # # #
# # # ]
1118 : : end if
1119 : 0 : call xmpi_win_fence(XMPI_MODE_NOSUCCEED, epsm1%epsm1_win, ierr) ! Close the RMA epoch.
1120 : 0 : ABI_CHECK_MPI(ierr, "")
1121 : 0 : call xcomm%free()
1122 : : end if
1123 : :
1124 [ + + ]: 26 : do iq_ibz=1,epsm1%nqibz
1125 [ + + ]: 22 : is_qeq0=0; if (normv(epsm1%qibz(:,iq_ibz),gmet,'G')<GW_TOLQ0) is_qeq0=1
1126 : :
1127 : 22 : dim_wing=0; if (is_qeq0==1) dim_wing=3 ! FIXME
1128 [ - + + - : 110 : ABI_MALLOC(dummy_lwing,(npwe*epsm1%nI,epsm1%nomega,dim_wing))
- + + - +
+ + - + -
+ + - + -
+ - + ]
1129 [ - + + - : 110 : ABI_MALLOC(dummy_uwing,(npwe*epsm1%nJ,epsm1%nomega,dim_wing))
- + + - +
+ + - + -
+ + - + -
+ - + ]
1130 [ + + + - : 110 : ABI_MALLOC(dummy_head,(dim_wing,dim_wing,epsm1%nomega))
+ + + - -
+ + - + -
+ + - + -
+ - + ]
1131 : :
1132 : 22 : comm__ = comm
1133 [ - + ]: 22 : if (epsm1%use_mpi_shared_win) comm__ = epsm1%shared_comm%value
1134 : :
1135 [ + - ]: 22 : if (approx_type<2 .or. approx_type>3) then
1136 : 22 : ABI_COMMENT('Entering in-core RPA and Kxc branch')
1137 : : call make_epsm1_driver(iq_ibz,dim_wing,npwe,epsm1%nI,epsm1%nJ,epsm1%nomega,epsm1%omega,&
1138 : : approx_type,option_test,Vcp,nfftot,ngfft,nkxc,kxcg,gvec,dummy_head,&
1139 [ - + - - : 22 : dummy_lwing,dummy_uwing,epsm1%epsm1(:,:,:,iq_ibz),spectra,comm__, epsm1%epsm1_win)
- - - - -
- - + - -
- - - - ]
1140 : : else
1141 : 0 : ABI_COMMENT('Entering in-core fxc_ADA branch')
1142 : : call make_epsm1_driver(iq_ibz,dim_wing,npwe,epsm1%nI,epsm1%nJ,epsm1%nomega,epsm1%omega,&
1143 : : approx_type,option_test,Vcp,nfftot,ngfft,nkxc,kxcg,gvec,dummy_head,&
1144 : : dummy_lwing,dummy_uwing,epsm1%epsm1(:,:,:,iq_ibz),spectra,comm__, epsm1%epsm1_win, &
1145 [ # # # # : 0 : fxc_ADA=fxc_ADA(:,:,iq_ibz))
# # # # #
# # # # #
# # # # #
# ]
1146 : : end if
1147 : :
1148 [ - + ]: 22 : ABI_FREE(dummy_lwing)
1149 [ - + ]: 22 : ABI_FREE(dummy_uwing)
1150 [ - + ]: 22 : ABI_FREE(dummy_head)
1151 : :
1152 [ + + ]: 22 : if (is_qeq0==1) then
1153 : 4 : call spectra%repr(msg)
1154 : 4 : call wrtout(units, msg)
1155 : : end if
1156 : :
1157 : 26 : call spectra%free()
1158 : : end do
1159 : :
1160 : 4 : epsm1%id = id_required
1161 [ + + ]: 8 : call epsm1%print([std_out], 0)
1162 : : end if
1163 : : end if
1164 : :
1165 [ - + - + : 498 : end subroutine epsm1_mkdump
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - - -
- - - - -
- - - - -
- - + - +
- + ]
1166 : : !!***
1167 : :
1168 : 132 : subroutine epsm1_malloc_epsm1_qbz(epsm1, npwc, nomega)
1169 : :
1170 : : !Arguments ------------------------------------
1171 : : class(epsm1_t),intent(inout) :: epsm1
1172 : : integer,intent(in) :: npwc, nomega
1173 : :
1174 : : !Local variables-------------------------------
1175 : : integer :: ierr
1176 : : integer(kind=XMPI_ADDRESS_KIND) :: count
1177 : : type(c_ptr) :: void_ptr
1178 : : ! *********************************************************************
1179 : :
1180 [ - + ]: 132 : ABI_CHECK(epsm1%npwe >= npwc, 'Too many G-vectors required')
1181 : :
1182 [ + + ]: 132 : if (.not. epsm1%use_mpi_shared_win) then
1183 [ - + + - : 420 : ABI_MALLOC_OR_DIE(epsm1%epsm1_qbz, (npwc, npwc, nomega), ierr)
- + + - -
+ + - + -
+ - + - +
- - + ]
1184 : : else
1185 : 48 : count = _MOK(2 * npwc) * _MOK(npwc) * _MOK(nomega)
1186 : 48 : call epsm1%shared_comm%allocate_shared_master(count, gwp, xmpi_info_null, void_ptr, epsm1%epsm1_qbz_win)
1187 [ + + ]: 192 : call c_f_pointer(void_ptr, epsm1%epsm1_qbz, shape=[npwc, npwc, nomega])
1188 : : end if
1189 : :
1190 : 132 : end subroutine epsm1_malloc_epsm1_qbz
1191 : : !!***
1192 : :
1193 : : !----------------------------------------------------------------------
1194 : :
1195 : : !!****f* m_screening/epsm1_free_epsm1_qbz
1196 : : !! NAME
1197 : : !! epsm1_free_epsm1_qbz
1198 : : !!
1199 : : !! FUNCTION
1200 : : !! Free the internal buffer %epsm1_qbz
1201 : :
1202 : 460 : subroutine epsm1_free_epsm1_qbz(epsm1)
1203 : :
1204 : : !Arguments ------------------------------------
1205 : : class(epsm1_t),intent(inout) :: epsm1
1206 : :
1207 : : !Local variables-------------------------------
1208 : : integer :: ierr
1209 : : ! *********************************************************************
1210 : :
1211 [ + + ]: 460 : if (.not. epsm1%use_mpi_shared_win) then
1212 [ + + ]: 304 : ABI_SFREE_PTR(epsm1%epsm1_qbz)
1213 : : else
1214 [ + + ]: 156 : if (epsm1%epsm1_qbz_win /= xmpi_undefined) then
1215 : 48 : call xmpi_win_free(epsm1%epsm1_qbz_win, ierr)
1216 : 48 : nullify(epsm1%epsm1_qbz)
1217 : : end if
1218 : : end if
1219 : :
1220 : 460 : end subroutine epsm1_free_epsm1_qbz
1221 : : !!***
1222 : :
1223 : : !----------------------------------------------------------------------
1224 : :
1225 : : !!****f* m_screening/get_epsm1
1226 : : !! NAME
1227 : : !! get_epsm1
1228 : : !!
1229 : : !! FUNCTION
1230 : : !! Work in progress but the main is idea is as follows:
1231 : : !!
1232 : : !! Return the symmetrized inverse dielectric matrix.
1233 : : !! This method implements both in-core and the out-of-core solution
1234 : : !! In the later, epsilon^-1 or chi0 are read from file.
1235 : : !! It is possible to specify options to retrieve (RPA |TDDDT, [TESTCHARGE|TESTPARTICLE]).
1236 : : !! All dimensions are already initialized in the epsm1% object, this method
1237 : : !! should act as a wrapper around rdscr and make_epsm1_driver. A better
1238 : : !! implementation will be done in the following once the coding of file handlers is completed.
1239 : : !!
1240 : : !! INPUTS
1241 : : !! Vcp<vcoul_t>=Structure gathering data on the Coulombian interaction
1242 : : !! iqibzA[optional]=Index of the q-point to be read from file (only for out-of-memory solutions)
1243 : : !! iomode=option definig the file format.
1244 : : !! option_test
1245 : : !! comm=MPI communicator.
1246 : : !!
1247 : : !! OUTPUT
1248 : : !! epsm1%epsm1
1249 : : !!
1250 : : !! TODO
1251 : : !! Remove this routine. Now everything should be done with epsm1_mkdump
1252 : : !!
1253 : : !! SOURCE
1254 : :
1255 : 290 : subroutine get_epsm1(epsm1,Vcp,approx_type,option_test,iomode,comm,iqibzA)
1256 : :
1257 : : !Arguments ------------------------------------
1258 : : !scalars
1259 : : class(epsm1_t),intent(inout) :: epsm1
1260 : : integer,intent(in) :: iomode,option_test,approx_type,comm
1261 : : integer,optional,intent(in) :: iqibzA
1262 : : type(vcoul_t),intent(in) :: Vcp
1263 : :
1264 : : !Local variables-------------------------------
1265 : : !scalars
1266 : : integer :: my_approx_type,my_option_test,ng,ierr
1267 : : ! *********************************************************************
1268 : :
1269 : 290 : my_approx_type = approx_type; my_option_test = option_test
1270 : :
1271 : : ! Vcp not yet used.
1272 : 290 : ng = Vcp%ng
1273 : :
1274 : 580 : select case (epsm1%mqmem)
1275 : : case (0)
1276 : : ! Out-of-core solution
1277 [ + + ]: 290 : ABI_SFREE_PTR(epsm1%epsm1)
1278 [ - + + - : 1450 : ABI_MALLOC_OR_DIE(epsm1%epsm1,(epsm1%npwe,epsm1%npwe,epsm1%nomega,1), ierr)
- + + - -
+ + - + -
+ - + - +
- - + ]
1279 : :
1280 [ + - - - : 290 : call read_screening(em1_ncname,epsm1%fname,epsm1%npwe,epsm1%nqibz,epsm1%nomega,epsm1%epsm1,iomode,comm,iqiA=iqibzA)
- - - - -
- - - - +
- - - - -
- - - ]
1281 : :
1282 [ - + ]: 290 : if (epsm1%id == 4) then
1283 : : ! If q-slice of epsilon^-1 has been read then return
1284 : : !call epsm1%print()
1285 : : return
1286 : : else
1287 : 0 : ABI_ERROR(sjoin('Wrong epsm1%ID', itoa(epsm1%id)))
1288 : : end if
1289 : :
1290 : : case default
1291 : : ! In-core solution.
1292 [ + - ]: 290 : ABI_ERROR("you should not be here")
1293 : : end select
1294 : :
1295 : : end subroutine get_epsm1
1296 : : !!***
1297 : :
1298 : : !----------------------------------------------------------------------
1299 : :
1300 : : !!****f* m_screening/decompose_epsm1
1301 : : !! NAME
1302 : : !! decompose_epsm1
1303 : : !!
1304 : : !! FUNCTION
1305 : : !! Decompose the complex symmetrized dielectric
1306 : : !!
1307 : : !! INPUTS
1308 : : !!
1309 : : !! OUTPUT
1310 : : !!
1311 : : !! SOURCE
1312 : :
1313 : 0 : subroutine decompose_epsm1(epsm1, iq_ibz, eigs)
1314 : :
1315 : : !Arguments ------------------------------------
1316 : : !scalars
1317 : : class(epsm1_t),intent(in) :: epsm1
1318 : : integer,intent(in) :: iq_ibz
1319 : : !arrays
1320 : : complex(dp),intent(out) :: eigs(epsm1%npwe,epsm1%nomega)
1321 : :
1322 : : !Local variables-------------------------------
1323 : : !scalars
1324 : : integer :: info,lwork,iw,negw,ig1,ig2,idx,sdim,npwe,ierr
1325 : : character(len=500) :: msg
1326 : : !arrays
1327 : 0 : real(dp),allocatable :: ww(:),rwork(:)
1328 : 0 : complex(dp),allocatable :: work(:),Adpp(:),eigvec(:,:),Afull(:,:),vs(:,:),wwc(:)
1329 : 0 : logical,allocatable :: bwork(:)
1330 : : logical :: sortcplx !BUG in abilint
1331 : : ! *********************************************************************
1332 : :
1333 [ # # ]: 0 : ABI_CHECK(epsm1%mqmem /= 0, 'mqmem==0 not implemented')
1334 : :
1335 : 0 : npwe = epsm1%npwe
1336 : :
1337 [ # # ]: 0 : do iw=1,epsm1%nomega
1338 : :
1339 [ # # ]: 0 : if (ABS(REAL(epsm1%omega(iw)))>0.00001) then
1340 : : ! Eigenvalues for a generic complex matrix
1341 : 0 : lwork=4*2*npwe
1342 [ # # # # : 0 : ABI_MALLOC(wwc,(npwe))
# # # # #
# # # ]
1343 [ # # # # : 0 : ABI_MALLOC(work,(lwork))
# # # # #
# # # ]
1344 [ # # # # : 0 : ABI_MALLOC(rwork,(npwe))
# # # # #
# # # ]
1345 [ # # # # : 0 : ABI_MALLOC(bwork,(npwe))
# # # # #
# # # ]
1346 [ # # # # : 0 : ABI_MALLOC(vs,(npwe,npwe))
# # # # #
# # # # #
# # ]
1347 [ # # # # : 0 : ABI_MALLOC(Afull,(npwe,npwe))
# # # # #
# # # ]
1348 : :
1349 [ # # # # : 0 : Afull=epsm1%epsm1(:,:,iw,iq_ibz)
# # # # #
# # # ]
1350 : :
1351 : : !for the moment no sort, maybe here I should sort using the real part?
1352 : 0 : call ZGEES('V','N',sortcplx,npwe,Afull,npwe,sdim,wwc,vs,npwe,work,lwork,rwork,bwork,info)
1353 [ # # ]: 0 : if (info/=0) then
1354 : 0 : ABI_ERROR(sjoin("ZGEES returned info:",itoa(info)))
1355 : : end if
1356 : :
1357 [ # # ]: 0 : eigs(:,iw)=wwc(:)
1358 : :
1359 [ # # ]: 0 : ABI_FREE(wwc)
1360 [ # # ]: 0 : ABI_FREE(work)
1361 [ # # ]: 0 : ABI_FREE(rwork)
1362 [ # # ]: 0 : ABI_FREE(bwork)
1363 [ # # ]: 0 : ABI_FREE(vs)
1364 [ # # ]: 0 : ABI_FREE(Afull)
1365 : :
1366 : : else
1367 : : ! Hermitian version.
1368 : 0 : lwork=2*npwe-1
1369 [ # # # # : 0 : ABI_MALLOC(ww,(npwe))
# # # # #
# # # ]
1370 [ # # # # : 0 : ABI_MALLOC(work,(lwork))
# # # # #
# # # ]
1371 [ # # # # : 0 : ABI_MALLOC(rwork,(3*npwe-2))
# # # # #
# # # #
# ]
1372 [ # # # # : 0 : ABI_MALLOC(eigvec,(npwe,npwe))
# # # # #
# # # # #
# # ]
1373 [ # # # # : 0 : ABI_MALLOC_OR_DIE(Adpp,(npwe*(npwe+1)/2), ierr)
# # # # #
# # # # #
# # ]
1374 : :
1375 : : idx=0 ! Pack the matrix
1376 [ # # ]: 0 : do ig2=1,npwe
1377 [ # # ]: 0 : do ig1=1,ig2
1378 : 0 : idx=idx+1
1379 : 0 : Adpp(idx)=epsm1%epsm1(ig1,ig2,iw,iq_ibz)
1380 : : end do
1381 : : end do
1382 : :
1383 : : ! For the moment we require also the eigenvectors.
1384 : 0 : call ZHPEV('V','U',npwe,Adpp,ww,eigvec,npwe,work,rwork,info)
1385 [ # # ]: 0 : if (info/=0) then
1386 : 0 : ABI_ERROR(sjoin('ZHPEV returned info=', itoa(info)))
1387 : : end if
1388 : :
1389 [ # # # # ]: 0 : negw=(COUNT((REAL(ww)<tol6)))
1390 [ # # ]: 0 : if (negw/=0) then
1391 : : write(msg,'(a,i5,a,i3,a,f8.4)')&
1392 [ # # # # : 0 : 'Found negative eigenvalues. No. ',negw,' at iq_ibz= ',iq_ibz,' minval= ',MINVAL(REAL(ww))
# # # # ]
1393 : 0 : ABI_WARNING(msg)
1394 : : end if
1395 : :
1396 [ # # ]: 0 : eigs(:,iw)=ww(:)
1397 : :
1398 [ # # ]: 0 : ABI_FREE(ww)
1399 [ # # ]: 0 : ABI_FREE(work)
1400 [ # # ]: 0 : ABI_FREE(rwork)
1401 [ # # ]: 0 : ABI_FREE(eigvec)
1402 [ # # ]: 0 : ABI_FREE(Adpp)
1403 : : end if
1404 : : end do !iw
1405 : :
1406 : : ! contains
1407 : : ! function sortcplx(carg) result(res)
1408 : : ! implicit none
1409 : : ! complex(dp),intent(in) :: carg
1410 : : ! logical :: res
1411 : : ! res=.TRUE.
1412 : : ! end function sortcplx
1413 : :
1414 : 0 : end subroutine decompose_epsm1
1415 : : !!***
1416 : :
1417 : : !----------------------------------------------------------------------
1418 : :
1419 : : !!****f* m_screening/make_epsm1_driver
1420 : : !! NAME
1421 : : !! make_epsm1_driver
1422 : : !!
1423 : : !! FUNCTION
1424 : : !! Driver routine to calculate the inverse symmetrical dielectric matrix starting
1425 : : !! from the irreducible polarizability. The routine considers a single q-point, and
1426 : : !! performs the following tasks:
1427 : : !!
1428 : : !! 1) Calculate $\tilde\epsilon^{-1}$ using different approximations:
1429 : : !! * RPA
1430 : : !! * ALDA within TDDFT
1431 : : !!
1432 : : !! 2) Use a special treatment of non-Analytic behavior of heads and wings in reciprocal space
1433 : : !! calculating these quantities for different small q-directions specified by the user
1434 : : !! (Not yet operative)
1435 : : !!
1436 : : !! 3) Output the electron energy loss function and the macroscopic dielectric function with and
1437 : : !! without local field effects (only if non-zero real frequencies are available)
1438 : : !!
1439 : : !! INPUTS
1440 : : !! iq_ibz=index of the q-point in the array Vcp%qibz where epsilon^-1 has to be calculated
1441 : : !! dim_wing=Dimension of the wings (0 or 3 if q-->0)
1442 : : !! npwe=Number of G-vectors in chi0.
1443 : : !! nI,nJ=Number of rows/columns in chi0_ij (1,1 if collinear case)
1444 : : !! nomega=Number of frequencies.
1445 : : !! omega(nomega)=Frequencines in Hartree
1446 : : !! approx_type=Integer flag defining the type of approximation
1447 : : !! == 0 for RPA ==
1448 : : !! == 1 for TDDFT ==
1449 : : !! option_test=Only for TDDFT:
1450 : : !! == 0 for TESTPARTICLE ==
1451 : : !! == 1 for TESTELECTRON ==
1452 : : !! Vcp<vcoul_t>=Structure gathering data on the Coulombian interaction
1453 : : !! nfftot=Total number of points in the FFT mesh.
1454 : : !! ngfft(18)=Info on the FFT mesh.
1455 : : !! nkxc=Dimension of the kernel in reciprocal space. 0 if kernel is not needed
1456 : : !! kxcg(nfftot,nkxc)=TDDFT kernel in reciprocal space on the FFT mesh. Used only if approx_type==1
1457 : : !! gvec(3,npwe)=G-vectors
1458 : : !! comm=MPI communicator.
1459 : : !! chi0_lwing(npwe*nI,nomega,dim_wing)=Lower wings of chi0 (only for q-->0)
1460 : : !! chi0_uwing(npwe*nJ,nomega,dim_wing)=Upper wings of chi0 (only for q-->0)
1461 : : !! chi0_head(dim_wing,dim_wing,nomega)=Head of of chi0 (only for q-->0)
1462 : : !!
1463 : : !! OUTPUT
1464 : : !! spectra<spectra_t>Object containing e_macro(w) and EELS(w)
1465 : : !!
1466 : : !! SIDE EFFECTS
1467 : : !! chi0(npwe*nI,npwe*nJ,nomega): in input the irreducible polarizability, in output
1468 : : !! the symmetrized inverse dielectric matrix.
1469 : : !!
1470 : : !! SOURCE
1471 : :
1472 [ - + - + : 950 : subroutine make_epsm1_driver(iq_ibz, dim_wing, npwe, nI, nJ, nomega, omega,&
- + - + ]
1473 : 950 : approx_type, option_test, Vcp, nfftot, ngfft, nkxc, kxcg, gvec, &
1474 [ - + - + : 950 : chi0_head, chi0_lwing, chi0_uwing, chi0, spectra, comm, epsm1_win, &
- + ]
1475 : 950 : fxc_ADA, rhor) ! optional argument
1476 : :
1477 : : !Arguments ------------------------------------
1478 : : !scalars
1479 : : integer,intent(in) :: iq_ibz,nI,nJ,npwe,nomega,dim_wing,approx_type,option_test,nkxc,nfftot,comm, epsm1_win
1480 : : real(dp),intent(in),optional :: rhor
1481 : : type(vcoul_t),target,intent(in) :: Vcp
1482 : : type(spectra_t),intent(out) :: Spectra
1483 : : !arrays
1484 : : integer,intent(in) :: ngfft(18),gvec(3,npwe)
1485 : : complex(gwp),intent(in) :: kxcg(nfftot,nkxc)
1486 : : complex(dp),intent(in) :: omega(nomega)
1487 : : complex(dp),intent(inout) :: chi0_lwing(:,:,:) !(npwe*nI,nomega,dim_wing)
1488 : : complex(dp),intent(inout) :: chi0_uwing(:,:,:) !(npwe*nJ,nomega,dim_wing)
1489 : : complex(dp),intent(inout) :: chi0_head(:,:,:) !(dim_wing,dim_wing,nomega)
1490 : : complex(gwp),intent(inout) :: chi0(npwe*nI,npwe*nJ,nomega)
1491 : : complex(gwp),intent(in),optional :: fxc_ADA(npwe*nI,npwe*nJ)
1492 : :
1493 : : !Local variables-------------------------------
1494 : : !scalars
1495 : : integer,parameter :: master=0
1496 : : integer :: i1,i2,ig1,ig2,io,ierr,irank,my_nqlwl !iqlwl
1497 : : integer :: nor,my_rank,nprocs,g1mg2_idx
1498 : : real(dp) :: ucvol
1499 : : logical :: is_qeq0, use_MPI, use_mpi_shared_win
1500 : : character(len=500) :: msg
1501 : : !arrays
1502 : 1900 : integer :: omega_distrb(nomega)
1503 : 950 : integer,allocatable :: istart(:),istop(:)
1504 : : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
1505 : 950 : real(dp),allocatable :: eelf(:,:),tmp_eelf(:)
1506 : 950 : complex(dp),allocatable :: epsm_lf(:,:),epsm_nlf(:,:),tmp_lf(:),tmp_nlf(:)
1507 : 950 : complex(dp),allocatable :: buffer_lwing(:,:),buffer_uwing(:,:)
1508 : 950 : complex(gwp),allocatable :: kxcg_mat(:,:)
1509 : :
1510 : : !bootstrap and LR
1511 : : integer :: istep,nstep
1512 : : real(dp) :: conv_err, alpha, Zr, qpg2(3), qpg2_nrm, cpu, wall, gflops
1513 : : real(gwp) :: chi00_head, fxc_head
1514 : 950 : complex(gwp),allocatable :: vfxc_boot(:,:), vfxc_boot0(:,:), vfxc_lr(:,:), vfxc_tmp(:,:), chi0_tmp(:,:), chi0_save(:,:,:)
1515 : 950 : complex(gwp), contiguous, pointer :: vc_sqrt(:)
1516 : : ! *************************************************************************
1517 : :
1518 [ + - - + ]: 950 : if (nI/=1.or.nJ/=1) then
1519 : 0 : ABI_ERROR("nI or nJ=/1 not yet implemented")
1520 : : end if
1521 : :
1522 : 950 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
1523 : 950 : call cwtime(cpu, wall, gflops, "start")
1524 : :
1525 : 950 : call metric(gmet,gprimd,-1,rmet,Vcp%rprimd,ucvol)
1526 : :
1527 : 950 : is_qeq0 = normv(Vcp%qibz(:,iq_ibz),gmet,'G') < GW_TOLQ0
1528 : :
1529 [ + + ]: 6713 : omega_distrb = my_rank
1530 : 950 : use_mpi_shared_win = (epsm1_win /= xmpi_undefined)
1531 : 950 : use_MPI = .FALSE.
1532 : 950 : use_MPI = nprocs >= nomega ! Parallelism is not used
1533 : :
1534 : : ! FIXME: MPI mode is temporarily disabled here because we need to know if
1535 : : ! screening is allocated in shared memory or not.
1536 : : ! Perhaps now it makes mores sense to use Scalapack/ELPA instead of parallelizing the loop over frequencies
1537 : : !use_MPI = .FALSE.
1538 : :
1539 [ + + ]: 950 : if (use_MPI) then
1540 : : ! Initialize distribution table for frequencies.
1541 [ - + + - : 891 : ABI_MALLOC(istart,(nprocs))
+ - + - -
+ - + ]
1542 [ - + + - : 594 : ABI_MALLOC(istop,(nprocs))
- + - + ]
1543 : 297 : call xmpi_split_work2_i4b(nomega,nprocs,istart,istop)
1544 [ + + ]: 852 : omega_distrb(:)=xmpi_undefined_rank
1545 [ + + ]: 1196 : do irank=0,nprocs-1
1546 : 899 : i1 = istart(irank+1)
1547 : 899 : i2 = istop (irank+1)
1548 [ + + + + ]: 1751 : if (i1<=i2) omega_distrb(i1:i2) = irank
1549 : : end do
1550 [ - + ]: 297 : ABI_FREE(istart)
1551 [ - + ]: 297 : ABI_FREE(istop)
1552 : : end if
1553 : :
1554 : : ! Initialize container for spectral results
1555 [ + + ]: 3677 : do nor=1,nomega
1556 [ + + ]: 3677 : if (ABS(AIMAG(omega(nor)))>1.e-3) EXIT
1557 : : end do
1558 [ + + ]: 950 : nor=nor-1; if (nor==0) nor = 1 ! only imag !?
1559 : :
1560 [ + + ]: 950 : if (dim_wing==3) then
1561 : 151 : call wrtout(std_out,' Analyzing long wavelength limit for several q')
1562 [ + - + + : 564 : call spectra_init(Spectra,nor,REAL(omega(1:nor)),Vcp%nqlwl,Vcp%qlwl)
- + ]
1563 : 151 : my_nqlwl = 1
1564 : : !my_nqlwl = dim_wing ! TODO
1565 : : !ABI_CHECK(dim_wing==SIZE(Vcp%vcqlwl_sqrt,DIM=2),"WRONG DIMS")
1566 : : else
1567 [ + - + + : 3119 : call spectra_init(Spectra,nor,REAL(omega(1:nor)),1,Vcp%qibz(:,iq_ibz))
- + ]
1568 : 799 : my_nqlwl = 1
1569 : : end if
1570 : : !
1571 : : ! NOTE: all processors have to perform this operation in order to have the
1572 : : ! epsm1 matrix when performing a sigma calculation starting with the file _SUS
1573 : : !
1574 : : ! Temporary arrays to store spectra.
1575 [ - + + - : 3800 : ABI_MALLOC(epsm_lf,(nomega,my_nqlwl))
- + + - +
- + - - +
- + ]
1576 [ - + + - : 2850 : ABI_MALLOC(epsm_nlf,(nomega,my_nqlwl))
- + + - -
+ - + ]
1577 [ - + + - : 3800 : ABI_MALLOC(eelf,(nomega,my_nqlwl))
- + + - +
- + - - +
- + ]
1578 [ + + + + : 21089 : epsm_lf=czero; epsm_nlf=czero; eelf=zero
+ + + + +
+ + + ]
1579 : :
1580 : : ! Temporary arrays used to store output results.
1581 [ - + + - : 2850 : ABI_MALLOC(tmp_lf, (my_nqlwl))
+ - + - -
+ - + ]
1582 [ - + + - : 1900 : ABI_MALLOC(tmp_nlf, (my_nqlwl))
- + - + ]
1583 [ - + + - : 2850 : ABI_MALLOC(tmp_eelf, (my_nqlwl))
+ - + - -
+ - + ]
1584 : :
1585 [ - + ]: 950 : if (use_mpi_shared_win) call xmpi_win_fence(XMPI_MODE_NOPRECEDE, epsm1_win, ierr) ! Start the RMA epoch.
1586 : :
1587 : 1872 : select case (approx_type)
1588 : :
1589 : : case (0)
1590 : : ! RPA: \tepsilon = 1 - Vc^{1/2} chi0 Vc^{1/2}
1591 : : ! vc_sqrt contains vc^{1/2}(q,G), complex-valued to allow for a possible cutoff.
1592 [ + + ]: 6273 : do io=1,nomega
1593 [ + + ]: 6273 : if (omega_distrb(io) == my_rank) then
1594 : : !write(std_out,*)"dim_wing",dim_wing
1595 : : call rpa_symepsm1(iq_ibz,Vcp,npwe,nI,nJ,chi0(:,:,io),my_nqlwl,dim_wing, &
1596 : : chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:), &
1597 : 5007 : tmp_lf,tmp_nlf,tmp_eelf,xmpi_comm_self)
1598 : :
1599 : : ! Store results.
1600 [ + + ]: 10014 : epsm_lf(io,:) = tmp_lf
1601 [ + + ]: 10014 : epsm_nlf(io,:) = tmp_nlf
1602 [ + + ]: 10014 : eelf(io,:) = tmp_eelf
1603 : : end if
1604 : : end do ! nomega
1605 : :
1606 : : case (1)
1607 : : ! Vertex correction from Adiabatic TDDFT. chi_{G1,G2} = [\delta -\chi0 (vc+kxc)]^{-1}_{G1,G3} \chi0_{G3,G2}
1608 [ - + ]: 12 : ABI_CHECK(Vcp%nqlwl==1,"nqlwl/=1 not coded")
1609 [ - + ]: 12 : ABI_CHECK(nkxc==1,"nkxc/=1 not coded")
1610 : :
1611 : : ! Make kxcg_mat(G1,G2) = kxcg(G1-G2) from kxcg defined on the FFT mesh.
1612 [ - + + - : 48 : ABI_MALLOC_OR_DIE(kxcg_mat,(npwe,npwe), ierr)
- + + - +
- + - + -
- + - + ]
1613 : :
1614 : 12 : ierr=0
1615 [ + + ]: 264 : do ig2=1,npwe
1616 [ + + ]: 5988 : do ig1=1,npwe
1617 [ + + ]: 22896 : g1mg2_idx = g2ifft(gvec(:,ig1)-gvec(:,ig2),ngfft)
1618 [ + - ]: 5976 : if (g1mg2_idx>0) then
1619 : 5724 : kxcg_mat(ig1,ig2) = kxcg(g1mg2_idx,1)
1620 : : else
1621 : 0 : ierr=ierr+1
1622 : 0 : kxcg_mat(ig1,ig2) = czero
1623 : : end if
1624 : : end do
1625 : : end do
1626 : :
1627 [ - + ]: 12 : if (ierr/=0) then
1628 : : write(msg,'(a,i4,3a)')&
1629 : 0 : 'Found ',ierr,' G1-G2 vectors falling outside the FFT box. ',ch10,&
1630 : 0 : 'Enlarge the FFT mesh to get rid of this problem. '
1631 : 0 : ABI_WARNING(msg)
1632 : : end if
1633 : :
1634 : : !FIXME "recheck TDDFT code and parallel"
1635 [ - + ]: 12 : ABI_CHECK(nkxc==1,"nkxc/=1 not coded")
1636 [ + + ]: 264 : do io=1,nomega
1637 [ + - ]: 264 : if (omega_distrb(io) == my_rank) then
1638 : : call atddft_symepsm1(iq_ibz,Vcp,npwe,nI,nJ,chi0(:,:,io),kxcg_mat,option_test,my_nqlwl,dim_wing,omega(io),&
1639 [ - + - - : 21404 : chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf, xmpi_comm_self)
- - - - +
- + + + +
+ + + - +
+ + + + +
- + - - -
- + - + +
+ + + - +
+ + + ]
1640 : :
1641 : : ! Store results.
1642 [ + + ]: 504 : epsm_lf(io,:) = tmp_lf
1643 [ + + ]: 504 : epsm_nlf(io,:) = tmp_nlf
1644 [ + + ]: 504 : eelf(io,:) = tmp_eelf
1645 : : end if
1646 : : end do
1647 : :
1648 [ - + ]: 12 : ABI_FREE(kxcg_mat)
1649 [ + + ]: 264 : do io=1,nomega
1650 : 252 : write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]'
1651 : 252 : call wrtout(std_out, msg)
1652 [ + + ]: 516 : call print_arr([std_out], chi0(:,:,io))
1653 : : end do
1654 : :
1655 : : case (2)
1656 : : ! ADA nonlocal vertex correction contained in fxc_ADA
1657 : 0 : ABI_WARNING('Entered fxc_ADA branch: EXPERIMENTAL!')
1658 : : ! Test that argument was passed
1659 [ # # ]: 0 : if (.not.present(fxc_ADA)) then
1660 : 0 : ABI_ERROR('make_epsm1_driver was not called with optional argument fxc_ADA')
1661 : : end if
1662 [ # # ]: 0 : ABI_CHECK(Vcp%nqlwl==1,"nqlwl/=1 not coded")
1663 : :
1664 [ # # ]: 0 : do io=1,nomega
1665 [ # # ]: 0 : if (omega_distrb(io) == my_rank) then
1666 : : call atddft_symepsm1(iq_ibz,Vcp,npwe,nI,nJ,chi0(:,:,io),fxc_ADA,option_test,my_nqlwl,dim_wing,omega(io),&
1667 [ # # # # : 0 : chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,xmpi_comm_self)
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # ]
1668 : :
1669 : : ! Store results.
1670 [ # # ]: 0 : epsm_lf(io,:) = tmp_lf
1671 [ # # ]: 0 : epsm_nlf(io,:) = tmp_nlf
1672 [ # # ]: 0 : eelf(io,:) = tmp_eelf
1673 : : end if
1674 : : end do
1675 : :
1676 [ # # ]: 0 : do io=1,nomega
1677 : 0 : write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]'
1678 : 0 : call wrtout(std_out,msg)
1679 [ # # ]: 0 : call print_arr([std_out], chi0(:,:,io))
1680 : : end do
1681 : :
1682 : : case (4)
1683 : : ! Bootstrap vertex kernel by Sharma [[cite:Sharma2011]]
1684 [ - + + - : 32 : ABI_MALLOC_OR_DIE(vfxc_boot,(npwe*nI,npwe*nJ), ierr)
- + + - +
- + - + -
- + - + ]
1685 [ - + + - : 24 : ABI_MALLOC_OR_DIE(vfxc_boot0,(npwe*nI,npwe*nJ), ierr)
- + + - +
- - + -
+ ]
1686 [ - + + - : 24 : ABI_MALLOC_OR_DIE(chi0_tmp,(npwe*nI,npwe*nJ), ierr)
- + + - +
- - + -
+ ]
1687 [ - + + - : 40 : ABI_MALLOC_OR_DIE(chi0_save,(npwe*nI,npwe*nJ,nomega), ierr)
- + + - -
+ + - + -
+ - + - -
+ - + ]
1688 : :
1689 [ + + ]: 8 : if (iq_ibz==1) then
1690 : 2 : vc_sqrt => Vcp%vcqlwl_sqrt(:,1) ! Use Coulomb term for q-->0
1691 : : else
1692 : 6 : vc_sqrt => Vcp%vc_sqrt(:,iq_ibz)
1693 : : end if
1694 : :
1695 [ - + + + : 19296 : chi0_save = chi0 ! a copy of chi0 (ks)
+ + + + ]
1696 : 8 : nstep = 50 ! max iteration steps
1697 : 8 : alpha = 0.6 ! mixing
1698 : 8 : chi00_head = chi0(1,1,1)*vc_sqrt(1)**2
1699 [ + + + + : 3848 : fxc_head = czero; vfxc_boot = czero; chi0_tmp = czero
+ + + + ]
1700 [ + + + + : 272 : epsm_lf = czero; epsm_nlf = czero; eelf = zero
+ + + + +
+ + + ]
1701 : 8 : write(msg,'(a,2f10.6)') ' -> chi0_dft(head): ', chi00_head
1702 : 8 : call wrtout(std_out,msg)
1703 : : ! loop
1704 : 8 : conv_err = 0.1
1705 [ + - ]: 71 : do istep=1, nstep
1706 [ + + + + : 171181 : chi0 = chi0_save
+ + ]
1707 : 71 : io=1 ! for now only at omega=0
1708 : : call atddft_symepsm1(iq_ibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_boot,0,my_nqlwl,dim_wing,omega(io),&
1709 [ - + - - : 6151 : chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,xmpi_comm_self)
- - - - +
- + - + +
+ + + - +
- + + + +
- + - - -
- + - + +
+ + + - +
+ + + ]
1710 : 71 : conv_err = smallest_real
1711 [ + + ]: 1136 : do ig2=1,npwe*nJ
1712 [ + + ]: 17111 : do ig1=1,npwe*nI
1713 : 17040 : conv_err= MAX(conv_err, ABS(chi0(ig1,ig2,1) - chi0_tmp(ig1,ig2)))
1714 : : end do
1715 : : end do
1716 : 71 : write(msg,'(a,i4,a,f10.6)') ' => bootstrap itr# ', istep, ', eps^-1 max error: ', conv_err
1717 : 71 : call wrtout(std_out,msg)
1718 : 71 : write(msg,'(a,2f10.6)') ' eps^-1(head): ', chi0(1,1,1)
1719 : 71 : call wrtout(std_out,msg)
1720 : 71 : write(msg,'(a,2f10.6)') ' v^-1*fxc(head): ', fxc_head
1721 : 71 : call wrtout(std_out,msg)
1722 [ + + ]: 71 : if (conv_err <= tol4) exit
1723 : : !
1724 [ - + + + : 15246 : chi0_tmp = chi0(:,:,1)
+ + ]
1725 [ - + + + : 15246 : vfxc_boot = chi0(:,:,1)/chi00_head
+ + ]
1726 [ + + ]: 63 : if (istep > 1) then
1727 [ - + + + : 13310 : vfxc_boot = alpha*vfxc_boot0 + (one-alpha)*vfxc_boot
+ + ]
1728 : : end if
1729 [ - + + + : 15246 : vfxc_boot0 = vfxc_boot
+ + ]
1730 : 63 : fxc_head = vfxc_boot(1,1)
1731 [ + + ]: 1016 : do ig1=1,npwe
1732 [ + + ]: 15183 : vfxc_boot(ig1,:) = vc_sqrt(ig1)*vc_sqrt(:)*vfxc_boot(ig1,:)
1733 : : end do
1734 : : end do
1735 : : ! end loop
1736 [ + - ]: 8 : if (istep <= nstep) then
1737 : 8 : write(msg,'(a,i4,a)') ' => bootstrap fxc converged after ', istep, ' iterations'
1738 : 8 : call wrtout(std_out,msg)
1739 : : else
1740 : 0 : write(msg,'(a,i4,a)') ' -> bootstrap fxc not converged after ', nstep, ' iterations'
1741 : 0 : ABI_WARNING(msg)
1742 : : end if
1743 : : !
1744 [ + + + + : 19288 : chi0 = chi0_save
+ + ]
1745 [ + + ]: 88 : do io=1,nomega
1746 [ + - ]: 88 : if (omega_distrb(io) == my_rank) then
1747 : : call atddft_symepsm1(iq_ibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_boot,option_test,my_nqlwl,dim_wing,omega(io),&
1748 [ - + - - : 7760 : chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,xmpi_comm_self)
- - - - +
- + - + +
+ + + - +
- + + + +
- + - - -
- + - + +
+ + + - +
+ + + ]
1749 [ + + ]: 160 : epsm_lf(io,:) = tmp_lf
1750 [ + + ]: 160 : epsm_nlf(io,:) = tmp_nlf
1751 [ + + ]: 160 : eelf(io,:) = tmp_eelf
1752 : : end if
1753 : : end do
1754 : :
1755 [ - + ]: 8 : ABI_FREE(chi0_tmp)
1756 [ - + ]: 8 : ABI_FREE(chi0_save)
1757 [ - + ]: 8 : ABI_FREE(vfxc_boot)
1758 [ - + ]: 8 : ABI_FREE(vfxc_boot0)
1759 : :
1760 [ + + ]: 88 : do io=1,nomega
1761 : 80 : write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]'
1762 : 80 : call wrtout(std_out,msg)
1763 [ + + ]: 168 : call print_arr([std_out], chi0(:,:,io))
1764 : : end do
1765 : :
1766 : : case (5)
1767 : : ! One-shot scalar bootstrap approximation
1768 [ - + + - : 16 : ABI_MALLOC_OR_DIE(vfxc_boot,(npwe*nI,npwe*nJ), ierr)
- + + - +
- + - + -
- + - + ]
1769 [ - + + - : 20 : ABI_MALLOC_OR_DIE(chi0_save,(npwe*nI,npwe*nJ,nomega), ierr)
- + + - -
+ + - + -
+ - + - -
+ - + ]
1770 : :
1771 [ + + ]: 4 : if (iq_ibz==1) then
1772 : 1 : vc_sqrt => Vcp%vcqlwl_sqrt(:,1) ! Use Coulomb term for q-->0
1773 : : else
1774 : 3 : vc_sqrt => Vcp%vc_sqrt(:,iq_ibz)
1775 : : end if
1776 : :
1777 [ - + + + : 9648 : chi0_save = chi0 ! a copy of chi0
+ + + + ]
1778 [ + + + + ]: 964 : fxc_head = czero; vfxc_boot = czero;
1779 [ + + + + : 136 : epsm_lf = czero; epsm_nlf = czero; eelf = zero
+ + + + +
+ + + ]
1780 : 4 : chi00_head = chi0(1,1,1)*vc_sqrt(1)**2
1781 : 4 : write(msg,'(a,2f10.6)') ' -> chi0_dft(head): ',chi00_head
1782 : 4 : call wrtout(std_out,msg)
1783 : :
1784 : 4 : fxc_head = vc_sqrt(1)**2/chi00_head + vc_sqrt(1)**2/chi00_head - vc_sqrt(1)**2
1785 : 4 : fxc_head = 0.5*fxc_head + 0.5*sqrt(fxc_head**2 - 4.0*vc_sqrt(1)**4/(chi00_head*chi00_head))
1786 : 4 : vfxc_boot(1,1) = fxc_head
1787 : 4 : write(msg,'(a,2f10.6)') ' -> v^-1*fxc(head): ',fxc_head/vc_sqrt(1)**2
1788 : 4 : call wrtout(std_out,msg)
1789 : :
1790 [ + + + + : 9644 : chi0 = chi0_save
+ + ]
1791 [ + + ]: 44 : do io=1,nomega
1792 [ + - ]: 44 : if (omega_distrb(io) == my_rank) then
1793 : : call atddft_symepsm1(iq_ibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_boot,option_test,my_nqlwl,dim_wing,omega(io),&
1794 [ - + - - : 3880 : chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,xmpi_comm_self)
- - - - +
- + - + +
+ + + - +
- + + + +
- + - - -
- + - + +
+ + + - +
+ + + ]
1795 [ + + ]: 80 : epsm_lf(io,:) = tmp_lf
1796 [ + + ]: 80 : epsm_nlf(io,:) = tmp_nlf
1797 [ + + ]: 80 : eelf(io,:) = tmp_eelf
1798 : : end if
1799 : : end do
1800 : 4 : write(msg,'(a,2f10.6)') ' eps^-1(head): ',chi0(1,1,1)
1801 : 4 : call wrtout(std_out,msg)
1802 : :
1803 [ - + ]: 4 : ABI_FREE(chi0_save)
1804 : 4 : ABI_FREE(vfxc_boot)
1805 : :
1806 [ + + ]: 44 : do io=1,nomega
1807 : 40 : write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]'
1808 : 40 : call wrtout(std_out,msg)
1809 [ + + ]: 84 : call print_arr([std_out], chi0(:,:,io))
1810 : : end do
1811 : :
1812 : : case (6)
1813 : : ! RPA bootstrap by Rigamonti [[cite:Rigamonti2015]] and Berger [[cite:Berger2015]]
1814 [ # # # # : 0 : ABI_MALLOC_OR_DIE(vfxc_boot,(npwe*nI,npwe*nJ), ierr)
# # # # #
# # # # #
# # # # ]
1815 [ # # # # : 0 : ABI_MALLOC_OR_DIE(chi0_save,(npwe*nI,npwe*nJ,nomega), ierr)
# # # # #
# # # # #
# # # # #
# # # ]
1816 [ # # # # : 0 : ABI_MALLOC_OR_DIE(chi0_tmp,(npwe*nI,npwe*nJ), ierr)
# # # # #
# # # #
# ]
1817 : :
1818 [ # # ]: 0 : if (iq_ibz==1) then
1819 : 0 : vc_sqrt => Vcp%vcqlwl_sqrt(:,1) ! Use Coulomb term for q-->0
1820 : : else
1821 : 0 : vc_sqrt => Vcp%vc_sqrt(:,iq_ibz)
1822 : : end if
1823 : :
1824 [ # # # # : 0 : chi0_save = chi0 ! a copy of chi0
# # # # ]
1825 [ # # # # ]: 0 : fxc_head = czero; vfxc_boot = czero;
1826 [ # # # # : 0 : epsm_lf = czero; epsm_nlf = czero; eelf = zero
# # # # #
# # # ]
1827 : : !chi00_head = chi0(1,1,1)*vc_sqrt(1)**2
1828 : 0 : write(msg,'(a,2f10.6)') ' -> chi0_dft(head): ',chi00_head
1829 : 0 : call wrtout(std_out,msg)
1830 : :
1831 : 0 : io = 1 ! static
1832 : : call atddft_symepsm1(iq_ibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_boot,0,my_nqlwl,dim_wing,omega(io),&
1833 [ # # # # : 0 : chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,xmpi_comm_self)
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # ]
1834 [ # # ]: 0 : epsm_lf(1,:) = tmp_lf
1835 : :
1836 : : ! chi(RPA) = chi0 * (1 - chi0 * v_c)^-1
1837 [ # # # # : 0 : chi0 = chi0_save
# # ]
1838 [ # # ]: 0 : do ig2=2,npwe
1839 [ # # ]: 0 : do ig1=2,npwe
1840 : 0 : chi0(ig1,ig2,io)=-vc_sqrt(ig1)*chi0(ig1,ig2,io)*vc_sqrt(ig2)
1841 : : end do
1842 : 0 : chi0(ig2,ig2,io)=one+chi0(ig2,ig2,io)
1843 : : end do
1844 [ # # # # ]: 0 : chi0(1,:,io) = czero; chi0(:,1,io) = czero; chi0(1,1,io) = one
1845 [ # # # # : 0 : chi0_tmp = chi0(:,:,io)
# # ]
1846 : 0 : call xginv(chi0_tmp,npwe,comm=comm)
1847 [ # # # # : 0 : chi0 = chi0_save
# # ]
1848 [ # # # # : 0 : chi0_tmp = MATMUL(chi0(:,:,io), chi0_tmp(:,:)) ! chi(RPA)
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # ]
1849 [ # # ]: 0 : do ig1=1,npwe
1850 [ # # ]: 0 : chi0_tmp(ig1,:) = vc_sqrt(ig1)*vc_sqrt(:)*chi0_tmp(ig1,:)
1851 : : end do
1852 : : !call xginv(chi0_tmp,npwe,comm=comm) ! chi(RPA)^-1
1853 : : !vfxc_boot = chi0_tmp/epsm_lf(1,1)
1854 : : !
1855 : : !vfxc_boot(1,1) = chi0_tmp(1,1)/epsm_lf(1,1)
1856 : 0 : vfxc_boot(1,1) = one/chi0_tmp(1,1)/epsm_lf(1,1)
1857 : : !@WC: alternatively:
1858 : : !chi00_head = chi0(1,1,io)*vc_sqrt(1)**2
1859 : : !vfxc_boot(1,1) = one/chi00_head/epsm_lf(1,1)
1860 : 0 : fxc_head = vfxc_boot(1,1)
1861 [ # # ]: 0 : do ig1=1,npwe
1862 [ # # ]: 0 : vfxc_boot(ig1,:) = vc_sqrt(ig1)*vc_sqrt(:)*vfxc_boot(ig1,:)
1863 : : end do
1864 : 0 : write(msg,'(a,2f10.6)') ' -> v^-1*fxc(head): ',fxc_head
1865 : 0 : call wrtout(std_out,msg)
1866 : :
1867 [ # # # # : 0 : chi0 = chi0_save
# # ]
1868 [ # # ]: 0 : do io=1,nomega
1869 [ # # ]: 0 : if (omega_distrb(io) == my_rank) then
1870 : : call atddft_symepsm1(iq_ibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_boot,option_test,my_nqlwl,dim_wing,omega(io),&
1871 [ # # # # : 0 : chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,xmpi_comm_self)
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # ]
1872 [ # # ]: 0 : epsm_lf(io,:) = tmp_lf
1873 [ # # ]: 0 : epsm_nlf(io,:) = tmp_nlf
1874 [ # # ]: 0 : eelf(io,:) = tmp_eelf
1875 : : end if
1876 : : end do
1877 : 0 : write(msg,'(a,2f10.6)') ' eps^-1(head): ',chi0(1,1,1)
1878 : 0 : call wrtout(std_out,msg)
1879 : :
1880 [ # # ]: 0 : ABI_FREE(chi0_save)
1881 : 0 : ABI_FREE(vfxc_boot)
1882 : 0 : ABI_FREE(chi0_tmp)
1883 : :
1884 [ # # ]: 0 : do io=1,nomega
1885 : 0 : write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]'
1886 : 0 : call wrtout(std_out,msg)
1887 [ # # ]: 0 : call print_arr([std_out], chi0(:,:,io))
1888 : : end do
1889 : :
1890 : : case (7)
1891 : : ! LR+ALDA hybrid vertex kernel by Tal
1892 : : ! First ALDA
1893 [ - + ]: 4 : ABI_CHECK(Vcp%nqlwl==1,"nqlwl/=1 not coded")
1894 [ - + ]: 4 : ABI_CHECK(nkxc==1,"nkxc/=1 not coded")
1895 : : ! Make kxcg_mat(G1,G2) = kxcg(G1-G2) from kxcg defined on the FFT mesh.
1896 [ - + + - : 16 : ABI_MALLOC_OR_DIE(kxcg_mat,(npwe,npwe), ierr)
- + + - +
- + - + -
- + - + ]
1897 : 4 : ierr=0
1898 [ + + ]: 64 : do ig2=1,npwe
1899 [ + + ]: 964 : do ig1=1,npwe
1900 [ + + ]: 3600 : g1mg2_idx = g2ifft(gvec(:,ig1)-gvec(:,ig2),ngfft)
1901 [ + - ]: 960 : if (g1mg2_idx>0) then
1902 : 900 : kxcg_mat(ig1,ig2) = kxcg(g1mg2_idx,1)
1903 : : else
1904 : 0 : ierr=ierr+1
1905 : 0 : kxcg_mat(ig1,ig2) = czero
1906 : : end if
1907 : : end do
1908 : : end do
1909 [ - + ]: 4 : if (ierr/=0) then
1910 : : write(msg,'(a,i4,3a)')&
1911 : 0 : 'Found ',ierr,' G1-G2 vectors falling outside the FFT box. ',ch10,&
1912 : 0 : 'Enlarge the FFT mesh to get rid of this problem. '
1913 : 0 : ABI_WARNING(msg)
1914 : : end if
1915 : : !FIXME "recheck TDDFT code and parallel"
1916 [ - + ]: 4 : ABI_CHECK(nkxc==1,"nkxc/=1 not coded")
1917 : :
1918 : : ! Now LR: (1-Z)*chi0^-1
1919 [ - + + - : 16 : ABI_MALLOC_OR_DIE(vfxc_lr,(npwe*nI,npwe*nJ), ierr)
- + + - +
- + - + -
- + - + ]
1920 [ - + + - : 12 : ABI_MALLOC_OR_DIE(vfxc_tmp,(npwe*nI,npwe*nJ), ierr)
- + + - +
- - + -
+ ]
1921 [ - + + - : 12 : ABI_MALLOC_OR_DIE(chi0_tmp,(npwe*nI,npwe*nJ), ierr)
- + + - +
- - + -
+ ]
1922 : :
1923 [ + + ]: 4 : if (iq_ibz==1) then
1924 : 1 : vc_sqrt => Vcp%vcqlwl_sqrt(:,1) ! Use Coulomb term for q-->0
1925 : : else
1926 : 3 : vc_sqrt => Vcp%vc_sqrt(:,iq_ibz)
1927 : : end if
1928 : :
1929 : 4 : Zr = 0.78
1930 : 4 : chi00_head = chi0(1,1,1)*vc_sqrt(1)**2
1931 [ + + + + : 1924 : fxc_head = czero; vfxc_lr = czero; vfxc_tmp = czero
+ + + + ]
1932 [ + + + + : 136 : epsm_lf = czero; epsm_nlf = czero; eelf = zero
+ + + + +
+ + + ]
1933 : 4 : write(msg,'(a,2f10.6)') ' -> chi0_dft(head): ', chi00_head
1934 : 4 : call wrtout(std_out,msg)
1935 : : !
1936 [ - + + + : 968 : chi0_tmp = chi0(:,:,1)
+ + ]
1937 : 4 : call xginv(chi0_tmp,npwe,comm=comm)
1938 [ - + + + : 968 : vfxc_lr = (one-Zr)*chi0_tmp(:,:)
+ + ]
1939 : 4 : write(msg,'(a)') ' Constructing LR+ALDA fxc kernel'
1940 : 4 : call wrtout(std_out,msg)
1941 : : !
1942 [ + + ]: 64 : do ig1=1,npwe
1943 [ + + ]: 964 : do ig2=1,npwe
1944 [ + + ]: 3600 : qpg2 = Vcp%qibz(:,iq_ibz) + gvec(:,ig1)
1945 : 900 : qpg2_nrm = normv(qpg2,gmet,"G")
1946 [ + + ]: 3600 : qpg2 = Vcp%qibz(:,iq_ibz) + gvec(:,ig2)
1947 : 900 : qpg2_nrm = SQRT(qpg2_nrm * normv(qpg2,gmet,"G"))
1948 : : vfxc_tmp(ig1,ig2) = vfxc_lr(ig1,ig2)*exp(-(qpg2_nrm/k_thfermi(rhor))**2) + &
1949 : 960 : kxcg_mat(ig1,ig2)*(one - exp(-(qpg2_nrm/k_thfermi(rhor))**2))
1950 : : !write(std_out,*) ig1, qpg2_nrm, k_thfermi(rhor), vfxc_lr(ig1,ig1), kxcg_mat(ig1,ig1), vfxc_tmp(ig1,ig1)
1951 : : end do
1952 : : end do
1953 : : !
1954 [ - + + + : 968 : vfxc_lr = vfxc_tmp
+ + ]
1955 : :
1956 [ + + ]: 44 : do io=1,nomega
1957 [ + - ]: 44 : if (omega_distrb(io) == my_rank) then
1958 : : call atddft_hyb_symepsm1(iq_ibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_lr,kxcg_mat,option_test,my_nqlwl,dim_wing,omega(io),&
1959 [ - + - - : 3880 : chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,xmpi_comm_self)
- - - - +
- + - + +
+ + + - +
- + + + +
- + - - -
- + - + +
+ + + - +
+ + + ]
1960 [ + + ]: 80 : epsm_lf(io,:) = tmp_lf
1961 [ + + ]: 80 : epsm_nlf(io,:) = tmp_nlf
1962 [ + + ]: 80 : eelf(io,:) = tmp_eelf
1963 : : end if
1964 : : end do
1965 : :
1966 [ - + ]: 4 : ABI_FREE(kxcg_mat)
1967 [ - + ]: 4 : ABI_FREE(chi0_tmp)
1968 [ - + ]: 4 : ABI_FREE(vfxc_lr)
1969 [ - + ]: 4 : ABI_FREE(vfxc_tmp)
1970 : :
1971 [ + + ]: 44 : do io=1,nomega
1972 : 40 : write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]'
1973 : 40 : call wrtout(std_out,msg)
1974 [ + + ]: 84 : call print_arr([std_out], chi0(:,:,io))
1975 : : end do
1976 : :
1977 : : case default
1978 [ + + - + : 950 : ABI_BUG(sjoin('Wrong approx_type:',itoa(approx_type)))
+ - + - ]
1979 : : end select
1980 : :
1981 [ - + ]: 950 : if (use_mpi_shared_win) call xmpi_win_fence(XMPI_MODE_NOSUCCEED, epsm1_win, ierr) ! Close the RMA epoch.
1982 : :
1983 [ + + ]: 950 : if (use_MPI .and. .not. use_mpi_shared_win) then
1984 : : ! Collect results on each node.
1985 [ - + + - : 1188 : ABI_MALLOC(buffer_lwing, (size(chi0_lwing,dim=1), size(chi0_lwing, dim=3)))
- + + - +
- + - - +
- + - + ]
1986 [ - + + - : 1188 : ABI_MALLOC(buffer_uwing, (size(chi0_uwing,dim=1), size(chi0_uwing, dim=3)))
- + + - +
- + - - +
- + - + ]
1987 : :
1988 [ + + ]: 852 : do io=1,nomega
1989 [ + + ]: 555 : if (omega_distrb(io)/=my_rank) then
1990 : : ! Zero arrays.
1991 [ + + + + ]: 967208 : chi0(:,:,io) = czero_gw
1992 [ + + ]: 344 : if(dim_wing>0) then
1993 [ + + + + ]: 5904 : chi0_lwing(:,io,:) = zero
1994 [ + + + + ]: 5904 : chi0_uwing(:,io,:) = zero
1995 [ + + + + ]: 624 : chi0_head(:,:,io) = czero
1996 : : endif
1997 [ + + ]: 688 : epsm_lf(io,:) = czero
1998 [ + + ]: 688 : epsm_nlf(io,:) = czero
1999 [ + + ]: 688 : eelf(io,:) = zero
2000 : : end if
2001 : :
2002 : 555 : call xmpi_sum(chi0(:,:,io), comm,ierr)
2003 : :
2004 [ + + ]: 852 : if (dim_wing>0) then
2005 : : ! Build contiguous arrays
2006 [ + - + - : 10562 : buffer_lwing = chi0_lwing(:,io,:)
- + - - +
+ + + ]
2007 [ + - + - : 10562 : buffer_uwing = chi0_uwing(:,io,:)
- + - - +
+ + + ]
2008 : 79 : call xmpi_sum(buffer_lwing,comm,ierr)
2009 : 79 : call xmpi_sum(buffer_uwing,comm,ierr)
2010 [ + + + + ]: 10483 : chi0_lwing(:,io,:) = buffer_lwing
2011 [ + + + + ]: 10483 : chi0_uwing(:,io,:) = buffer_uwing
2012 [ + + + - ]: 237 : if (size(chi0_head(:,:,io))/= zero) then
2013 : 79 : call xmpi_sum(chi0_head(:,:,io),comm,ierr)
2014 : : end if
2015 : : end if
2016 : :
2017 : : end do ! iw
2018 : :
2019 : 297 : call xmpi_sum(epsm_lf, comm, ierr)
2020 : 297 : call xmpi_sum(epsm_nlf,comm, ierr)
2021 : 297 : call xmpi_sum(eelf, comm, ierr)
2022 [ - + ]: 297 : ABI_FREE(buffer_lwing)
2023 [ - + ]: 594 : ABI_FREE(buffer_uwing)
2024 : : end if
2025 : :
2026 : : ! Save results in Spectra%, mind the slicing.
2027 [ + + + + ]: 4633 : Spectra%emacro_nlf(:,:) = epsm_nlf(1:nor,:)
2028 [ + + + + ]: 4633 : Spectra%emacro_lf (:,:) = epsm_lf (1:nor,:)
2029 [ + + + + ]: 4633 : Spectra%eelf (:,:) = eelf (1:nor,:)
2030 : :
2031 [ - + ]: 950 : ABI_FREE(epsm_lf)
2032 [ - + ]: 950 : ABI_FREE(epsm_nlf)
2033 [ - + ]: 950 : ABI_FREE(eelf)
2034 [ - + ]: 950 : ABI_FREE(tmp_lf)
2035 [ - + ]: 950 : ABI_FREE(tmp_nlf)
2036 [ - + ]: 950 : ABI_FREE(tmp_eelf)
2037 : :
2038 : 950 : call cwtime_report(" make_epsm1_driver", cpu, wall, gflops)
2039 : :
2040 [ - + - + : 2850 : end subroutine make_epsm1_driver
- + ]
2041 : : !!***
2042 : :
2043 : : !----------------------------------------------------------------------
2044 : :
2045 : : !!****f* m_screening/rpa_symepsm1
2046 : : !! NAME
2047 : : !! rpa_symepsm1
2048 : : !!
2049 : : !! FUNCTION
2050 : : !! Calculate RPA $\tilde\epsilon^{-1}$
2051 : : !!
2052 : : !! Use a special treatment of non-Analytic behavior of heads and wings in reciprocal space
2053 : : !! calculating these quantities for different small q-directions specified by the user
2054 : : !! (Not yet operative)
2055 : : !!
2056 : : !! INPUTS
2057 : : !! iq_ibz=index of the q-point in the array Vcp%qibz where epsilon^-1 has to be calculated
2058 : : !! Vcp<vcoul_t>=Structure gathering data on the Coulombian interaction
2059 : : !! npwe=Number of G-vectors in chi0.
2060 : : !! nI,nJ=Number of rows/columns in chi0_ij (1,1 in collinear case)
2061 : : !! dim_wing=Dimension of the wings (0 or 3 if q-->0)
2062 : : !! chi0_head(dim_wing,dim_wing)=Head of of chi0 (only for q-->0)
2063 : : !! chi0_lwing(npwe*nI,dim_wing)=Lower wings of chi0 (only for q-->0)
2064 : : !! chi0_uwing(npwe*nJ,dim_wing)=Upper wings of chi0 (only for q-->0)
2065 : : !! comm=MPI communicator.
2066 : : !!
2067 : : !! SIDE EFFECTS
2068 : : !! chi0(npwe*nI,npwe*nJ): in input the irreducible polarizability,
2069 : : !! in output the symmetrized inverse dielectric matrix.
2070 : : !!
2071 : : !! SOURCE
2072 : :
2073 [ - + - + ]: 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)
2074 : :
2075 : : !Arguments ------------------------------------
2076 : : !scalars
2077 : : integer,intent(in) :: iq_ibz,nI,nJ,npwe,dim_wing,my_nqlwl,comm
2078 : : type(vcoul_t),target,intent(in) :: Vcp
2079 : : !arrays
2080 : : complex(gwp),intent(inout) :: chi0(npwe*nI,npwe*nJ)
2081 : : complex(dp),intent(inout) :: chi0_lwing(:,:) !(npwe*nI,dim_wing)
2082 : : complex(dp),intent(inout) :: chi0_uwing(:,:) !(npwe*nJ,dim_wing)
2083 : : complex(dp),intent(inout) :: chi0_head(:,:) !(dim_wing,dim_wing)
2084 : : real(dp),intent(out) :: eelf(my_nqlwl)
2085 : : complex(dp),intent(out) :: epsm_lf(my_nqlwl),epsm_nlf(my_nqlwl)
2086 : :
2087 : : !Local variables-------------------------------
2088 : : !scalars
2089 : : integer,parameter :: master=0,prtvol=0
2090 : : integer :: ig1,ig2,iqlwl,my_rank,nprocs
2091 : : real(dp) :: ucvol
2092 : : logical :: is_qeq0
2093 : : !character(len=500) :: msg
2094 : : !arrays
2095 : : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
2096 : 5007 : complex(gwp), contiguous, pointer :: vc_sqrt(:)
2097 : 5007 : complex(gwp),allocatable :: chi0_save(:,:)
2098 : : ! *************************************************************************
2099 : :
2100 : : ABI_UNUSED(chi0_head(1,1))
2101 : :
2102 [ - + ]: 5007 : if (nI/=1.or.nJ/=1) then
2103 : 0 : ABI_ERROR("nI or nJ=/1 not yet implemented")
2104 : : end if
2105 : :
2106 : 5007 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
2107 : :
2108 : 5007 : call metric(gmet,gprimd,-1,rmet,Vcp%rprimd,ucvol)
2109 : :
2110 : 5007 : is_qeq0 = (normv(Vcp%qibz(:,iq_ibz),gmet,'G')<GW_TOLQ0)
2111 [ + + ]: 5007 : if (is_qeq0) then
2112 [ - + ]: 790 : ABI_CHECK(iq_ibz==1, "q is 0 but iq_ibz /= 1")
2113 : : end if
2114 : :
2115 [ - + ]: 5007 : if (my_nqlwl>1) then
2116 [ # # # # : 0 : ABI_MALLOC(chi0_save,(npwe*nI,npwe*nJ))
# # # # #
# # # # #
# # ]
2117 [ # # # # : 0 : chi0_save = chi0
# # ]
2118 : : end if
2119 : : !
2120 : : ! Symmetrized RPA epsilon: 1 - Vc^{1/2} chi0 Vc^{1/2}
2121 : : ! vc_sqrt contains vc^{1/2}(q, G)
2122 : : !
2123 : : ! Loop over small q"s (if any) to treat the nonanalytical behavior.
2124 [ + + ]: 10014 : do iqlwl=my_nqlwl,1,-1
2125 : :
2126 [ - + ]: 5007 : if (my_nqlwl>1) then
2127 [ # # # # ]: 0 : chi0(:,:) = chi0_save ! restore pristine polarizability
2128 [ # # ]: 0 : chi0(:,1) = chi0_lwing(:,iqlwl) ! change the wings
2129 [ # # ]: 0 : chi0(1,:) = chi0_uwing(:,iqlwl)
2130 : : end if
2131 : :
2132 [ + + ]: 5007 : if (iq_ibz==1) then
2133 : 790 : vc_sqrt => Vcp%vcqlwl_sqrt(:,iqlwl) ! Use Coulomb term for q-->0
2134 : : else
2135 : 4217 : vc_sqrt => Vcp%vc_sqrt(:,iq_ibz)
2136 : : end if
2137 : :
2138 [ + + ]: 227716 : do ig2=1,npwe*nJ
2139 [ + + ]: 17303524 : do ig1=1,npwe*nI
2140 : 17303524 : chi0(ig1,ig2) = -vc_sqrt(ig1) * chi0(ig1,ig2) * vc_sqrt(ig2)
2141 : : end do
2142 : 227716 : chi0(ig2,ig2) = one + chi0(ig2,ig2)
2143 : : end do
2144 : : ! chi0, now contains \tepsilon.
2145 : :
2146 : 5007 : epsm_nlf(iqlwl)=chi0(1,1)
2147 : :
2148 : : if (prtvol > 0) then
2149 : : call wrtout(std_out,' Symmetrical epsilon(G,G'') ')
2150 : : call print_arr([std_out], chi0)
2151 : : end if
2152 : : !
2153 : : ! === Invert tepsilon and calculate macroscopic dielectric constant ===
2154 : : ! * epsm_lf(w) = 1 / epsm1(G=0,Gp=0,w).
2155 : : ! * Since G=Gp=0, there is no difference btw symmetrical and not symmetrical.
2156 : :
2157 : 5007 : call xginv(chi0,npwe,comm=comm)
2158 : :
2159 : 5007 : epsm_lf(iqlwl) = one/chi0(1,1)
2160 : 5007 : eelf(iqlwl) = -AIMAG(chi0(1,1))
2161 : :
2162 : : if (prtvol > 0) then
2163 : : call wrtout(std_out," Symmetrical epsilon^-1(G,G'')")
2164 : : call print_arr([std_out], chi0)
2165 : : end if
2166 : : !
2167 : : ! Save wings of e^-1 overwriting input values.
2168 : 5007 : if (dim_wing>0.and..FALSE.) then
2169 : : chi0_lwing(:,iqlwl) = chi0(:,1)
2170 : : chi0_uwing(:,iqlwl) = chi0(1,:)
2171 : : end if
2172 : :
2173 : : end do !iqlwl
2174 : :
2175 [ - + ]: 5007 : ABI_SFREE(chi0_save)
2176 : :
2177 [ + - ]: 5007 : end subroutine rpa_symepsm1
2178 : : !!***
2179 : :
2180 : : !----------------------------------------------------------------------
2181 : :
2182 : : !!****f* m_screening/atddft_symepsm1
2183 : : !! NAME
2184 : : !! atddft_symepsm1
2185 : : !!
2186 : : !! FUNCTION
2187 : : !! Calculate $\tilde\epsilon^{-1}$ using ALDA within TDDFT
2188 : : !!
2189 : : !! 2) Use a special treatment of non-Analytic behavior of heads and wings in reciprocal space
2190 : : !! calculating these quantities for different small q-directions specified by the user
2191 : : !! (Not yet operative)
2192 : : !!
2193 : : !! Output the electron energy loss function and the macroscopic dielectric function with and
2194 : : !! without local field effects (only if non-zero real frequencies are available)
2195 : : !!
2196 : : !! INPUTS
2197 : : !! iq_ibz=index of the q-point in the array Vcp%qibz where epsilon^-1 has to be calculated
2198 : : !! npwe=Number of G-vectors in chi0.
2199 : : !! nI,nJ=Number of rows/columns in chi0_ij (1,1 in collinear case)
2200 : : !! dim_wing=Dimension of the wings (0 or 3 if q-->0)
2201 : : !! option_test=Only for TDDFT:
2202 : : !! == 0 for TESTPARTICLE ==
2203 : : !! == 1 for TESTELECTRON ==
2204 : : !! Vcp<vcoul_t>=Structure gathering data on the Coulombian interaction
2205 : : !! %nqibz=Number of q-points.
2206 : : !! %qibz(3,nqibz)=q-points in the IBZ.
2207 : : !! comm=MPI communicator.
2208 : : !! chi0_lwing(npwe*nI,dim_wing)=Lower wings of chi0 (only for q-->0)
2209 : : !! chi0_uwing(npwe*nJ,dim_wing)=Upper wings of chi0 (only for q-->0)
2210 : : !! chi0_head(dim_wing,dim_wing)=Head of of chi0 (only for q-->0)
2211 : : !!
2212 : : !! OUTPUT
2213 : : !!
2214 : : !! SIDE EFFECTS
2215 : : !! chi0(npwe*nI,npwe*nJ): in input the irreducible polarizability, in output
2216 : : !! the symmetrized inverse dielectric matrix.
2217 : : !!
2218 : : !! SOURCE
2219 : :
2220 : 443 : subroutine atddft_symepsm1(iq_ibz,Vcp,npwe,nI,nJ,chi0,kxcg_mat,option_test,my_nqlwl,dim_wing,omega,&
2221 : 443 : & chi0_head,chi0_lwing,chi0_uwing,epsm_lf,epsm_nlf,eelf,comm)
2222 : :
2223 : : !Arguments ------------------------------------
2224 : : !scalars
2225 : : integer,intent(in) :: iq_ibz,nI,nJ,npwe,dim_wing,my_nqlwl
2226 : : integer,intent(in) :: option_test,comm
2227 : : type(vcoul_t),target,intent(in) :: Vcp
2228 : : !arrays
2229 : : complex(gwp),intent(in) :: kxcg_mat(npwe,npwe)
2230 : : complex(dp),intent(in) :: omega
2231 : : complex(dp),intent(inout) :: chi0_lwing(npwe*nI,dim_wing)
2232 : : complex(dp),intent(inout) :: chi0_uwing(npwe*nJ,dim_wing)
2233 : : complex(dp),intent(inout) :: chi0_head(dim_wing,dim_wing)
2234 : : complex(gwp),intent(inout) :: chi0(npwe*nI,npwe*nJ)
2235 : : real(dp),intent(out) :: eelf(my_nqlwl)
2236 : : complex(dp),intent(out) :: epsm_lf(my_nqlwl),epsm_nlf(my_nqlwl)
2237 : :
2238 : : !Local variables-------------------------------
2239 : : !scalars
2240 : : integer,parameter :: master=0,prtvol=0
2241 : : integer :: ig1,ig2,my_rank,nprocs,ierr
2242 : : real(dp) :: ucvol
2243 : : logical :: is_qeq0
2244 : : character(len=500) :: msg
2245 : : !arrays
2246 : : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
2247 : 443 : complex(gwp),allocatable :: chitmp(:,:)
2248 : 443 : complex(gwp), contiguous, pointer :: vc_sqrt(:)
2249 : : ! *************************************************************************
2250 : :
2251 : : ABI_UNUSED(chi0_head(1,1))
2252 : : ABI_UNUSED(chi0_lwing(1,1))
2253 : : ABI_UNUSED(chi0_uwing(1,1))
2254 : :
2255 [ - + ]: 443 : if (nI/=1.or.nJ/=1) then
2256 : 0 : ABI_ERROR("nI or nJ=/1 not yet implemented")
2257 : : end if
2258 : :
2259 [ - + ]: 443 : ABI_CHECK(Vcp%nqlwl==1,"nqlwl/=1 not coded")
2260 [ - + ]: 443 : ABI_CHECK(my_nqlwl==1,"my_nqlwl/=1 not coded")
2261 : :
2262 : 443 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
2263 : :
2264 : 443 : call metric(gmet,gprimd,-1,rmet,Vcp%rprimd,ucvol)
2265 : :
2266 : 443 : is_qeq0 = (normv(Vcp%qibz(:,iq_ibz),gmet,'G')<GW_TOLQ0)
2267 : :
2268 [ + + ]: 443 : if (iq_ibz==1) then
2269 : : !%vc_sqrt => Vcp%vcqlwl_sqrt(:,iqlwl) ! Use Coulomb term for q-->0
2270 : 84 : vc_sqrt => Vcp%vcqlwl_sqrt(:,1) ! TODO add treatment of non-Analytic behavior
2271 : : else
2272 : 359 : vc_sqrt => Vcp%vc_sqrt(:,iq_ibz)
2273 : : end if
2274 : :
2275 : 443 : write(msg,'(a,f8.2,a)')" chitmp requires: ",npwe**2*gwp*b2Mb," Mb"
2276 [ - + + - : 1772 : ABI_MALLOC_OR_DIE(chitmp, (npwe,npwe), ierr)
- + + - +
- + - + -
- + - + ]
2277 : : !
2278 : : ! Calculate chi0*fxc.
2279 [ + - + - : 3607552 : chitmp = MATMUL(chi0,kxcg_mat)
+ - - + -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - + +
+ + + + +
+ + + + -
+ - - - -
- - - ]
2280 : : ! First calculate the NLF contribution
2281 [ + + ]: 7232 : do ig1=1,npwe
2282 [ + + ]: 112512 : do ig2=1,npwe
2283 : 112512 : chitmp(ig1,ig2)=-chitmp(ig1,ig2)
2284 : : end do
2285 : 7232 : chitmp(ig1,ig1)=chitmp(ig1,ig1)+one
2286 : : end do
2287 : :
2288 : 443 : call xginv(chitmp,npwe,comm=comm)
2289 : :
2290 [ + - + - : 7102592 : chitmp = MATMUL(chitmp,chi0)
+ + + + +
- + - + -
- + - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- + + + +
+ + + + +
+ + - + -
+ - + - -
- - - - -
+ - ]
2291 : : !if (.not. ABS(REAL(omega))> tol3) call hermitianize(chitmp,"All")
2292 : 443 : chitmp(1,1)=-vc_sqrt(1)*chitmp(1,1)*vc_sqrt(1)
2293 : 443 : chitmp(1,1)=chitmp(1,1)+one
2294 : :
2295 : 443 : epsm_nlf(1)=chitmp(1,1)
2296 : :
2297 [ + - + - : 3607109 : chitmp = MATMUL(chi0,kxcg_mat)
- + - - -
- - - - -
- - - - -
- - - + +
+ + + + +
+ + + + -
+ - - - -
- ]
2298 : : ! * Calculate (1-chi0*Vc-chi0*Kxc) and put it in chitmp.
2299 [ + + ]: 7232 : do ig1=1,npwe
2300 [ + + ]: 112512 : do ig2=1,npwe
2301 : 112512 : chitmp(ig1,ig2)=-chitmp(ig1,ig2)-chi0(ig1,ig2)*vc_sqrt(ig2)**2
2302 : : end do
2303 : 7232 : chitmp(ig1,ig1)=chitmp(ig1,ig1)+one
2304 : : end do
2305 : :
2306 : : ! * Invert (1-chi0*Vc-chi0*Kxc) and Multiply by chi0.
2307 : 443 : call xginv(chitmp,npwe,comm=comm)
2308 [ + - + - : 7102592 : chitmp=MATMUL(chitmp,chi0)
+ + + + +
- + - + -
- + - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- + + + +
+ + + + +
+ + - + -
+ - + - -
- - - - -
+ - ]
2309 : :
2310 : : ! * Save result, now chi0 contains chi.
2311 [ + + + + ]: 112955 : chi0=chitmp
2312 : :
2313 : 874 : select case (option_test)
2314 : : case (0)
2315 : : ! Symmetrized TESTPARTICLE epsilon^-1
2316 : 431 : call wrtout(std_out,' Calculating TESTPARTICLE epsilon^-1(G,G") = 1 + Vc*chi')
2317 [ + + ]: 6896 : do ig1=1,npwe
2318 [ + + ]: 103440 : chi0(ig1,:)=(vc_sqrt(ig1)*vc_sqrt(:))*chi0(ig1,:)
2319 : 6896 : chi0(ig1,ig1)=one+chi0(ig1,ig1)
2320 : : end do
2321 : :
2322 : : case (1)
2323 : : ! Symmetrized TESTELECTRON epsilon^-1
2324 : 12 : call wrtout(std_out,' Calculating TESTELECTRON epsilon^-1(G,G") = 1 + (Vc + fxc)*chi')
2325 [ + - + - : 490560 : chitmp=MATMUL(kxcg_mat,chi0)
+ - - + -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - + +
+ + + + +
+ + + + -
+ - - - -
- - - ]
2326 : :
2327 : : ! Perform hermitianization, only valid along the imaginary axis.
2328 [ + - ]: 12 : if (.not. ABS(REAL(omega))> tol3) call hermitianize(chitmp,"All")
2329 : :
2330 [ + + ]: 336 : do ig1=1,npwe
2331 [ + + ]: 9072 : chi0(ig1,:)=(vc_sqrt(ig1)*vc_sqrt(:))*chi0(ig1,:)+chitmp(ig1,:)
2332 : 336 : chi0(ig1,ig1)=one+chi0(ig1,ig1)
2333 : : end do
2334 : :
2335 : : case default
2336 [ + + - ]: 443 : ABI_BUG(sjoin('Wrong option_test:',itoa(option_test)))
2337 : : end select
2338 : :
2339 [ - + ]: 443 : ABI_FREE(chitmp)
2340 : : !
2341 : : ! === chi0 now contains symmetrical epsm1 ===
2342 : : ! * Calculate macroscopic dielectric constant epsm_lf(w)=1/epsm1(G=0,Gp=0,w).
2343 : 443 : epsm_lf(1) = one/chi0(1,1)
2344 : 443 : eelf (1) = -AIMAG(chi0(1,1))
2345 : :
2346 : : if (prtvol > 0) then
2347 : : write(msg,'(a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at omega',omega*Ha_eV,' [eV]'
2348 : : call wrtout(std_out,msg)
2349 : : call print_arr([std_out], chi0)
2350 : : end if
2351 : :
2352 [ + - ]: 443 : end subroutine atddft_symepsm1
2353 : : !!***
2354 : :
2355 : : !----------------------------------------------------------------------
2356 : :
2357 : : !!****f* m_screening/atddft_hyb_symepsm1
2358 : : !! NAME
2359 : : !! atddft_hyb_symepsm1
2360 : : !!
2361 : : !! FUNCTION
2362 : : !! Calculate $\tilde\epsilon^{-1}$ using ALDA within TDDFT
2363 : : !!
2364 : : !! Based on atddft_symepsm1, modified for the LR+ALDA hybrid kernel
2365 : : !!
2366 : : !! Output the electron energy loss function and the macroscopic dielectric function with and
2367 : : !! without local field effects (only if non-zero real frequencies are available)
2368 : : !!
2369 : : !! INPUTS
2370 : : !! iq_ibz=index of the q-point in the array Vcp%qibz where epsilon^-1 has to be calculated
2371 : : !! npwe=Number of G-vectors in chi0.
2372 : : !! nI,nJ=Number of rows/columns in chi0_ij (1,1 in collinear case)
2373 : : !! dim_wing=Dimension of the wings (0 or 3 if q-->0)
2374 : : !! kxcg_mat_sr=Short-range fxc kernel used in the test-electron epsilon^-1
2375 : : !! option_test=Only for TDDFT:
2376 : : !! == 0 for TESTPARTICLE ==
2377 : : !! == 1 for TESTELECTRON ==
2378 : : !! Vcp<vcoul_t>=Structure gathering data on the Coulombian interaction
2379 : : !! %nqibz=Number of q-points.
2380 : : !! %qibz(3,nqibz)=q-points in the IBZ.
2381 : : !! comm=MPI communicator.
2382 : : !! chi0_lwing(npwe*nI,dim_wing)=Lower wings of chi0 (only for q-->0)
2383 : : !! chi0_uwing(npwe*nJ,dim_wing)=Upper wings of chi0 (only for q-->0)
2384 : : !! chi0_head(dim_wing,dim_wing)=Head of of chi0 (only for q-->0)
2385 : : !!
2386 : : !! OUTPUT
2387 : : !!
2388 : : !! SIDE EFFECTS
2389 : : !! chi0(npwe*nI,npwe*nJ): in input the irreducible polarizability, in output
2390 : : !! the symmetrized inverse dielectric matrix.
2391 : : !!
2392 : : !! SOURCE
2393 : :
2394 : 40 : subroutine atddft_hyb_symepsm1(iq_ibz,Vcp,npwe,nI,nJ,chi0,kxcg_mat,kxcg_mat_sr,option_test,my_nqlwl,dim_wing,omega,&
2395 : 40 : & chi0_head,chi0_lwing,chi0_uwing,epsm_lf,epsm_nlf,eelf,comm)
2396 : :
2397 : : !Arguments ------------------------------------
2398 : : !scalars
2399 : : integer,intent(in) :: iq_ibz,nI,nJ,npwe,dim_wing,my_nqlwl
2400 : : integer,intent(in) :: option_test,comm
2401 : : type(vcoul_t),target,intent(in) :: Vcp
2402 : : !arrays
2403 : : complex(gwp),intent(in) :: kxcg_mat(npwe,npwe), kxcg_mat_sr(npwe,npwe)
2404 : : complex(dp),intent(in) :: omega
2405 : : complex(dp),intent(inout) :: chi0_lwing(npwe*nI,dim_wing)
2406 : : complex(dp),intent(inout) :: chi0_uwing(npwe*nJ,dim_wing)
2407 : : complex(dp),intent(inout) :: chi0_head(dim_wing,dim_wing)
2408 : : complex(gwp),intent(inout) :: chi0(npwe*nI,npwe*nJ)
2409 : : real(dp),intent(out) :: eelf(my_nqlwl)
2410 : : complex(dp),intent(out) :: epsm_lf(my_nqlwl),epsm_nlf(my_nqlwl)
2411 : :
2412 : : !Local variables-------------------------------
2413 : : !scalars
2414 : : integer,parameter :: master=0,prtvol=0
2415 : : integer :: ig1,ig2,my_rank,nprocs,ierr
2416 : : real(dp) :: ucvol, Zr
2417 : : logical :: is_qeq0
2418 : : character(len=500) :: msg
2419 : : !arrays
2420 : : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
2421 : 40 : complex(gwp),allocatable :: chitmp(:,:)
2422 : 40 : complex(gwp), contiguous, pointer :: vc_sqrt(:)
2423 : : ! *************************************************************************
2424 : :
2425 : : ABI_UNUSED(chi0_head(1,1))
2426 : : ABI_UNUSED(chi0_lwing(1,1))
2427 : : ABI_UNUSED(chi0_uwing(1,1))
2428 : :
2429 [ - + ]: 40 : if (nI/=1.or.nJ/=1) then
2430 : 0 : ABI_ERROR("nI or nJ=/1 not yet implemented")
2431 : : end if
2432 : :
2433 [ - + ]: 40 : ABI_CHECK(Vcp%nqlwl==1,"nqlwl/=1 not coded")
2434 [ - + ]: 40 : ABI_CHECK(my_nqlwl==1,"my_nqlwl/=1 not coded")
2435 : :
2436 : 40 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
2437 : :
2438 : 40 : call metric(gmet,gprimd,-1,rmet,Vcp%rprimd,ucvol)
2439 : :
2440 : 40 : is_qeq0 = (normv(Vcp%qibz(:,iq_ibz),gmet,'G')<GW_TOLQ0)
2441 : :
2442 [ + + ]: 40 : if (iq_ibz==1) then
2443 : : !%vc_sqrt => Vcp%vcqlwl_sqrt(:,iqlwl) ! Use Coulomb term for q-->0
2444 : 10 : vc_sqrt => Vcp%vcqlwl_sqrt(:,1) ! TODO add treatment of non-Analytic behavior
2445 : : else
2446 : 30 : vc_sqrt => Vcp%vc_sqrt(:,iq_ibz)
2447 : : end if
2448 : :
2449 : 40 : write(msg,'(a,f8.2,a)')" chitmp requires: ",npwe**2*gwp*b2Mb," Mb"
2450 [ - + + - : 160 : ABI_MALLOC_OR_DIE(chitmp,(npwe,npwe), ierr)
- + + - +
- + - + -
- + - + ]
2451 : : !
2452 : : ! * Calculate chi0*fxc.
2453 [ + - + - : 289280 : chitmp = MATMUL(chi0,kxcg_mat)
+ - - + -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - + +
+ + + + +
+ + + + -
+ - - - -
- - - ]
2454 : : ! * First calculate the NLF contribution
2455 [ + + ]: 640 : do ig1=1,npwe
2456 [ + + ]: 9600 : do ig2=1,npwe
2457 : 9600 : chitmp(ig1,ig2)=-chitmp(ig1,ig2)
2458 : : end do
2459 : 640 : chitmp(ig1,ig1)=chitmp(ig1,ig1)+one
2460 : : end do
2461 : :
2462 : 40 : call xginv(chitmp,npwe,comm=comm)
2463 : :
2464 [ + - + - : 568960 : chitmp = MATMUL(chitmp,chi0)
+ + + + +
- + - + -
- + - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- + + + +
+ + + + +
+ + - + -
+ - + - -
- - - - -
+ - ]
2465 : : !if (.not. ABS(REAL(omega))> tol3) call hermitianize(chitmp,"All")
2466 : 40 : chitmp(1,1)=-vc_sqrt(1)*chitmp(1,1)*vc_sqrt(1)
2467 : 40 : chitmp(1,1)=chitmp(1,1)+one
2468 : :
2469 : 40 : epsm_nlf(1)=chitmp(1,1)
2470 : :
2471 [ + - + - : 289240 : chitmp = MATMUL(chi0,kxcg_mat)
- + - - -
- - - - -
- - - - -
- - - + +
+ + + + +
+ + + + -
+ - - - -
- ]
2472 : : ! * Calculate (1-chi0*Vc-chi0*Kxc) and put it in chitmp.
2473 [ + + ]: 640 : do ig1=1,npwe
2474 [ + + ]: 9600 : do ig2=1,npwe
2475 : 9600 : chitmp(ig1,ig2)=-chitmp(ig1,ig2)-chi0(ig1,ig2)*vc_sqrt(ig2)**2
2476 : : end do
2477 : 640 : chitmp(ig1,ig1)=chitmp(ig1,ig1)+one
2478 : : end do
2479 : :
2480 : : ! * Invert (1-chi0*Vc-chi0*Kxc) and Multiply by chi0.
2481 : 40 : call xginv(chitmp,npwe,comm=comm)
2482 [ + - + - : 568960 : chitmp=MATMUL(chitmp,chi0)
+ + + + +
- + - + -
- + - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- + + + +
+ + + + +
+ + - + -
+ - + - -
- - - - -
+ - ]
2483 : :
2484 : : ! * Save result, now chi0 contains chi.
2485 [ + + + + ]: 9640 : chi0=chitmp
2486 : :
2487 : 40 : select case (option_test)
2488 : : case (0)
2489 : : ! Symmetrized TESTPARTICLE epsilon^-1
2490 : 0 : call wrtout(std_out,' Calculating TESTPARTICLE epsilon^-1(G,G") = 1 + Vc*chi')
2491 [ # # ]: 0 : do ig1=1,npwe
2492 [ # # ]: 0 : chi0(ig1,:)=(vc_sqrt(ig1)*vc_sqrt(:))*chi0(ig1,:)
2493 : 0 : chi0(ig1,ig1)=one+chi0(ig1,ig1)
2494 : : end do
2495 : :
2496 : : case (1)
2497 : : ! Symmetrized TESTELECTRON epsilon^-1
2498 : 40 : call wrtout(std_out,' Calculating TESTELECTRON epsilon^-1(G,G") = 1 + Vc*chi + Zr*Kxc_sr*chi')
2499 [ + - + - : 289280 : chitmp=MATMUL(kxcg_mat_sr,chi0)
+ - - + -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - + +
+ + + + +
+ + + + -
+ - - - -
- - - ]
2500 : 40 : Zr = 0.78
2501 : :
2502 : : ! Perform hermitianization, only valid along the imaginary axis.
2503 [ + + ]: 40 : if (.not. ABS(REAL(omega))> tol3) call hermitianize(chitmp,"All")
2504 : :
2505 [ + + ]: 640 : do ig1=1,npwe
2506 [ + + ]: 9600 : chi0(ig1,:)=(vc_sqrt(ig1)*vc_sqrt(:))*chi0(ig1,:)+chitmp(ig1,:)*Zr
2507 : 640 : chi0(ig1,ig1)=one+chi0(ig1,ig1)
2508 : : end do
2509 : :
2510 : : case default
2511 [ - + - ]: 40 : ABI_BUG(sjoin('Wrong option_test:',itoa(option_test)))
2512 : : end select
2513 : :
2514 [ - + ]: 40 : ABI_FREE(chitmp)
2515 : : !
2516 : : ! === chi0 now contains symmetrical epsm1 ===
2517 : : ! * Calculate macroscopic dielectric constant epsm_lf(w)=1/epsm1(G=0,Gp=0,w).
2518 : 40 : epsm_lf(1) = one/chi0(1,1)
2519 : 40 : eelf (1) = -AIMAG(chi0(1,1))
2520 : :
2521 : : if (prtvol > 0) then
2522 : : write(msg,'(a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at omega',omega*Ha_eV,' [eV]'
2523 : : call wrtout(std_out,msg)
2524 : : call print_arr([std_out], chi0)
2525 : : end if
2526 : :
2527 [ + - ]: 40 : end subroutine atddft_hyb_symepsm1
2528 : : !!***
2529 : :
2530 : : !----------------------------------------------------------------------
2531 : :
2532 : : !!****f* m_screening/mkem1_q0
2533 : : !! NAME
2534 : : !! mkem1_q0
2535 : : !!
2536 : : !! FUNCTION
2537 : : !! This routine construct the microscopic dieletric matrix for q-->0 starting from the heads, wings and the body
2538 : : !! of the irreducible polarizability. Afterwards it calculates the symmetrized inverse dieletric matrix
2539 : : !! via a block wise inversion thus obtaining the heads and the wings of e^{-1} that can be
2540 : : !! used to describe the non-analytic behavior for q-->0.
2541 : : !!
2542 : : !! INPUTS
2543 : : !! npwe=Number of Gs used to describe chi0
2544 : : !! nomega=Number of frequencies in chi0.
2545 : : !! n1,n2=Factors used to define the same of the chi0 matrix (1,1 if collinear, the typical case)
2546 : : !! Cryst<crystal_t>=Structure describing the crystal structure.
2547 : : !! Vcp<vcoul_t>=datatypes gathering info on the Coulomb term
2548 : : !! gvec(3,npwe)=G-vector for chi0 in reduced coordinates.
2549 : : !! comm=MPI communicator
2550 : : !!
2551 : : !! OUTPUT
2552 : : !! eps_head(3,3,nomega)=The macroscopic dieletric tensor in reduced coordinates.
2553 : : !! The dieletric matrix along versor \hat q can be obtained with
2554 : : !! e(\hat q) = \hat q.eps_head \hat q if all quantities are given in Cartesian coordinates.
2555 : : !!
2556 : : !! SIDE EFFECTS
2557 : : !! chi0(npwe*n1,npwe*n2,nomega)= Input: polarizability. output: inverse dieletric matrix (only the body is used)
2558 : : !! chi0_lwing(npwe*n1,nomega,3)
2559 : : !! chi0_uwing(npwe*n2,nomega,3) Input: the lower and upper wings of the polarizability
2560 : : !! Output: the "lower" and "upper" wings of the inverse dieletric matrix. See notes below.
2561 : : !! chi0_head(3,3,nomega)= Input: the polarizability tensor in Cartesian coordinates.
2562 : : !! Output: The "head" of the inverse dieletric matrix. See notes below.
2563 : : !!
2564 : : !! NOTES
2565 : : !! Matrix inversion in block form.
2566 : : !!
2567 : : !! 1 n-1
2568 : : !! M = | c u^t| 1 ==> M^{-1} = | 1/k -u^t A^{-1}/k |
2569 : : !! | v A | n-1 | -A^{-1} v/k A^{-1} + (A^{-1}v u^t A^{-1})/k |
2570 : : !!
2571 : : !! where k = c - u^t A^{-1} v
2572 : : !!
2573 : : !! Let q be a versor in reciprocal space, the symmetrized dielectric matrix with bare coulomb interaction
2574 : : !! can be written as
2575 : : !!
2576 : : !! \tilde\epsilon = | q.Eq q.Wl(G2) | where E_ij = \delta_ij -4\pi chi0_head_ij
2577 : : !! | q.Wl(G1) B(G1,G2 | Wl(G1) = -4\pi chi0_lwing(G1)
2578 : : !! Wu(G2) = -4\pi chi0_uwing(G1)
2579 : : !! therefore, in Cartesian coordinates, we have:
2580 : : !!
2581 : : !! 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)
2582 : : !!
2583 : : !! 2) e^{-1}_{0G'}(q) = -e^{-1}_{00}(q) [ \sum_{iG} q_i Wu_i(G)a B_{GG'}^{-1} ] = (q.Su) /(q.Lq)
2584 : : !!
2585 : : !! 3) e^{-1}_{G0}(q) = -e^{-1}_{00}(q) [ \sum_{iG'} q_i B_{GG'}^{-1} Wl_i(G') ] = (q.Sl) /(q.Lq)
2586 : : !!
2587 : : !! 4) e^{-1}_{GG'}(q) = B_{GG'}^{-1} +
2588 : : !! [ \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)
2589 : : !!
2590 : : !! where Su(G,3) and Sl(G,3) are the "upper" and "lower" wings of the inverse dielectric matrix and
2591 : : !! L is the inverse dielectric tensor. Similar equations hold even if vectors and tensors are given in terms
2592 : : !! of the reciprocal lattice vectors provided that the metric tensor is taken into account.
2593 : : !! The main difference is in the expression for the tensor as only one metric tensor can be
2594 : : !! absorbed in the scalar product, the second metric multiplies one of the wings.
2595 : : !!
2596 : : !! *) The present implementation assumes that no cutoff technique is used in the Coulomb term.
2597 : : !!
2598 : : !! *) Once the tensor in know it is possible to average the quadratic form on the sphere exactly.
2599 : : !! In Cartesian coordinates one obtains.
2600 : : !!
2601 : : !! \dfrac{1}{4\pi} \int v.Tv d\Omega = Trace(T)/3
2602 : : !!
2603 : : !! For the inverse dielectric matrix we have to resort to a numerical integration
2604 : : !!
2605 : : !! SOURCE
2606 : :
2607 : 0 : subroutine mkem1_q0(npwe,n1,n2,nomega,Cryst,Vcp,gvec,chi0_head,chi0_lwing,chi0_uwing,chi0,eps_head,comm)
2608 : :
2609 : : !Arguments ------------------------------------
2610 : : !scalars
2611 : : integer,intent(in) :: npwe,nomega,n1,n2,comm
2612 : : type(crystal_t),intent(in) :: Cryst
2613 : : type(vcoul_t),intent(in) :: Vcp
2614 : : !arrays
2615 : : integer,intent(in) :: gvec(3,npwe)
2616 : : complex(gwp),intent(in) :: chi0(npwe*n1,npwe*n2,nomega)
2617 : : complex(dp),intent(inout) :: chi0_lwing(npwe*n1,nomega,3)
2618 : : complex(dp),intent(inout) :: chi0_uwing(npwe*n2,nomega,3)
2619 : : complex(dp),intent(inout) :: chi0_head(3,3,nomega)
2620 : : complex(dp),intent(out) :: eps_head(3,3,nomega)
2621 : :
2622 : : !Local variables ------------------------------
2623 : : !scalars
2624 : : integer :: iw,ig,ig1,ig2,idir,jdir
2625 : : !arrays
2626 : 0 : real(dp),allocatable :: modg_inv(:)
2627 : 0 : complex(dp),allocatable :: eps_lwing(:,:),eps_uwing(:,:),eps_body(:,:),cvec(:)
2628 : : !************************************************************************
2629 : :
2630 [ # # ]: 0 : ABI_CHECK(npwe /= 1, "npwe must be >1")
2631 : : ABI_UNUSED(comm)
2632 : :
2633 : : ! Precompute 1/|G|.
2634 [ # # # # : 0 : ABI_MALLOC(modg_inv,(npwe-1))
# # # # #
# # # #
# ]
2635 [ # # ]: 0 : do ig=1,npwe-1
2636 : 0 : modg_inv(ig) = one/normv(gvec(:,ig+1),Cryst%gmet,'G')
2637 : : end do
2638 : :
2639 [ # # # # : 0 : ABI_MALLOC(eps_uwing,((npwe-1)*n1,3))
# # # # #
# # # ]
2640 [ # # # # : 0 : ABI_MALLOC(eps_lwing,((npwe-1)*n2,3))
# # # # #
# # # ]
2641 [ # # # # : 0 : ABI_MALLOC(eps_body,(npwe-1,npwe-1))
# # # # #
# # # # #
# # ]
2642 [ # # # # : 0 : ABI_MALLOC(cvec,(npwe-1))
# # # # #
# # # ]
2643 : :
2644 [ # # ]: 0 : do iw=1,nomega
2645 : : !
2646 : : ! Head and wings of the symmetrized epsilon.
2647 [ # # # # ]: 0 : eps_head(:,:,iw) = -four_pi*chi0_head(:,:,iw)
2648 [ # # ]: 0 : do idir=1,3
2649 : 0 : eps_head(idir,idir,iw) = one + eps_head(idir,idir,iw)
2650 [ # # ]: 0 : eps_lwing(:,idir) = -four_pi * modg_inv * chi0_lwing(2:,iw,idir)
2651 [ # # ]: 0 : eps_uwing(:,idir) = -four_pi * modg_inv * chi0_uwing(2:,iw,idir)
2652 : : !eps_lwing(:,idir) = -chi0_lwing(2:,iw,idir) * SQRT(four_pi) * Vcp%vcqlwl_sqrt(2:npwe,1)
2653 : : !eps_uwing(:,idir) = -chi0_uwing(2:,iw,idir) * SQRT(four_pi) * Vcp%vcqlwl_sqrt(2:npwe,1)
2654 : : end do
2655 : :
2656 : 0 : write(std_out,*)" espilon head"
2657 [ # # ]: 0 : call print_arr([std_out], eps_head(:,:,iw))
2658 : : !
2659 : : ! Construct the body of the symmetrized epsilon then invert it.
2660 [ # # ]: 0 : do ig2=1,npwe-1
2661 [ # # ]: 0 : do ig1=1,npwe-1
2662 : 0 : eps_body(ig1,ig2) = -four_pi * modg_inv(ig1)*chi0(ig1+1,ig2+1,iw )*modg_inv(ig2)
2663 : : !eps_body(ig1,ig2) = -Vcp%vcqlwl_sqrt(ig1+1,1)*chi0(ig1+1,ig2+1,iw)* Vcp%vcqlwl_sqrt(ig2+1,1)
2664 : : end do
2665 : 0 : eps_body(ig2,ig2) = one + eps_body(ig2,ig2)
2666 : : end do
2667 : :
2668 : 0 : call xginv(eps_body,npwe-1)
2669 : : !
2670 : : ! Overwrite chi0_head and chi0_wings with the head and the wings of the inverse dielectric matrix.
2671 [ # # ]: 0 : do jdir=1,3
2672 : : !
2673 : : ! Head.
2674 [ # # ]: 0 : cvec=czero
2675 [ # # ]: 0 : do idir=1,3
2676 [ # # # # : 0 : cvec = cvec + two_pi*Cryst%gmet(jdir,idir)*MATMUL(eps_body,eps_lwing(:,idir)) ! as we work in reciprocal coords.
# # # # #
# # # # #
# # # # #
# # # ]
2677 : : end do
2678 : : !cvec = MATMUL(eps_body,eps_lwing(:,jdir))
2679 [ # # ]: 0 : do idir=1,3
2680 : 0 : chi0_head(idir,jdir,iw) = eps_head(idir,jdir,iw) - xdotu(npwe-1,eps_uwing(:,idir),1,cvec,1)
2681 : : end do
2682 : : !
2683 : : ! Now the wings.
2684 [ # # # # : 0 : chi0_uwing(2:,iw,jdir) = -MATMUL(eps_uwing(:,jdir),eps_body)
# # # # #
# # # # #
# # # # #
# # # ]
2685 [ # # # # : 0 : chi0_lwing(2:,iw,jdir) = -MATMUL(eps_body,eps_lwing(:,jdir))
# # # # #
# # # # #
# # # # #
# ]
2686 : : !
2687 : : end do !jdir
2688 : :
2689 : 0 : call wrtout(std_out, "espilon^1 head after block inversion")
2690 [ # # ]: 0 : call print_arr([std_out], chi0_head(:,:,iw))
2691 : : !
2692 : : ! Change the body but do not add the corrections due to the head and the wings.
2693 : : ! since they can be obtained on the fly from eps_body and the wings of eps^{-1}.
2694 : : !%chi0(2:,2:,iw) = eps_body
2695 : : end do !iw
2696 : :
2697 [ # # ]: 0 : ABI_FREE(modg_inv)
2698 [ # # ]: 0 : ABI_FREE(cvec)
2699 [ # # ]: 0 : ABI_FREE(eps_lwing)
2700 [ # # ]: 0 : ABI_FREE(eps_uwing)
2701 [ # # ]: 0 : ABI_FREE(eps_body)
2702 : :
2703 : : RETURN
2704 : : ABI_UNUSED(Vcp%ng)
2705 : :
2706 : : end subroutine mkem1_q0
2707 : : !!***
2708 : :
2709 : : !----------------------------------------------------------------------
2710 : :
2711 : : !!****f* m_screening/lebedev_laikov_int
2712 : : !! NAME
2713 : : !! lebedev_laikov_int
2714 : : !!
2715 : : !! FUNCTION
2716 : : !!
2717 : : !! INPUTS
2718 : : !!
2719 : : !! OUTPUT
2720 : : !!
2721 : : !! SOURCE
2722 : :
2723 : : subroutine lebedev_laikov_int()
2724 : :
2725 : : !Arguments ------------------------------------
2726 : :
2727 : : !Local variables-------------------------------
2728 : : !scalars
2729 : : integer :: on,npts,ii,ll,mm,lmax,leb_idx !ierr,
2730 : : real(dp) :: accuracy
2731 : : complex(dp) :: ang_int
2732 : : !arrays
2733 : : real(dp) :: cart_vpt(3) !,real_pars(0)
2734 : : real(dp),allocatable :: vx(:),vy(:),vz(:),ww(:)
2735 : : complex(dp) :: tensor(3,3),cplx_pars(9)
2736 : : complex(dp),allocatable :: ref_func(:),expd_func(:) !tmp_momenta(:)
2737 : : ! *************************************************************************
2738 : :
2739 : : ABI_ERROR("lebedev_laikov_int is still under development")
2740 : :
2741 : : !tensor=RESHAPE((/4.0,2.0,4.0,0.5,2.1,0.0,5.4,2.1,5.0/),(/3,3/))
2742 : : tensor=RESHAPE((/4.0,0.0,0.0,0.0,4.0,0.0,0.0,0.0,5.0/),(/3,3/))
2743 : : !tensor=RESHAPE((/1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0/),(/3,3/))
2744 : :
2745 : : npts=26
2746 : : ABI_MALLOC(vx,(npts))
2747 : : ABI_MALLOC(vy,(npts))
2748 : : ABI_MALLOC(vz,(npts))
2749 : : ABI_MALLOC(ww,(npts))
2750 : :
2751 : : !call LD0026(vx,vy,vz,ww,on)
2752 : :
2753 : : ang_int=czero
2754 : : do ii=1,npts
2755 : : cart_vpt = [vx(ii),vy(ii),vz(ii)]
2756 : : ang_int = ang_int + ww(ii)*DOT_PRODUCT(cart_vpt,MATMUL(tensor,cart_vpt))
2757 : : end do
2758 : :
2759 : : !write(std_out,*)"quadratic form associated to tensor=",tensor
2760 : : write(std_out,*)"on ang_int",on,ang_int
2761 : :
2762 : : ABI_FREE(vx)
2763 : : ABI_FREE(vy)
2764 : : ABI_FREE(vz)
2765 : : ABI_FREE(ww)
2766 : :
2767 : : !call init_lebedev_gridset()
2768 : : cplx_pars = RESHAPE(tensor,(/9/)); accuracy=tol10
2769 : :
2770 : : ! This is the function to be expanded evaluated on the lebedev_laikov grid of index leb_idx
2771 : : leb_idx=3; npts=lebedev_npts(leb_idx)
2772 : : ABI_MALLOC(ref_func,(npts))
2773 : : do ii=1,npts
2774 : : !cart_vpt = Lgridset(leb_idx)%versor(:,ii)
2775 : : ref_func(ii) = one/DOT_PRODUCT(cart_vpt,MATMUL(tensor,cart_vpt))
2776 : : end do
2777 : :
2778 : : ! Calculate the expansion in angular momenta of 1/{q.Tq}.
2779 : : ! Only even l-components contribute thanks to the parity of the integrand.
2780 : : ! tol6 seems to be an acceptable error, convergence wrt lmax is very slow even for simple tensors.
2781 : : ABI_MALLOC(expd_func,(npts))
2782 : : expd_func=czero
2783 : : lmax=10
2784 : : do ll=0,lmax,2
2785 : : !allocate(tmp_momenta(-ll:ll))
2786 : : do mm=-ll,ll
2787 : : ! MG: Commented because it causes problems with the new version of abilint
2788 : : !call lebedev_quadrature(ylmstar_over_qTq,(/ll,mm/),real_pars,cplx_pars,ang_int,ierr,accuracy)
2789 : : write(std_out,*)ll,mm,ang_int
2790 : : !tmp_momenta(mm) = ang_int
2791 : : do ii=1,npts
2792 : : !cart_vpt = Lgridset(leb_idx)%versor(:,ii)
2793 : : expd_func(ii) = expd_func(ii) + four_pi*ang_int*ylmc(ll,mm,cart_vpt)
2794 : : end do
2795 : : end do
2796 : : !deallocate(tmp_momenta)
2797 : : write(std_out,*)"Error in angular expansion at l=",ll," is ",MAXVAL(ABS(expd_func-ref_func))
2798 : : end do
2799 : :
2800 : : !BEGINDEBUG
2801 : : ! do ii=1,npts
2802 : : ! write(777,*)ref_func(ii)
2803 : : ! write(778,*)expd_func(ii)
2804 : : ! end do
2805 : : !ENDDEBUG
2806 : :
2807 : : ABI_FREE(expd_func)
2808 : : ABI_FREE(ref_func)
2809 : :
2810 : : ABI_ERROR("Exiting from lebedev_laikov_int")
2811 : :
2812 : : end subroutine lebedev_laikov_int
2813 : : !!***
2814 : :
2815 : : !----------------------------------------------------------------------
2816 : :
2817 : : !!****f* m_screening/ylmstar_over_qTq
2818 : : !! NAME
2819 : : !! ylmstar_over_qTq
2820 : : !!
2821 : : !! FUNCTION
2822 : : !! Return Ylm(q)^*/(q.Tq) where q is a versor in Cartesian coordinates.
2823 : : !! and Ylm is a complex spherical Harmonics whose index (l,m) are
2824 : : !! passed via int_pars(1:2). T is a tensore in Cartesian coordinates
2825 : : !! passed via cplx_pars(1:9).
2826 : : !!
2827 : : !! INPUTS
2828 : : !! cart_vers(3)=Cartesian components of the versor
2829 : : !! int_pars(1:2)=(l,m) indices in Ylm. l>=0 and m \in [-l,l]
2830 : : !! cplx_pars(1:9)=Tensor T in Cartesian coordinates.
2831 : : !! real_pars=Not used.
2832 : : !!
2833 : : !! OUTPUT
2834 : : !! Value of Ylm(q)^*/(q.Tq)
2835 : : !!
2836 : : !! SOURCE
2837 : :
2838 : : function ylmstar_over_qTq(cart_vers,int_pars,real_pars,cplx_pars)
2839 : :
2840 : : !Arguments ------------------------------------
2841 : : !scalars
2842 : : real(dp),intent(in) :: cart_vers(3)
2843 : : integer,intent(in) :: int_pars(:)
2844 : : real(dp),intent(in) :: real_pars(:)
2845 : : complex(dp),intent(in) :: cplx_pars(:)
2846 : : complex(dp) :: ylmstar_over_qTq
2847 : : !arrays
2848 : :
2849 : : !Local variables-------------------------------
2850 : : !scalars
2851 : : integer :: ll,mm
2852 : : !arrays
2853 : : complex(dp) :: tensor(3,3)
2854 : : ! *************************************************************************
2855 : :
2856 : : tensor = RESHAPE(cplx_pars(1:9),(/3,3/))
2857 : : ll = int_pars(1) ! ll starts from zero.
2858 : : mm = int_pars(2) ! m \in [-l,l]
2859 : :
2860 : : ylmstar_over_qTq = CONJG(ylmc(ll,mm,cart_vers))/DOT_PRODUCT(cart_vers,MATMUL(tensor,cart_vers))
2861 : :
2862 : : RETURN
2863 : : ABI_UNUSED(real_pars(1))
2864 : :
2865 : : end function ylmstar_over_qTq
2866 : : !!***
2867 : :
2868 : : !----------------------------------------------------------------------
2869 : :
2870 : : !!****f* m_screening/ylmstar_wtq_over_qTq
2871 : : !! NAME
2872 : : !! ylmstar_wtq_over_qTq
2873 : : !!
2874 : : !! FUNCTION
2875 : : !! Return Ylm(q)^* weight(q)/(q.Tq) where q is a versor in Cartesian coordinates.
2876 : : !! Ylm is a complex spherical Harmonics whose index (l,m) are
2877 : : !! passed via int_pars(1:2). T is a tensor in Cartesian coordinates
2878 : : !! passed via cplx_pars(1:9). weight(q) is the weighting function giving
2879 : : !! the length of the vector parallel to versor q that connects the origin
2880 : : !! of the lattice to one of the boundaries of the small cell centered at Gamma
2881 : : !!
2882 : : !! INPUTS
2883 : : !! cart_vers(3)=Cartesian components of the versor
2884 : : !! int_pars(1:2)=(l,m) indices in Ylm. l>=0 and m \in [-l,l]
2885 : : !! cplx_pars(1:9)=Tensor T in Cartesian coordinates.
2886 : : !! real_pars(1:9)=The Cartesian vectors defining the small box centered around gamma point
2887 : : !! when referred to this vectors the points in the box are given by {(x,y,z) | x,y,z \in [-1,1]}.
2888 : : !!
2889 : : !! OUTPUT
2890 : : !! Value of Ylm(q)^* weigh(q)/(q.Tq)
2891 : : !!
2892 : : !! SOURCE
2893 : :
2894 : : function ylmstar_wtq_over_qTq(cart_vers,int_pars,real_pars,cplx_pars)
2895 : :
2896 : : !Arguments ------------------------------------
2897 : : !scalars
2898 : : real(dp),intent(in) :: cart_vers(3)
2899 : : integer,intent(in) :: int_pars(:)
2900 : : real(dp),intent(in) :: real_pars(:)
2901 : : complex(dp),intent(in) :: cplx_pars(:)
2902 : : complex(dp) :: ylmstar_wtq_over_qTq
2903 : : !arrays
2904 : :
2905 : : !Local variables-------------------------------
2906 : : !scalars
2907 : : integer :: ll,mm
2908 : : real(dp) :: wtq
2909 : : !arrays
2910 : : real(dp) :: gprimd(3,3),rprimd(3,3),red_vers(3)
2911 : : complex(dp) :: tensor(3,3)
2912 : : ! *************************************************************************
2913 : :
2914 : : ABI_ERROR("Work in progress")
2915 : : ! box_len has to be tested
2916 : :
2917 : : gprimd = RESHAPE(real_pars(1:9),(/3,3/))
2918 : : red_vers = MATMUL(rprimd,cart_vers)
2919 : : wtq = box_len(red_vers,gprimd)
2920 : :
2921 : : tensor = RESHAPE(cplx_pars(1:9),(/3,3/))
2922 : : ll = int_pars(1) ! true ll i.e. not shifted
2923 : : mm = int_pars(2)
2924 : :
2925 : : ylmstar_wtq_over_qTq = CONJG(ylmc(ll,mm,cart_vers))*wtq/DOT_PRODUCT(cart_vers,MATMUL(tensor,cart_vers))
2926 : :
2927 : : end function ylmstar_wtq_over_qTq
2928 : : !!***
2929 : :
2930 : : !----------------------------------------------------------------------
2931 : :
2932 : : !!****f* m_screening/mdielf_bechstedt
2933 : : !! NAME
2934 : : !! mdielf_bechstedt
2935 : : !!
2936 : : !! FUNCTION
2937 : : !! Calculates the model dielectric function for the homogeneous system
2938 : : !! as proposed by F. Bechstedt, in Solid State Commun. 84, 765 1992.
2939 : : !!
2940 : : !! INPUTS
2941 : : !! eps_inf=Dielectric constant of the material
2942 : : !! qnrm=The modulus of the q-point.
2943 : : !! rhor=The local value of the density
2944 : : !!
2945 : : !! SOURCE
2946 : :
2947 : 140216832 : elemental function mdielf_bechstedt(eps_inf, qnrm, rhor) result(mdielf)
2948 : :
2949 : : !Arguments ------------------------------------
2950 : : !scalars
2951 : : real(dp),intent(in) :: eps_inf,qnrm,rhor
2952 : : real(dp) :: mdielf
2953 : : ! *************************************************************************
2954 : :
2955 : : mdielf = one + &
2956 : 140216832 : one / ( one/(eps_inf-one) + (qnrm/k_thfermi(rhor))**2 + (three*qnrm**4)/(four*k_fermi(rhor)**2 * k_thfermi(rhor)**2) )
2957 : :
2958 : 140216832 : end function mdielf_bechstedt
2959 : : !!***
2960 : :
2961 : : !----------------------------------------------------------------------
2962 : :
2963 : : !!****f* m_screening/screen_mdielf
2964 : : !! NAME
2965 : : !! screen_mdielf
2966 : : !!
2967 : : !! FUNCTION
2968 : : !! Calculates W_{G,G'}(q,w) for a given q-point in the BZ using a model dielectric function.
2969 : : !!
2970 : : !! INPUTS
2971 : : !! iq_bz=The index of the q-point in the BZ where W(q) is calculated.
2972 : : !! npw=Number of plane waves for W
2973 : : !! nomega=Number of frequency points.
2974 : : !! model_type=Flag defining the model.
2975 : : !! eps_inf=Dielectric constant of the material.
2976 : : !! Cryst<crystal_t>=Info on the unit cell
2977 : : !! Qmesh<kmesh_t>=Info on the set of q-points.
2978 : : !! Vcp<vcoul_t datatype>= containing information on the cutoff technique
2979 : : !! Gsph<Gsphere>=The G-sphere for W.
2980 : : !! nspden=Number of spin density components of the density.
2981 : : !! nfft=Number of FFT points on the dense FFT mesh
2982 : : !! ngfft(18)=contain all needed information about 3D FFT.
2983 : : !! rhor(nfft,nspden)=Electron density in real space (The PAW AE term is included)
2984 : : !! which= Set it to "EM1" if the symmetrized inverse dielectric matrix is wanted.
2985 : : !! By default the routines returns W.
2986 : : !! comm=MPI communicator.
2987 : : !!
2988 : : !! OUTPUT
2989 : : !! w_qbz(npw,npw,nomega)
2990 : : !!
2991 : : !! NOTES
2992 : : !! W_{G1,G2} = 1/2 {
2993 : : !! v(q+G1) \int em1(|q+G1|,r) e^{-i(G1-G2).r} dr +
2994 : : !! v(q+G2) \int em1(|q+G2|,r) e^{-i(G1-G2).r} dr } / \Omega
2995 : : !!
2996 : : !! SOURCE
2997 : :
2998 : 631 : subroutine screen_mdielf(iq_bz,npw,nomega,model_type,eps_inf,Cryst,Qmesh,Vcp,Gsph,nspden,nfft,ngfft,rhor,which,w_qbz,comm)
2999 : :
3000 : : !Arguments ------------------------------------
3001 : : !scalars
3002 : : integer,intent(in) :: npw,nomega,nfft,nspden,iq_bz,comm,model_type
3003 : : real(dp),intent(in) :: eps_inf
3004 : : character(len=*),intent(in) :: which
3005 : : type(kmesh_t),intent(in) :: Qmesh
3006 : : type(crystal_t),intent(in) :: Cryst
3007 : : type(vcoul_t),target,intent(in) :: Vcp
3008 : : type(gsphere_t),intent(in) :: Gsph
3009 : : !arrays
3010 : : integer,intent(in) :: ngfft(18)
3011 : : real(dp),intent(in) :: rhor(nfft,nspden)
3012 : : complex(gwp),intent(out) :: w_qbz(npw,npw,nomega)
3013 : :
3014 : : !Local variables-------------------------------
3015 : : !scalars
3016 : : integer,parameter :: tim_fourdp0=0,paral_kgb0=0,cplex1=1
3017 : : integer :: my_gstart,my_gstop,iq_ibz,ig,itim_q,isym_q
3018 : : integer :: ig1,ig2,g1mg2_fft,iw,ii,ierr,nprocs,isg,ifft !,row ,col
3019 : : real(dp) :: qpg2_nrm
3020 : : complex(dp) :: ph_mqbzt
3021 : : logical :: is_qeq0,isirred
3022 : : !character(len=500) :: msg
3023 : 631 : type(MPI_type) :: MPI_enreg_seq
3024 : : !arrays
3025 : : integer :: umklp(3)
3026 : 631 : integer,allocatable :: igfft(:),g1mg2(:,:)
3027 : : real(dp) :: qpg2(3),qpt_bz(3)
3028 : 631 : real(dp),allocatable :: em1_qpg2r(:),fofg(:,:)
3029 : 631 : complex(gwp),contiguous, pointer :: vc_sqrt_ibz(:)
3030 : 631 : complex(gwp),allocatable :: vc_qbz(:),ctmp(:,:)
3031 : 631 : logical,allocatable :: mask(:)
3032 : : ! *************************************************************************
3033 : :
3034 : 0 : ABI_CHECK(nomega==1,"screen_mdielf does not support nomega>1")
3035 : :
3036 : : ! Fake MPI_type for the sequential part.
3037 : 631 : call initmpi_seq(MPI_enreg_seq)
3038 : 631 : call MPI_enreg_seq%distribfft%init_seq('c',ngfft(2),ngfft(3),'all')
3039 : :
3040 : 631 : nprocs = xmpi_comm_size(comm)
3041 : 631 : call xmpi_split_work(npw,comm,my_gstart,my_gstop)
3042 : :
3043 : 631 : call qmesh%get_bz_item(iq_bz,qpt_bz,iq_ibz,isym_q,itim_q,ph_mqbzt,umklp,isirred)
3044 : :
3045 : : !if (itim_q/=1.or.isym_q/=1.or.ANY(umklp/=0) ) then
3046 : : ! ABI_ERROR("Bug in mdielf_bechstedt")
3047 : : !end if
3048 : : !
3049 : : ! Symmetrize Vc in the full BZ.
3050 : 631 : is_qeq0 = (normv(qpt_bz,Cryst%gmet,'G')<GW_TOLQ0) ! Check if q==0
3051 [ + + ]: 631 : if (is_qeq0) then
3052 : 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?
3053 : : else
3054 : 582 : vc_sqrt_ibz => Vcp%vc_sqrt(:,iq_ibz)
3055 : : end if
3056 : :
3057 [ - + + - : 1893 : ABI_MALLOC(vc_qbz,(npw))
+ - + - -
+ - + ]
3058 [ + + ]: 33260 : do ig=1,npw
3059 : 32629 : isg = Gsph%rottb(ig,itim_q,isym_q)
3060 : 33260 : vc_qbz(isg) = vc_sqrt_ibz(ig)**2
3061 : : end do
3062 : :
3063 [ - + + - : 1893 : ABI_MALLOC(igfft,(npw))
+ - + - -
+ - + ]
3064 [ - + + - : 1893 : ABI_MALLOC(g1mg2,(3,npw))
+ - + - -
+ - + ]
3065 [ - + + - : 1893 : ABI_MALLOC(fofg,(2,nfft))
+ - + - -
+ - + ]
3066 [ - + + - : 1893 : ABI_MALLOC(em1_qpg2r,(nfft))
+ - + - -
+ - + ]
3067 [ - + + - : 1262 : ABI_MALLOC(mask,(npw))
- + - + ]
3068 : :
3069 [ + + + + : 1724402 : w_qbz=czero
+ + ]
3070 [ + + ]: 31372 : do ig2=my_gstart,my_gstop
3071 : : !
3072 : : ! Compute the index of G-G2 wave in the FFT grid.
3073 [ + + ]: 1609860 : do ii=1,npw
3074 [ + + ]: 6347217 : g1mg2(:,ii) = Gsph%gvec(:,ii) - Gsph%gvec(:,ig2)
3075 : : end do
3076 : 30741 : call kgindex(igfft,g1mg2,mask,MPI_enreg_seq,ngfft,npw)
3077 : :
3078 : : ! TODO can use zero-padding FFT to speed up the transform.
3079 : : !call sphereboundary(gbound,istwfk1,g1mg2,mgfft,npw)
3080 : :
3081 : : ! Evaluate em1_qpg2r = \int em1(|q+G2|,r) e^{-i(G1-G2).r} dr }.
3082 [ + + ]: 122964 : qpg2 = qpt_bz + Gsph%gvec(:,ig2)
3083 : 30741 : qpg2_nrm = normv(qpg2,Cryst%gmet,"G")
3084 : :
3085 [ + + ]: 62113 : do iw=1,nomega
3086 : : !
3087 : 30741 : select case (model_type)
3088 : : case (1)
3089 [ + + ]: 140247573 : do ifft=1,nfft
3090 : 140247573 : em1_qpg2r(ifft) = one / mdielf_bechstedt(eps_inf,qpg2_nrm,rhor(ifft,1))
3091 : : end do
3092 : : case default
3093 [ + - ]: 30741 : ABI_ERROR(sjoin("Unknown model_type:",itoa(model_type)))
3094 : : end select
3095 : :
3096 : 30741 : call fourdp(cplex1,fofg,em1_qpg2r,-1,MPI_enreg_seq,nfft,1,ngfft,tim_fourdp0)
3097 : : !
3098 : : ! Here, unlike the other parts of the code, the unsymmetrized e^{-1} is used.
3099 [ + + ]: 1640601 : do ig1=1,npw
3100 : 1579119 : g1mg2_fft = igfft(ig1)
3101 : 1609860 : w_qbz(ig1,ig2,iw) = DCMPLX(fofg(1,g1mg2_fft), fofg(2,g1mg2_fft)) * vc_qbz(ig2) !/ Cryst%ucvol
3102 : : end do
3103 : : end do ! iw
3104 : : end do ! ig2
3105 : :
3106 [ - + ]: 631 : ABI_FREE(em1_qpg2r)
3107 [ - + ]: 631 : ABI_FREE(fofg)
3108 [ - + ]: 631 : ABI_FREE(igfft)
3109 [ - + ]: 631 : ABI_FREE(g1mg2)
3110 [ - + ]: 631 : ABI_FREE(mask)
3111 : : !
3112 : : ! W = 1/2 * (A + A^H)
3113 : : ! The MPI sum is done inside the loop to avoid problems with the size of the packet.
3114 [ - + + - : 2524 : ABI_MALLOC_OR_DIE(ctmp,(npw,npw), ierr)
- + + - +
- + - + -
- + - + ]
3115 : :
3116 [ + + ]: 1262 : do iw=1,nomega
3117 : : !ctmp = TRANSPOSE(CONJG(w_qbz(:,:,iw)))
3118 [ - + + + : 1724402 : ctmp = GWPC_CONJG(w_qbz(:,:,iw))
+ + ]
3119 : 631 : call sqmat_itranspose(npw,ctmp)
3120 [ + + + + ]: 1723771 : w_qbz(:,:,iw) = half * (ctmp + w_qbz(:,:,iw))
3121 : 1262 : call xmpi_sum(w_qbz(:,:,iw),comm,ierr)
3122 : : end do
3123 : : !
3124 : : ! Calculate the symmetrized Em1. W = vc(G1)^{1/2} \tilde Em1 vc(G2)^{1/2} -------------------------
3125 [ + - ]: 631 : if (toupper(which)=="EM1") then
3126 [ + + ]: 33260 : do ig=1,npw
3127 : 32629 : isg = Gsph%rottb(ig,itim_q,isym_q)
3128 : 33260 : vc_qbz(isg) = vc_sqrt_ibz(ig) ! Workspace storing vc*{1/2}(q_BZ,G).
3129 : : end do
3130 : :
3131 [ + + ]: 33260 : do ig2=1,npw
3132 [ + + ]: 1723771 : do ig1=1,npw
3133 : 1723140 : ctmp(ig1,ig2) = one / (vc_qbz(ig1) * vc_qbz(ig2))
3134 : : end do
3135 : : end do
3136 : :
3137 [ + + ]: 1893 : do iw=1,nomega
3138 [ + + + + ]: 1724402 : w_qbz(:,:,iw) = w_qbz(:,:,iw) * ctmp(:,:)
3139 : : end do
3140 : : end if
3141 : :
3142 : 631 : call destroy_mpi_enreg(MPI_enreg_seq)
3143 : :
3144 [ - + ]: 631 : ABI_FREE(vc_qbz)
3145 [ + - ]: 631 : ABI_SFREE(ctmp)
3146 : :
3147 [ - + - + : 1893 : end subroutine screen_mdielf
- + - + -
+ - + - +
- + - + -
+ - + -
+ ]
3148 : : !!***
3149 : :
3150 : : !----------------------------------------------------------------------
3151 : :
3152 : : !!****f* m_screening/chi_new
3153 : : !! NAME
3154 : : !! chi_new
3155 : : !!
3156 : : !! FUNCTION
3157 : : !!
3158 : : !! INPUTS
3159 : : !!
3160 : : !! OUTPUT
3161 : : !!
3162 : : !! SOURCE
3163 : :
3164 : 0 : type(chi_t) function chi_new(npwe, nomega) result(chi)
3165 : :
3166 : : !Arguments ------------------------------------
3167 : : integer,intent(in) :: npwe,nomega
3168 : : ! *************************************************************************
3169 : :
3170 : 0 : chi%nomega = nomega; chi%npwe = npwe
3171 : :
3172 [ # # # # : 0 : ABI_MALLOC(chi%mat, (npwe,npwe,nomega))
# # # # #
# # # # #
# # # # ]
3173 : :
3174 [ # # # # : 0 : ABI_MALLOC(chi%head, (3,3,nomega))
# # # # #
# # # ]
3175 [ # # # # : 0 : ABI_MALLOC(chi%lwing, (npwe, nomega,3))
# # # # #
# # # # #
# # # # ]
3176 [ # # # # : 0 : ABI_MALLOC(chi%uwing, (npwe, nomega,3))
# # # # #
# # # ]
3177 : :
3178 [ # # ]: 0 : end function chi_new
3179 : : !!***
3180 : :
3181 : : !----------------------------------------------------------------------
3182 : :
3183 : : !!****f* m_screening/chi_free
3184 : : !! NAME
3185 : : !! chi_free
3186 : : !!
3187 : : !! FUNCTION
3188 : : !!
3189 : : !! INPUTS
3190 : : !!
3191 : : !! OUTPUT
3192 : : !!
3193 : : !! SOURCE
3194 : :
3195 : 0 : subroutine chi_free(chi)
3196 : :
3197 : : !Arguments ------------------------------------
3198 : : !scalars
3199 : : class(chi_t),intent(inout) :: chi
3200 : : ! *************************************************************************
3201 : :
3202 [ # # ]: 0 : ABI_SFREE(chi%mat)
3203 [ # # ]: 0 : ABI_SFREE(chi%head)
3204 [ # # ]: 0 : ABI_SFREE(chi%lwing)
3205 [ # # ]: 0 : ABI_SFREE(chi%uwing)
3206 : :
3207 : 0 : end subroutine chi_free
3208 : : !!***
3209 : :
3210 : : !----------------------------------------------------------------------
3211 : :
3212 : : !!****f* m_screening/lwl_write
3213 : : !! NAME
3214 : : !! lwl_write
3215 : : !!
3216 : : !! FUNCTION
3217 : : !!
3218 : : !! INPUTS
3219 : : !!
3220 : : !! OUTPUT
3221 : : !!
3222 : : !! SOURCE
3223 : :
3224 : 0 : subroutine lwl_write(path, cryst, vcp, npwe, nomega, gvec, chi0, chi0_head, chi0_lwing, chi0_uwing, comm)
3225 : :
3226 : : !Arguments ------------------------------------
3227 : : !scalars
3228 : : integer,intent(in) :: npwe,nomega,comm
3229 : : character(len=*),intent(in) :: path
3230 : : type(crystal_t),intent(in) :: cryst
3231 : : type(vcoul_t),intent(in) :: Vcp
3232 : : !arrays
3233 : : integer,intent(in) :: gvec(3,npwe)
3234 : : complex(gwp),intent(in) :: chi0(npwe,npwe,nomega)
3235 : : complex(dp),intent(inout) :: chi0_head(3,3,nomega),chi0_lwing(npwe,nomega,3),chi0_uwing(npwe,nomega,3)
3236 : :
3237 : : !Local variables-------------------------------
3238 : : !scalars
3239 : : integer,parameter :: master=0,prtvol=0
3240 : : integer :: iw,ii,iomode,unt,my_rank
3241 : : character(len=500) :: msg
3242 : : real(dp) :: length
3243 : : complex(dp) :: wng(3),em1_00
3244 : : ! type(hscr_t),intent(out) :: hscr
3245 : : !arrays
3246 : 0 : complex(dp),allocatable :: wtest(:),eps_head(:,:,:)
3247 : : ! *************************************************************************
3248 : :
3249 : : !if (xmpi_comm_rank(comm) /= master) goto 100
3250 : 0 : my_rank = xmpi_comm_rank(comm)
3251 : :
3252 [ # # # # : 0 : ABI_MALLOC(wtest,(npwe))
# # # # #
# # # ]
3253 : :
3254 : : if (prtvol > 0 .and. my_rank == master) then
3255 : : call wrtout(std_out, "head of chi0")
3256 : : do iw=1,nomega
3257 : : call print_arr([std_out], chi0_head(:,:,iw), max_r=3, max_c=3)
3258 : : end do
3259 : :
3260 : : do iw=1,nomega
3261 : : call wrtout(std_out, "symmetrized e_00 via tensor")
3262 : : wng = MATMUL(chi0_head(:,:,iw),GW_Q0_DEFAULT)
3263 : : write(std_out,*) one - vdotw(GW_Q0_DEFAULT,wng,cryst%gmet,"G") * Vcp%vcqlwl_sqrt(1,1)*Vcp%vcqlwl_sqrt(1,1)
3264 : :
3265 : : call wrtout(std_out, "symmetrized e_0G via tensor")
3266 : : do ii=1,npwe
3267 : : wng = chi0_uwing(ii,iw,:)
3268 : : wtest(ii) = - vdotw(GW_Q0_DEFAULT,wng,cryst%gmet,"G") * Vcp%vcqlwl_sqrt(1,1) * Vcp%vcqlwl_sqrt(ii,1)
3269 : : end do
3270 : : call print_arr([std_out], wtest,max_r=9)
3271 : :
3272 : : call wrtout(std_out, "symmetrized e_G0 via tensor")
3273 : : do ii=1,npwe
3274 : : wng = chi0_lwing(ii,iw,:)
3275 : : wtest(ii) = - vdotw(GW_Q0_DEFAULT,wng,cryst%gmet,"G") * Vcp%vcqlwl_sqrt(1,1) * Vcp%vcqlwl_sqrt(ii,1)
3276 : : end do
3277 : : call print_arr([std_out], wtest,max_r=9)
3278 : : end do
3279 : : end if
3280 : :
3281 : : ! Write chi0 data
3282 [ # # ]: 0 : iomode = IO_MODE_FORTRAN; if (endswith(path, ".nc")) iomode = IO_MODE_ETSF
3283 : :
3284 [ # # ]: 0 : if (my_rank == master) then
3285 [ # # ]: 0 : if (iomode == IO_MODE_FORTRAN) then
3286 [ # # ]: 0 : if (open_file(path,msg,newunit=unt,form="unformatted", action="write") /= 0) then
3287 : 0 : ABI_ERROR(msg)
3288 : : end if
3289 : : !call epsm1%hscr%io(fform,rdwr,unt,comm,master,iomode)
3290 [ # # ]: 0 : do iw=1,nomega
3291 : 0 : write(unt)chi0_head(:,:,iw)
3292 : : end do
3293 : : !do iw=1,nomega
3294 : : ! write(unt)chi0_lwing(:,iw,:)
3295 : : !end do
3296 : : !do iw=1,nomega
3297 : : ! write(unt)chi0_uwing(:,iw,:)
3298 : : !end do
3299 : :
3300 : : else
3301 : 0 : ABI_ERROR(sjoin("iomode", itoa(iomode), "is not supported"))
3302 : : end if
3303 : : end if
3304 : :
3305 [ # # # # : 0 : ABI_MALLOC(eps_head,(3,3,nomega))
# # # # #
# # # ]
3306 : 0 : call mkem1_q0(npwe,1,1,nomega,cryst,Vcp,gvec,chi0_head,chi0_lwing,chi0_uwing,chi0,eps_head,comm)
3307 : :
3308 [ # # ]: 0 : if (my_rank == master) then
3309 [ # # ]: 0 : if (iomode == IO_MODE_FORTRAN) then
3310 [ # # ]: 0 : do iw=1,nomega
3311 : 0 : write(unt)eps_head(:,:,iw)
3312 : : end do
3313 : : !do iw=1,nomega
3314 : : ! write(unt)chi0_lwing(:,iw,:)
3315 : : !end do
3316 : : !do iw=1,nomega
3317 : : ! write(unt)chi0_uwing(:,iw,:)
3318 : : !end do
3319 : : else
3320 : 0 : ABI_ERROR(sjoin("iomode:", itoa(iomode), "is not supported"))
3321 : : end if
3322 : : end if
3323 : :
3324 [ # # ]: 0 : ABI_FREE(eps_head)
3325 : :
3326 : : if (prtvol > 0 .and. my_rank == master) then
3327 : : length = normv(GW_Q0_DEFAULT,cryst%gmet,"G")
3328 : :
3329 : : do iw=1,nomega
3330 : : em1_00 = one / vdotw(GW_Q0_DEFAULT/length, MATMUL(chi0_head(:,:,iw),GW_Q0_DEFAULT/length),cryst%gmet,"G")
3331 : : call wrtout(std_out, "e^1_{00} from tensor")
3332 : : write(std_out,*) em1_00
3333 : :
3334 : : call wrtout(std_out, "symmetrized e^-1_0G via tensor")
3335 : : do ii=1,npwe
3336 : : wng = chi0_uwing(ii,iw,:)
3337 : : wtest(ii) = em1_00*vdotw(GW_Q0_DEFAULT/length,wng,cryst%gmet,"G")
3338 : : end do
3339 : : wtest(1) = em1_00
3340 : : call print_arr([std_out], wtest,max_r=9)
3341 : :
3342 : : call wrtout(std_out, "symmetrized e^-1_G0 via tensor")
3343 : : do ii=1,npwe
3344 : : wng = chi0_lwing(ii,iw,:)
3345 : : wtest(ii) = em1_00*vdotw(GW_Q0_DEFAULT/length,wng,cryst%gmet,"G")
3346 : : end do
3347 : : wtest(1) = em1_00
3348 : : call print_arr([std_out], wtest,max_r=9)
3349 : : end do !iw
3350 : : end if
3351 : :
3352 [ # # ]: 0 : ABI_FREE(wtest)
3353 : :
3354 [ # # ]: 0 : if (my_rank == master) then
3355 [ # # ]: 0 : if (iomode == IO_MODE_FORTRAN) then
3356 : 0 : close(unt)
3357 : : else
3358 [ # # ]: 0 : NCF_CHECK(nf90_close(unt))
3359 : : end if
3360 : : end if
3361 : :
3362 : : !100 call xmpi_barrier(comm)
3363 : :
3364 : 0 : end subroutine lwl_write
3365 : : !!***
3366 : :
3367 : : !----------------------------------------------------------------------
3368 : :
3369 : : !!****f* m_screening/lwl_init
3370 : : !! NAME
3371 : : !! lwl_init
3372 : : !!
3373 : : !! FUNCTION
3374 : : !!
3375 : : !! INPUTS
3376 : : !!
3377 : : !! OUTPUT
3378 : : !!
3379 : : !! SOURCE
3380 : :
3381 [ # # # # : 0 : subroutine lwl_init(lwl, path, method, cryst, vcp, npwe, gvec, comm)
# # ]
3382 : :
3383 : : !Arguments ------------------------------------
3384 : : !scalars
3385 : : integer,intent(in) :: comm,method,npwe
3386 : : character(len=*),intent(in) :: path
3387 : : type(crystal_t),intent(in) :: cryst
3388 : : type(vcoul_t),intent(in) :: vcp
3389 : : type(lwl_t),intent(out) :: lwl
3390 : : !arrays
3391 : : integer,intent(in) :: gvec(3,npwe)
3392 : :
3393 : : !Local variables-------------------------------
3394 : : !scalars
3395 : : integer,parameter :: master=0
3396 : : integer :: iomode,my_rank,nproc,unt
3397 : : character(len=500) :: msg
3398 : : ! *************************************************************************
3399 : :
3400 : : ABI_UNUSED((/cryst%natom, gvec(1,1), vcp%ng/))
3401 : :
3402 : 0 : my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
3403 : :
3404 [ # # ]: 0 : lwl%fname = path
3405 : 0 : lwl%method = method
3406 [ # # ]: 0 : ABI_CHECK(any(method == [1,2,3]), sjoin("Wrong method:", itoa(method)))
3407 : :
3408 [ # # ]: 0 : iomode = IO_MODE_FORTRAN; if (endswith(path, ".nc")) iomode = IO_MODE_ETSF
3409 : :
3410 : : ! Only master reads.
3411 [ # # ]: 0 : if (my_rank == master) then
3412 : :
3413 : 0 : select case (iomode)
3414 : : case (IO_MODE_FORTRAN)
3415 [ # # ]: 0 : if (open_file(path, msg, newunit=unt, action="read", form="unformatted", status="old") /= 0) then
3416 : 0 : ABI_ERROR(msg)
3417 : : end if
3418 : :
3419 : 0 : close(unt)
3420 : :
3421 : : case default
3422 [ # # ]: 0 : ABI_ERROR(sjoin("iomode:", itoa(iomode), "is not coded"))
3423 : : end select
3424 : : end if
3425 : :
3426 : : ! Broad cast data
3427 : : if (nproc > 1) then
3428 : : end if
3429 : :
3430 [ # # ]: 0 : end subroutine lwl_init
3431 : : !!***
3432 : :
3433 : : !----------------------------------------------------------------------
3434 : :
3435 : : !!****f* m_screening/lwl_free
3436 : : !! NAME
3437 : : !! lwl_free
3438 : : !!
3439 : : !! FUNCTION
3440 : : !!
3441 : : !! INPUTS
3442 : : !!
3443 : : !! OUTPUT
3444 : : !!
3445 : : !! SOURCE
3446 : :
3447 : 0 : subroutine lwl_free(lwl)
3448 : :
3449 : : !Arguments ------------------------------------
3450 : : !scalars
3451 : : type(lwl_t),intent(inout) :: lwl
3452 : : ! *************************************************************************
3453 : :
3454 [ # # ]: 0 : ABI_SFREE(lwl%head)
3455 [ # # ]: 0 : ABI_SFREE(lwl%lwing)
3456 [ # # ]: 0 : ABI_SFREE(lwl%uwing)
3457 [ # # ]: 0 : ABI_SFREE(lwl%body)
3458 : :
3459 : 0 : end subroutine lwl_free
3460 : : !!***
3461 : :
3462 [ # # # # : 0 : end module m_screening
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # ]
3463 : : !!***
|
