Branch data Line data Source code
1 : : !!****m* ABINIT/m_ppmodel
2 : : !! NAME
3 : : !! m_ppmodel
4 : : !!
5 : : !! FUNCTION
6 : : !! Module containing the definition of the ppmodel_t used to deal with the plasmonpole technique.
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2008-2026 ABINIT group (MG, GMR, VO, LR, RWG, RS)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! SOURCE
15 : :
16 : : #if defined HAVE_CONFIG_H
17 : : #include "config.h"
18 : : #endif
19 : :
20 : : #include "abi_common.h"
21 : :
22 : : module m_ppmodel
23 : :
24 : : use defs_basis
25 : : use m_errors
26 : : use m_abicore
27 : : use m_array
28 : : use m_linalg_interfaces
29 : : use m_distribfft
30 : : use m_yaml
31 : :
32 : : use defs_abitypes, only : MPI_type
33 : : use m_fstrings, only : sjoin, itoa, ktoa
34 : : use m_hide_lapack, only : xhegv
35 : : use m_gwdefs, only : GW_Q0_DEFAULT, czero_gw
36 : : use m_crystal, only : crystal_t
37 : : use m_bz_mesh, only : kmesh_t
38 : : use m_gsphere, only : gsphere_t
39 : : use m_vcoul, only : vcoul_t
40 : : use m_qplusg, only : cmod_qpg
41 : : use m_fft_mesh, only : g2ifft
42 : : use m_fft, only : fourdp
43 : : use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq
44 : : use m_pstat, only : pstat_proc
45 : :
46 : : implicit none
47 : :
48 : : private
49 : : !!***
50 : :
51 : : integer,public,parameter :: PPM_NONE = 0
52 : : integer,public,parameter :: PPM_GODBY_NEEDS = 1
53 : : integer,public,parameter :: PPM_HYBERTSEN_LOUIE = 2
54 : : integer,public,parameter :: PPM_LINDEN_HORSH = 3
55 : : integer,public,parameter :: PPM_ENGEL_FARID = 4
56 : :
57 : : ! Flags giving the status of the plasmon-pole tables
58 : : integer,public,parameter :: PPM_NOTAB = 0
59 : : integer,public,parameter :: PPM_TAB_ALLOCATED = 1
60 : : integer,public,parameter :: PPM_TAB_STORED = 2
61 : :
62 : : !----------------------------------------------------------------------
63 : :
64 : : !!****t* m_ppmodel/ppmodel_t
65 : : !! NAME
66 : : !! ppmodel_t
67 : : !!
68 : : !! FUNCTION
69 : : !! This datatype gathers all the information on the Plasmonpole technique used in the calculations
70 : : !!
71 : : !! SOURCE
72 : :
73 : : type,public :: ppmodel_t
74 : :
75 : : integer :: dm2_botsq
76 : : ! =npwc if ppmodel=1,2
77 : : ! =1 if ppmodel=3,4
78 : :
79 : : integer :: dm_eig
80 : : ! =0 if ppmodel=1,2,4
81 : : ! =npwc if ppmodel=3
82 : :
83 : : integer :: dm2_otq
84 : : ! =npwc if ppmodel=1,2
85 : : ! =1 if ppmodel=3,4
86 : :
87 : : integer :: invalid_freq
88 : : ! what to do when PPM frequencies are invalid.
89 : :
90 : : integer :: model
91 : : ! The type of Plasmonpole model.
92 : :
93 : : integer :: mqmem
94 : : ! =nqibz if in-core solution.
95 : : ! =0 for out-of-core for which the last dimension in the ppm arrays has size 1.
96 : :
97 : : integer :: nqibz
98 : : ! Number of q-points in the IBZ
99 : :
100 : : integer :: npwc
101 : : ! Number of G vectors in $\tilde \epsilon $
102 : :
103 : : integer :: userho
104 : : ! 1 if the ppmodel requires rho(G).
105 : :
106 : : integer :: iq_bz = 0
107 : : ! The index of the q-point in the BZ that is referenced by the internal pointer.
108 : : ! when we perform a symmetrization from the IBZ to the BZ.
109 : :
110 : : real(dp) :: drude_plsmf
111 : : ! Drude plasma frequency
112 : :
113 : : real(dp) :: force_plsmf = zero
114 : : ! Force plasma frequency to that set by ppmfreq (not used if set zero).
115 : :
116 : : ! arrays
117 : : logical,allocatable :: keep_qibz(:)
118 : : ! (nqibz)
119 : : ! .TRUE. if the ppmodel tables for this q in the IBZ are kept in memory.
120 : :
121 : : integer,allocatable :: has_qibz(:)
122 : : ! Flag defining the status of the tables for the different q. See the PPM_TAB flags.
123 : :
124 : : complex(gwp),allocatable :: bigomegatwsq_qbz_vals(:,:)
125 : : ! (Points|Stores) the symmetrized plasmon pole parameters $\tilde\Omega^2_{G Gp}(q_bz)$.
126 : :
127 : : complex(gwp),allocatable :: omegatw_qbz_vals(:,:)
128 : : ! (Points|Stores) the symmetrized plasmon pole parameters $\tilde\omega_{G Gp}(q_bz)$.
129 : :
130 : : complex(gwp),allocatable :: eigpot_qbz_vals(:,:)
131 : : ! (Points|Stores) the eigvectors of the symmetrized inverse dielectric matrix.
132 : :
133 : : type(array2_gwpc_t),allocatable :: bigomegatwsq(:)
134 : : ! (nqibz)%value(npwc,dm2_botsq)
135 : : ! Plasmon pole parameters $\tilde\Omega^2_{G Gp}(q)$.
136 : :
137 : : type(array2_gwpc_t),allocatable :: omegatw(:)
138 : : ! (nqibz)%value(npwc,dm2_otq)
139 : : ! Plasmon pole parameters $\tilde\omega_{G Gp}(q)$.
140 : :
141 : : type(array2_gwpc_t),allocatable :: eigpot(:)
142 : : ! (nqibz)%value(dm_eig,dm_eig)
143 : : ! Eigvectors of the symmetrized inverse dielectric matrix.
144 : :
145 : : contains
146 : :
147 : : procedure :: get_qbz => ppm_get_qbz
148 : : ! Symmetrize the ppm parameters in the BZ.
149 : :
150 : : procedure :: init => ppm_init
151 : : ! Initialize dimensions and pointers
152 : :
153 : : procedure :: free => ppm_free
154 : : ! Free dynamic memory.
155 : :
156 : : procedure :: setup => ppm_setup
157 : : ! Main Driver
158 : :
159 : : procedure :: new_setup => ppm_new_setup
160 : : ! New Main Driver
161 : :
162 : : procedure :: print => ppm_print
163 : : ! Print info on object
164 : :
165 : : procedure :: calc_sigc => ppm_calc_sigc
166 : : ! Matrix elements of the correlated self-energy with ppmodel.
167 : :
168 : : procedure :: calc_dsigc => ppm_calc_dsigc
169 : : ! Frequency derivative of the above (double-pole term for GWPT at q=Gamma).
170 : :
171 : : procedure :: rotate_iqbz => ppm_rotate_iqbz
172 : :
173 : : procedure :: malloc_iqibz => ppm_malloc_iqibz
174 : :
175 : : procedure :: table_free_iqibz => ppm_table_free_iqibz
176 : :
177 : : procedure :: get_eigenvalues => ppm_get_eigenvalues
178 : :
179 : : procedure :: getem1 => ppm_getem1
180 : : ! Reconstruct e^{-1}(w) from ppm.
181 : :
182 : : procedure :: getem1_one_ggp => ppm_getem1_one_ggp
183 : : ! Reconstruct e^{-1}(w) from ppm for one (G,G') pair
184 : :
185 : : end type ppmodel_t
186 : :
187 : : public :: cqratio
188 : : !!***
189 : :
190 : : contains
191 : : !!***
192 : :
193 : : !!****f* m_ppmodel/ppm_get_qbz
194 : : !! NAME
195 : : !! ppm_get_qbz
196 : : !!
197 : : !! FUNCTION
198 : : !! Compute plasmon-pole matrix elements for q in the BZ from the symmetrical image in the IBZ
199 : : !!
200 : : !! INPUTS
201 : : !! Gsph<gsphere_t>=data related to the G-sphere
202 : : !! Qmesh<kmesh_t>=Info on the q-mesh
203 : : !! iq_bz=Index of the q-point in the BZ where ppmodel parameters have to be symmetrized
204 : : !!
205 : : !! OUTPUT
206 : : !! botsq
207 : : !! otq
208 : : !! eig (only if ppm%ppmodel==3)
209 : : !!
210 : : !! NOTES
211 : : !! In the present implementation we are not considering a possible umklapp vector G0.
212 : : !! In this case,indeed, the equation is different since we have to consider G-G0.
213 : : !! There is however a check in sigma
214 : : !!
215 : : !! * Remember the symmetry properties of \tilde\espilon^{-1}
216 : : !!
217 : : !! If q_bz=Sq_ibz+G0:
218 : : !!
219 : : !! $\epsilon^{-1}_{SG1-G0,SG2-G0}(q_bz) = e^{+iS(G2-G1).\tau}\epsilon^{-1}_{G1,G2)}(q)
220 : : !!
221 : : !! If time-reversal symmetry can be used then:
222 : : !!
223 : : !! $\epsilon^{-1}_{G1,G2}(-q_bz) = e^{+i(G1-G2).\tau}\epsilon^{-1}_{-S^{-1}(G1+Go),-S^{-1}(G2+G0)}^*(q)
224 : : !!
225 : : !! * Note that eig is used only if ppm%model==3
226 : : !!
227 : : !! SOURCE
228 : :
229 : 14848 : subroutine ppm_get_qbz(ppm, Gsph, Qmesh, iq_bz, botsq, otq, eig)
230 : :
231 : : !Arguments ------------------------------------
232 : : !scalars
233 : : class(ppmodel_t),target,intent(inout) :: ppm
234 : : integer,intent(in) :: iq_bz
235 : : type(gsphere_t),target,intent(in) :: Gsph
236 : : type(kmesh_t),intent(in) :: Qmesh
237 : : complex(gwp),allocatable,intent(out) :: botsq(:,:),otq(:,:),eig(:,:)
238 : :
239 : : !Local variables-------------------------------
240 : : !scalars
241 : : integer :: ii,jj,iq_ibz,itim_q,isym_q,iq_curr,isg1,isg2
242 : : !arrays
243 : 14848 : integer, contiguous, pointer :: grottb(:)
244 : 14848 : complex(gwp),contiguous, pointer :: phsgt(:),bigomegatwsq(:,:),omegatw(:,:)
245 : : ! *********************************************************************
246 : :
247 : : ! Save the index of the q-point for checking purpose.
248 : 14848 : ppm%iq_bz = iq_bz
249 : :
250 [ + - - + : 59392 : ABI_MALLOC(botsq, (ppm%npwc, ppm%dm2_botsq))
+ - + - +
- - + - +
- + ]
251 [ - + + - : 59392 : ABI_MALLOC(otq, (ppm%npwc, ppm%dm2_otq))
- + + - +
- + - - +
- + - + ]
252 [ - + + - : 59392 : ABI_MALLOC(eig, (ppm%dm_eig, ppm%dm_eig))
- + + - +
- + - - +
- + - + ]
253 : :
254 : : ! Here there is a problem with the small q, still cannot use BZ methods
255 : 14848 : iq_ibz = Qmesh%tab(iq_bz); isym_q = Qmesh%tabo(iq_bz); itim_q = (3-Qmesh%tabi(iq_bz))/2
256 : :
257 [ + + + - : 59392 : ABI_CHECK(all(abs(qmesh%umklp(:, iq_bz)) < tol6), "umklapp in get_qbz are not supported!")
- + ]
258 : :
259 : : !call Qmesh%get_bz_item(iq_bz,qbz,iq_ibz,isym_q,itim_q,isirred=q_isirred)
260 [ + + ]: 14848 : iq_curr = iq_ibz; if (ppm%mqmem == 0) iq_curr = 1
261 : :
262 : 14848 : grottb => Gsph%rottb (1:ppm%npwc, itim_q, isym_q)
263 : 14848 : phsgt => Gsph%phmSGt(1:ppm%npwc, isym_q)
264 : 14848 : bigomegatwsq => ppm%bigomegatwsq(iq_curr)%vals
265 : 14848 : omegatw => ppm%omegatw(iq_curr)%vals
266 : :
267 : : ! Symmetrize the PPM parameters
268 : 14848 : select case (ppm%model)
269 : : case (PPM_GODBY_NEEDS, PPM_HYBERTSEN_LOUIE)
270 : : ! Plasmon pole frequencies otq are invariant under symmetry
271 : : !$omp parallel do private(isg1, isg2)
272 [ + + ]: 474726 : do jj=1,ppm%npwc
273 : 460070 : isg2 = grottb(jj)
274 [ + + ]: 21095446 : do ii=1,ppm%npwc
275 : 20620720 : isg1 = grottb(ii)
276 : 20620720 : botsq(isg1,isg2) = bigomegatwsq(ii,jj)*phsgt(ii)*CONJG(phsgt(jj))
277 : 21080790 : otq (isg1,isg2) = omegatw(ii,jj)
278 : : end do
279 : : end do
280 : :
281 : : case (PPM_LINDEN_HORSH)
282 : : ! For notations see page 22 of Quasiparticle Calculations in solid (Aulbur et al)
283 : : ! If q_bz = Sq_ibz + G0 then:
284 : : !
285 : : ! $\omega^2_{ii}(q_bz) = \omega^2_{ii}(q)$ (otq array)
286 : : ! $\alpha_{ii}(q_bz) = \alpha_{ii}(q)$ (botq array
287 : : ! $\Phi_{SG-G0}(q_bz) = \Phi_{G}(q) e^{-iSG.t}$ (eigenvectors of e^{-1}, eig array)
288 : : !
289 [ + + ]: 2688 : do ii=1,ppm%npwc ! DM bands index
290 : 2592 : botsq(ii,1) = bigomegatwsq(ii,1)
291 : 2592 : otq (ii,1) = omegatw (ii,1)
292 [ + + ]: 72672 : do jj=1,ppm%npwc
293 : 72576 : eig(grottb(jj),ii) = ppm%eigpot(iq_curr)%vals(jj,ii) * phsgt(jj)
294 : : end do
295 : : end do
296 [ - + - - : 96 : if (itim_q==2) eig=CONJG(eig) ! Time-reversal
- - ]
297 : :
298 : : case (PPM_ENGEL_FARID)
299 : : ! For notations see page 23 of Quasiparticle Calculations in solid (Aulbur et al.)
300 : : ! If q_bz = Sq_ibz + G0 then:
301 : : !
302 : : ! $\omega^2_{ii}(q_bz) = \omega^2_{ii}(q)$ (otq array)
303 : : ! $y_{SG-G0}(q_bz) = y_{G}(q) e^{-iSG.t}$ (y=Lx)
304 : : !
305 [ + + ]: 2688 : do ii=1,ppm%npwc ! DM bands index
306 : 2592 : otq(ii,1) = omegatw(ii,1)
307 [ + + ]: 72672 : do jj=1,ppm%npwc
308 : 72576 : botsq(grottb(jj),ii) = bigomegatwsq(jj,ii)*phsgt(jj)
309 : : end do
310 : : end do
311 : :
312 : : case default
313 [ + + + - ]: 14848 : ABI_BUG(sjoin('Wrong ppm%model:',itoa(ppm%model)))
314 : : end select
315 : :
316 : : ! Take into account time-reversal symmetry.
317 [ + + ]: 14752 : if (itim_q == 2) then
318 : : !$omp parallel workshare
319 [ + + + + ]: 809858 : botsq=CONJG(botsq)
320 : : !$omp end parallel workshare
321 : : end if
322 : :
323 [ - + ]: 14848 : end subroutine ppm_get_qbz
324 : : !!***
325 : :
326 : : !----------------------------------------------------------------------
327 : :
328 : : !!****f* m_ppmodel/ppm_free
329 : : !! NAME
330 : : !! ppm_free
331 : : !!
332 : : !! FUNCTION
333 : : !! Deallocate all associated pointers defined in a variable of type ppmodel_t.
334 : : !!
335 : : !! SOURCE
336 : :
337 : 238 : subroutine ppm_free(ppm)
338 : :
339 : : !Arguments ------------------------------------
340 : : class(ppmodel_t),intent(inout) :: ppm
341 : :
342 : : !Local variables-------------------------------
343 : : integer :: dim_q,iq_ibz
344 : : ! *********************************************************************
345 : :
346 [ - + ]: 238 : ABI_SFREE(ppm%bigomegatwsq_qbz_vals)
347 [ - + ]: 238 : ABI_SFREE(ppm%omegatw_qbz_vals)
348 [ - + ]: 238 : ABI_SFREE(ppm%eigpot_qbz_vals)
349 : :
350 [ + + ]: 238 : dim_q = ppm%nqibz; if (ppm%mqmem==0) dim_q=1
351 : :
352 [ + + ]: 238 : if (allocated(ppm%bigomegatwsq)) then
353 [ + + ]: 1010 : do iq_ibz=1,dim_q
354 : 1010 : call ppm%bigomegatwsq(iq_ibz)%free()
355 : : end do
356 [ + - + + : 1010 : ABI_FREE(ppm%bigomegatwsq)
- + - + ]
357 : : end if
358 [ + + ]: 238 : if (allocated(ppm%omegatw)) then
359 [ + + ]: 1010 : do iq_ibz=1,dim_q
360 : 1010 : call ppm%omegatw(iq_ibz)%free()
361 : : end do
362 [ + - + + : 1010 : ABI_FREE(ppm%omegatw)
- + - + ]
363 : : end if
364 [ + + ]: 238 : if (allocated(ppm%eigpot)) then
365 [ + + ]: 1010 : do iq_ibz=1,dim_q
366 : 1010 : call ppm%eigpot(iq_ibz)%free()
367 : : end do
368 [ + - + + : 1010 : ABI_FREE(ppm%eigpot)
- + - + ]
369 : : end if
370 : :
371 : : ! logical flags must be deallocated here.
372 [ + + ]: 238 : ABI_SFREE(ppm%keep_qibz)
373 [ + + ]: 238 : ABI_SFREE(ppm%has_qibz)
374 : :
375 : 238 : end subroutine ppm_free
376 : : !!***
377 : :
378 : : !----------------------------------------------------------------------
379 : :
380 : : !!****f* m_ppmodel/ppm_malloc_iqibz
381 : : !! NAME
382 : : !! ppm_malloc_iqibz
383 : : !!
384 : : !! FUNCTION
385 : : !! Allocate the ppmodel tables for the selected q-point in the IBZ.
386 : : !!
387 : : !! INPUT
388 : : !! iq_ibz=Index of the q-point in the IBZ.
389 : : !!
390 : : !! SOURCE
391 : :
392 : 872 : subroutine ppm_malloc_iqibz(ppm, iq_ibz)
393 : :
394 : : !Arguments ------------------------------------
395 : : class(ppmodel_t),intent(inout) :: ppm
396 : : integer,intent(in) :: iq_ibz
397 : :
398 : : !Local variables-------------------------------
399 : : integer :: ierr
400 : : ! *********************************************************************
401 : :
402 [ - + ]: 872 : ABI_CHECK(allocated(ppm%bigomegatwsq), "bigomegatwsq is not allocated")
403 [ - + ]: 872 : ABI_CHECK(allocated(ppm%omegatw), "omegatwsq is not allocated")
404 [ - + ]: 872 : ABI_CHECK(allocated(ppm%eigpot), "eigpot is not allocated")
405 : :
406 [ - + ]: 872 : ABI_CHECK_IGEQ(size(ppm%bigomegatwsq), iq_ibz, "bigomegatwsq too small")
407 [ - + ]: 872 : ABI_CHECK_IGEQ(size(ppm%omegatw), iq_ibz, "omegatwsq too small")
408 [ - + ]: 872 : ABI_CHECK_IGEQ(size(ppm%eigpot), iq_ibz, "eigpot too small")
409 : :
410 [ - + + - : 3488 : ABI_MALLOC_OR_DIE(ppm%bigomegatwsq(iq_ibz)%vals, (ppm%npwc, ppm%dm2_botsq), ierr)
- + + - +
- + - + -
+ - - + +
- - + ]
411 [ - + + - : 3488 : ABI_MALLOC_OR_DIE(ppm%omegatw(iq_ibz)%vals, (ppm%npwc, ppm%dm2_otq), ierr)
- + + - +
- + - + -
+ - - + +
- - + ]
412 [ - + + - : 3488 : ABI_MALLOC_OR_DIE(ppm%eigpot(iq_ibz)%vals, (ppm%dm_eig, ppm%dm_eig), ierr)
- + + - +
- + - + -
+ - - + +
- - + ]
413 : :
414 : 872 : ppm%has_qibz(iq_ibz) = PPM_TAB_ALLOCATED
415 : :
416 : 872 : end subroutine ppm_malloc_iqibz
417 : : !!***
418 : :
419 : : !----------------------------------------------------------------------
420 : :
421 : : !!****f* m_ppmodel/ppm_table_free_iqibz
422 : : !! NAME
423 : : !! ppm_table_free_iqibz
424 : : !!
425 : : !! FUNCTION
426 : : !! Free the ppmodel tables for the selected q-point in the IBZ.
427 : : !!
428 : : !! INPUT
429 : : !! iq_ibz = Index of the q-point in the IBZ.
430 : : !!
431 : : !! SOURCE
432 : :
433 : 0 : subroutine ppm_table_free_iqibz(ppm, iq_ibz)
434 : :
435 : : !Arguments ------------------------------------
436 : : class(ppmodel_t),intent(inout) :: ppm
437 : : integer,intent(in) :: iq_ibz
438 : : ! *********************************************************************
439 : :
440 [ # # ]: 0 : if (allocated(ppm%bigomegatwsq)) call ppm%bigomegatwsq(iq_ibz)%free()
441 [ # # ]: 0 : if (allocated(ppm%omegatw)) call ppm%omegatw(iq_ibz)%free()
442 [ # # ]: 0 : if (allocated(ppm%eigpot)) call ppm%eigpot(iq_ibz)%free()
443 : :
444 : 0 : ppm%has_qibz(iq_ibz) = PPM_NOTAB
445 : :
446 : 0 : end subroutine ppm_table_free_iqibz
447 : : !!***
448 : :
449 : : !----------------------------------------------------------------------
450 : :
451 : : !!****f* m_ppmodel/ppm_init
452 : : !! NAME
453 : : !! ppm_init
454 : : !!
455 : : !! FUNCTION
456 : : !! Initialize dimensions and other important variables related to the ppmodel
457 : : !!
458 : : !! INPUTS
459 : : !! ppmodel=
460 : : !! drude_plsmf=
461 : : !!
462 : : !! SOURCE
463 : :
464 [ + - ]: 138 : subroutine ppm_init(ppm, mqmem, nqibz, npwe, ppmodel, drude_plsmf, invalid_freq)
465 : :
466 : : !Arguments ------------------------------------
467 : : class(ppmodel_t),intent(out) :: ppm
468 : : integer,intent(in) :: mqmem, nqibz, npwe, ppmodel, invalid_freq
469 : : real(dp),intent(in) :: drude_plsmf
470 : :
471 : : !Local variables-------------------------------
472 : : !scalars
473 : : integer :: dim_q,iq_ibz
474 : : logical :: ltest
475 : : !character(len=500) :: msg
476 : : ! *********************************************************************
477 : :
478 : 138 : ppm%nqibz = nqibz; ppm%mqmem = mqmem; ppm%invalid_freq = invalid_freq
479 : : !call wrtout(std_out, sjoin(' ppm%mqmem:', itoa(ppm%mqmem), 'ppm%nqibz:', itoa(ppm%nqibz)))
480 : : ltest = (ppm%mqmem == 0 .or. ppm%mqmem == ppm%nqibz)
481 : : ! ABI_CHECK(ltest,'Wrong value for mqmem')
482 : :
483 : 138 : ppm%npwc = npwe
484 : 138 : ppm%model = ppmodel
485 : 138 : ppm%drude_plsmf = drude_plsmf
486 : 138 : ppm%userho = 0
487 [ + + ]: 138 : if (any(ppmodel == [PPM_HYBERTSEN_LOUIE, PPM_LINDEN_HORSH, PPM_ENGEL_FARID])) ppm%userho = 1
488 : :
489 [ - + + - : 414 : ABI_MALLOC(ppm%keep_qibz, (nqibz))
+ - + - -
+ - + -
+ ]
490 [ + + + + : 1884 : ppm%keep_qibz = .FALSE.; if (ppm%mqmem > 0) ppm%keep_qibz = .TRUE.
+ + ]
491 : :
492 [ - + + - : 276 : ABI_MALLOC(ppm%has_qibz, (nqibz))
- + - + -
+ ]
493 [ + + ]: 1029 : ppm%has_qibz = PPM_NOTAB
494 : :
495 : : ! Full q-mesh is stored or out-of-memory solution.
496 [ + + ]: 138 : dim_q = ppm%nqibz; if (ppm%mqmem == 0) dim_q=1
497 : :
498 [ - + + - : 1286 : ABI_MALLOC(ppm%bigomegatwsq, (dim_q))
+ - + - -
+ - + - +
+ + ]
499 [ - + + - : 1148 : ABI_MALLOC(ppm%omegatw, (dim_q))
- + - + -
+ + + ]
500 [ - + + - : 1148 : ABI_MALLOC(ppm%eigpot, (dim_q))
- + - + -
+ + + ]
501 : :
502 : 138 : select case (ppm%model)
503 : : case (PPM_NONE)
504 : 0 : ABI_WARNING("Called with ppmodel == 0")
505 : 0 : ppm%dm2_botsq = 0
506 : 0 : ppm%dm2_otq = 0
507 : 0 : ppm%dm_eig = 0
508 : 0 : RETURN
509 : :
510 : : case (PPM_GODBY_NEEDS, PPM_HYBERTSEN_LOUIE)
511 : 132 : ppm%dm2_botsq = ppm%npwc
512 : 132 : ppm%dm2_otq = ppm%npwc
513 : 132 : ppm%dm_eig = 1 ! Should be set to 0, but g95 does not like zero-sized arrays
514 : :
515 : : case (PPM_LINDEN_HORSH)
516 : 3 : ppm%dm2_botsq = 1
517 : 3 : ppm%dm2_otq = 1
518 : 3 : ppm%dm_eig = ppm%npwc
519 : :
520 : : case (PPM_ENGEL_FARID)
521 : 3 : ppm%dm2_botsq = ppm%npwc
522 : 3 : ppm%dm2_otq = 1
523 : 3 : ppm%dm_eig = 1 ! Should be set to 0, but g95 does not like zero-sized arrays
524 : :
525 : : case default
526 [ - + + + : 138 : ABI_BUG(sjoin('Wrong ppm%model:', itoa(ppm%model)))
- ]
527 : : end select
528 : :
529 : : ! Allocate tables depending on the value of keep_qibz.
530 [ + + ]: 1010 : do iq_ibz=1,dim_q
531 : 1010 : call ppm%malloc_iqibz(iq_ibz)
532 : : end do
533 : :
534 : 138 : call pstat_proc%print(_PSTAT_ARGS_)
535 : :
536 : : end subroutine ppm_init
537 : : !!***
538 : :
539 : : !----------------------------------------------------------------------
540 : :
541 : : !!****f* m_ppmodel/ppm_setup
542 : : !! NAME
543 : : !! ppm_setup
544 : : !!
545 : : !! FUNCTION
546 : : !! Initialize some values of several arrays of the ppm datastructure
547 : : !! that are used in case of plasmonpole calculations
548 : : !! This is a wrapper around different plasmonpole routines.
549 : : !!
550 : : !! INPUTS
551 : : !! Cryst<crystal_t>=Info on the unit cell and crystal symmetries.
552 : : !! Qmesh<kmesh_t>=the q-mesh used for the inverse dielectric matrix
553 : : !! npwe=number of G vectors for the correlation part
554 : : !! nomega=number of frequencies in $\epsilon^{-1}$
555 : : !! omega=frequencies in epsm1
556 : : !! epsm1=the inverse dielctric matrix
557 : : !! ngfftf(18)=contain all needed information about the 3D fine FFT mesh, see ~abinit/doc/variables/vargs.htm#ngfft
558 : : !! gprimd(3,3)=dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$)
559 : : !! nfftf=the number of points in the FFT mesh (for this processor)
560 : : !! rhor_tot(nfftf)=the total charge in real space
561 : : !!
562 : : !! SIDE EFFECTS
563 : : !! == if ppmodel 1 or 2 ==
564 : : !! %omegatw and %bigomegatwsq
565 : : !! == if ppmodel 3 ==
566 : : !! %omegatw, %bigomegatwsq and %eigpot
567 : : !! == if ppmodel 4 ==
568 : : !! %omegatw and %bigomegatwsq
569 : : !!
570 : : !! NOTES
571 : : !! * FFT parallelism not implemented.
572 : : !! * TODO: rhor_tot should be replaced by rhog_tot
573 : : !!
574 : : !! SOURCE
575 : :
576 [ - + ]: 408 : subroutine ppm_setup(ppm, Cryst, Qmesh, npwe, nomega, omega, epsm1, nfftf, gvec, ngfftf, rhor_tot, &
577 : : iqiA) ! Optional
578 : :
579 : : !Arguments ------------------------------------
580 : : !scalars
581 : : class(ppmodel_t),intent(inout) :: ppm
582 : : integer,intent(in) :: nfftf,npwe,nomega
583 : : integer,intent(in),optional :: iqiA
584 : : type(kmesh_t),intent(in) :: Qmesh
585 : : type(crystal_t),intent(in) :: Cryst
586 : : !arrays
587 : : integer,intent(in) :: gvec(3,npwe),ngfftf(18)
588 : : real(dp),intent(in) :: rhor_tot(nfftf)
589 : : complex(dp),intent(in) :: omega(nomega)
590 : : complex(gwp),intent(in) :: epsm1(:,:,:,:)
591 : :
592 : : !Local variables-------------------------------
593 : : !scalars
594 : : integer :: nqiA,iq_ibz
595 : : real(dp) :: n_at_G_zero
596 : : logical :: single_q
597 : : character(len=500) :: msg
598 : : !scalars
599 : : real(dp) :: qpt(3)
600 : : ! *************************************************************************
601 : :
602 : : !@ppmodel_t
603 : : !
604 : : ! === if iqiA is present, then consider only one qpoint to save memory ===
605 : : ! * This means the object has been already initialized
606 : 408 : nqiA = Qmesh%nibz; single_q = .FALSE.
607 [ + + ]: 408 : if (PRESENT(iqiA)) then
608 : 290 : nqiA = 1; single_q = .TRUE.
609 : : end if
610 : :
611 : : ! Allocate plasmonpole parameters
612 : : ! TODO ppmodel==1 by default, should be set to 0 if AC and CD
613 : 408 : select case (ppm%model)
614 : :
615 : : case (PPM_NONE)
616 : 0 : ABI_COMMENT(' Skipping Plasmompole model calculation')
617 : :
618 : : case (PPM_GODBY_NEEDS)
619 : : ! Note: the q-dependency enters only through epsilon^-1.
620 [ + + ]: 1464 : do iq_ibz=1,nqiA
621 : : call cppm1par(npwe,nomega,omega,ppm%drude_plsmf,&
622 [ - + - - : 1464 : epsm1(:,:,:,iq_ibz),ppm%omegatw(iq_ibz)%vals,ppm%bigomegatwsq(iq_ibz)%vals)
- - - - -
- - + ]
623 : : end do
624 : :
625 : : case (PPM_HYBERTSEN_LOUIE)
626 [ + + ]: 28 : do iq_ibz=1,nqiA
627 [ + + - + : 96 : qpt = Qmesh%ibz(:,iq_ibz); if (single_q) qpt=Qmesh%ibz(:,iqiA)
- - ]
628 : :
629 : : call cppm2par(qpt,npwe,epsm1(:,:,1,iq_ibz),ngfftf,gvec,Cryst%gprimd,rhor_tot,nfftf,Cryst%gmet,&
630 [ - + - - : 28 : ppm%bigomegatwsq(iq_ibz)%vals,ppm%omegatw(iq_ibz)%vals,ppm%invalid_freq)
- - - - -
+ ]
631 : : end do
632 : :
633 : : ! Quick-and-dirty change of the plasma frequency. Never executed in standard runs.
634 [ - + ]: 4 : if (ppm%force_plsmf>tol6) then ! Integrate the real-space density
635 [ # # ]: 0 : n_at_G_zero = SUM(rhor_tot(:))/nfftf
636 : : ! Change the prefactor
637 : 0 : write(msg,'(2(a,es16.8))') 'Forced ppmfreq:',ppm%force_plsmf*Ha_eV,' nelec/ucvol:',n_at_G_zero
638 : 0 : ABI_WARNING(msg)
639 : 0 : ppm%force_plsmf = (ppm%force_plsmf**2)/(four_pi*n_at_G_zero)
640 [ # # ]: 0 : do iq_ibz=1,ppm%nqibz
641 [ # # # # ]: 0 : ppm%bigomegatwsq(iq_ibz)%vals = ppm%force_plsmf * ppm%bigomegatwsq(iq_ibz)%vals
642 [ # # # # ]: 0 : ppm%omegatw(iq_ibz)%vals = ppm%force_plsmf * ppm%omegatw(iq_ibz)%vals
643 : : end do
644 : 0 : write(msg,'(a,es16.8)') 'Plasma frequency forced in HL ppmodel, new prefactor is:',ppm%force_plsmf
645 : 0 : ABI_WARNING(msg)
646 : : end if
647 : :
648 : : case (PPM_LINDEN_HORSH) ! TODO Check better double precision, this routine is in a messy state
649 [ + + ]: 21 : do iq_ibz=1,nqiA
650 [ + + - + : 72 : qpt = Qmesh%ibz(:,iq_ibz); if (single_q) qpt=Qmesh%ibz(:,iqiA)
- - ]
651 : : call cppm3par(qpt,npwe,epsm1(:,:,1,iq_ibz),ngfftf,gvec,Cryst%gprimd,rhor_tot,nfftf,&
652 [ - + - - : 21 : ppm%bigomegatwsq(iq_ibz)%vals,ppm%omegatw(iq_ibz)%vals(:,1),ppm%eigpot(iq_ibz)%vals)
- - - - -
+ ]
653 : : end do
654 : :
655 : : case (PPM_ENGEL_FARID) ! TODO Check better double precision, this routine is in a messy state
656 [ + + ]: 21 : do iq_ibz=1,nqiA
657 [ + + - + : 72 : qpt = Qmesh%ibz(:,iq_ibz); if (single_q) qpt=Qmesh%ibz(:,iqiA)
- - ]
658 [ + + + + : 30 : if ((ALL(ABS(qpt)<1.0e-3))) qpt = GW_Q0_DEFAULT ! FIXME
+ + ]
659 : : call cppm4par(qpt,npwe,epsm1(:,:,1,iq_ibz),ngfftf,gvec,Cryst%gprimd,rhor_tot,nfftf,&
660 [ - + - - : 21 : ppm%bigomegatwsq(iq_ibz)%vals,ppm%omegatw(iq_ibz)%vals(:,1))
- - - - -
+ ]
661 : : end do
662 : :
663 : : case default
664 [ - + + + : 408 : ABI_BUG(sjoin('Wrong ppm%model:',itoa(ppm%model)))
+ - ]
665 : : end select
666 : :
667 : 408 : end subroutine ppm_setup
668 : : !!***
669 : :
670 : : !----------------------------------------------------------------------
671 : :
672 : : !!****f* m_ppmodel/ppm_getem1
673 : : !! NAME
674 : : !! ppm_getem1
675 : : !!
676 : : !! FUNCTION
677 : : !! Calculate the symmetrized inverse dielectric matrix from the parameters of the plasmon-pole model.
678 : : !!
679 : : !! INPUTS
680 : : !!
681 : : !! OUTPUT
682 : : !!
683 : : !! SOURCE
684 : :
685 : 0 : subroutine ppm_getem1(ppm, mpwc, iqibz, zcut, nomega, omega, Vcp, em1q, &
686 : : only_ig1, only_ig2) ! Optional
687 : :
688 : : !Arguments ------------------------------------
689 : : !scalars
690 : : class(ppmodel_t),intent(in) :: ppm
691 : : integer,intent(in) :: mpwc,iqibz,nomega
692 : : type(vcoul_t),intent(in) :: Vcp
693 : : real(dp),intent(in) :: zcut
694 : : integer,optional,intent(in) :: only_ig1,only_ig2
695 : : !arrays
696 : : complex(dp),intent(in) :: omega(nomega)
697 : : complex(dp),intent(out) :: em1q(mpwc,mpwc,nomega)
698 : :
699 : : !Local variables-------------------------------
700 : : !scalars
701 : : integer :: ig1,ig2,io,idm,ig1_min,ig2_min,ig1_max,ig2_max
702 : : real(dp) :: den
703 : : complex(dp) :: qpg1,qpg2,ug1,ug2
704 : : complex(dp) :: delta,em1ggp,otw,zzpq,yg1,yg2,bot1,bot2,chig1g2
705 : : !character(len=500) :: msg
706 : : ! *************************************************************************
707 : :
708 [ # # ]: 0 : ABI_CHECK(ppm%mqmem/=0,'mqmem==0 not implemented')
709 : :
710 : : !TODO zcut should be an entry in ppm
711 : 0 : delta=CMPLX(zero,zcut)
712 : :
713 : : ! To save memory, a particular combination of
714 : : ! ig1 and ig2 can be selected
715 : 0 : ig1_min = 1
716 : 0 : ig2_min = 1
717 : 0 : ig1_max = ppm%npwc
718 : 0 : ig2_max = ppm%npwc
719 [ # # ]: 0 : if (present(only_ig1)) then
720 : 0 : ig1_min = only_ig1
721 : 0 : ig1_max = only_ig1
722 : : end if
723 [ # # ]: 0 : if (present(only_ig2)) then
724 : 0 : ig2_min = only_ig2
725 : 0 : ig2_max = only_ig2
726 : : end if
727 : :
728 : 0 : select case (ppm%model)
729 : : case (PPM_GODBY_NEEDS, PPM_HYBERTSEN_LOUIE)
730 [ # # ]: 0 : do io=1,nomega
731 [ # # ]: 0 : do ig2=ig2_min,ig2_max
732 [ # # ]: 0 : do ig1=ig1_min,ig1_max
733 : : !den = omega(io)**2-REAL(ppm%omegatw(iqibz)%vals(ig1,ig2)**2)
734 : : !if (den**2<zcut**2) den = omega(io)**2-REAL( (ppm%omegatw(iqibz)%vals(ig1,ig2)-delta)**2 )
735 : 0 : den = omega(io)**2 - REAL( (ppm%omegatw(iqibz)%vals(ig1,ig2)-delta)**2 )
736 : 0 : em1ggp = ppm%bigomegatwsq(iqibz)%vals(ig1,ig2)/den
737 [ # # ]: 0 : if (ig1==ig2) em1ggp=em1ggp+one
738 : 0 : em1q(ig1,ig2,io)=em1ggp
739 : : !em1q(ig1,ig2,io)=em1ggp*Vcp%vc_sqrt(ig1,iqibz)*Vcp%vc_sqrt(ig2,iqibz)
740 : : end do
741 : : end do
742 : : !
743 : : end do !io
744 : :
745 : : case (PPM_LINDEN_HORSH)
746 : : !TODO Check coefficients
747 [ # # ]: 0 : do io=1,nomega
748 [ # # ]: 0 : do ig2=ig2_min,ig2_max
749 [ # # ]: 0 : do ig1=ig1_min,ig1_max
750 : : !
751 : : em1ggp=czero
752 [ # # ]: 0 : do idm=1,ppm%npwc
753 : : !den=omega(io)**2-(ppm%omegatw(iqibz)%vals(ig1,ig2)-delta)**2
754 : : !em1w(io)=em1w(io)+eigvec(ig1,idm,iqibz)*conjg(eigvec(ig2,idm,iqibz))*bigomegatwsq(ig1,ig2,iqibz)/den
755 : 0 : ug1 = ppm%eigpot(iqibz)%vals(ig1,idm)
756 : 0 : ug2 = ppm%eigpot(iqibz)%vals(ig2,idm)
757 : 0 : otw = ppm%bigomegatwsq(iqibz)%vals(idm,1)*ppm%omegatw(iqibz)%vals(idm,1)
758 : 0 : zzpq=ppm%bigomegatwsq(iqibz)%vals(idm,1)
759 : 0 : den=half*REAL(zzpq*otw*( one/(omega(io)-otw+delta) - one/(omega(io)+otw-delta) ))
760 : 0 : em1ggp=em1ggp+ug1*CONJG(ug2)*den
761 : : !eigenvalues(idm,io)=one + half*REAL(zzpq*otw*( one/(omega(io)-otw+delta) - one/(omega(io)+otw-delta) ))
762 : : end do
763 [ # # ]: 0 : if (ig2==ig1) em1ggp=em1ggp+one
764 : 0 : em1q(ig1,ig2,io)=em1ggp
765 : : end do !ig1
766 : : end do !ig2
767 : : !
768 : : end do !iomega
769 : :
770 : : case (PPM_ENGEL_FARID)
771 : : ! Make e^-1
772 [ # # ]: 0 : do io=1,nomega
773 [ # # ]: 0 : do ig2=ig2_min,ig2_max
774 : 0 : qpg2=one/Vcp%vc_sqrt(ig2,iqibz)
775 [ # # ]: 0 : do ig1=ig1_min,ig1_max
776 : 0 : qpg1=one/Vcp%vc_sqrt(ig1,iqibz)
777 : :
778 : 0 : chig1g2=czero
779 [ # # ]: 0 : do idm=1,ppm%npwc
780 : 0 : otw =ppm%omegatw(iqibz)%vals(idm,1)
781 : 0 : bot1=ppm%bigomegatwsq(iqibz)%vals(ig1,idm)
782 : 0 : bot2=ppm%bigomegatwsq(iqibz)%vals(ig2,idm)
783 : 0 : yg1=SQRT(otw/four_pi)*qpg1*bot1
784 : 0 : yg2=SQRT(otw/four_pi)*qpg2*bot2
785 : 0 : chig1g2=chig1g2 + yg1*CONJG(yg2)/(omega(io)**2-(otw-delta)**2)
786 : : end do
787 : :
788 : 0 : em1ggp=four_pi*chig1g2/(qpg1*qpg2)
789 [ # # ]: 0 : if (ig1==ig2) em1ggp=em1ggp+one
790 : 0 : em1q(ig1,ig2,io)=em1ggp !*Vcp%vc_sqrt(ig1,iqibz)*Vcp%vc_sqrt(ig2,iqibz)
791 : : end do !ig1
792 : : end do !ig2
793 : : end do !iomega
794 : :
795 : : case default
796 [ # # # # ]: 0 : ABI_BUG(sjoin('Wrong ppm%model:',itoa(ppm%model)))
797 : : end select
798 : :
799 : 0 : end subroutine ppm_getem1
800 : : !!***
801 : :
802 : : !----------------------------------------------------------------------
803 : :
804 : : !!****f* m_ppmodel/ppm_getem1_one_ggp
805 : : !! NAME
806 : : !! ppm_getem1_one_ggp
807 : : !!
808 : : !! FUNCTION
809 : : !! Same as ppm_getem1, but does it for a single set of G,G' vectors
810 : : !!
811 : : !! INPUTS
812 : : !!
813 : : !! OUTPUT
814 : : !!
815 : : !! SOURCE
816 : :
817 : 0 : subroutine ppm_getem1_one_ggp(ppm, iqibz, zcut, nomega, omega, Vcp, em1q, ig1, ig2)
818 : :
819 : : !Arguments ------------------------------------
820 : : !scalars
821 : : class(ppmodel_t),intent(in) :: ppm
822 : : integer,intent(in) :: iqibz,nomega
823 : : type(vcoul_t),intent(in) :: Vcp
824 : : real(dp),intent(in) :: zcut
825 : : integer, intent(in) :: ig1,ig2
826 : : !arrays
827 : : complex(dp),intent(in) :: omega(nomega)
828 : : complex(dp),intent(out) :: em1q(nomega)
829 : :
830 : : !Local variables-------------------------------
831 : : !scalars
832 : : integer :: io,idm !,ig1_min,ig2_min,ig2_max
833 : : real(dp) :: den
834 : : complex(dp) :: qpg1,qpg2,ug1,ug2
835 : : complex(dp) :: delta,em1ggp,otw,zzpq,yg1,yg2,bot1,bot2,chig1g2
836 : : !character(len=500) :: msg
837 : : ! *************************************************************************
838 : :
839 [ # # ]: 0 : ABI_CHECK(ppm%mqmem /= 0, 'mqmem==0 not implemented')
840 : :
841 : : !TODO zcut should be an entry in ppm
842 : 0 : delta=CMPLX(zero,zcut)
843 : :
844 : 0 : select case (ppm%model)
845 : :
846 : : case (PPM_GODBY_NEEDS, PPM_HYBERTSEN_LOUIE)
847 [ # # ]: 0 : do io=1,nomega
848 : : !den = omega(io)**2-REAL(ppm%omegatw(iqibz)%vals(ig1,ig2)**2)
849 : : !if (den**2<zcut**2) den = omega(io)**2-REAL( (ppm%omegatw(iqibz)%value(ig1,ig2)-delta)**2 )
850 : 0 : den = omega(io)**2-REAL( (ppm%omegatw(iqibz)%vals(ig1,ig2)-delta)**2 )
851 : 0 : em1ggp = ppm%bigomegatwsq(iqibz)%vals(ig1,ig2)/den
852 [ # # ]: 0 : if (ig1==ig2) em1ggp=em1ggp+one
853 : 0 : em1q(io)=em1ggp
854 : : !em1q(io)=em1ggp*Vcp%vc_sqrt(ig1,iqibz)*Vcp%vc_sqrt(ig2,iqibz)
855 : : end do !io
856 : :
857 : : case (PPM_LINDEN_HORSH)
858 : : !TODO Check coefficients
859 [ # # ]: 0 : do io=1,nomega
860 : 0 : em1ggp=czero
861 [ # # ]: 0 : do idm=1,ppm%npwc
862 : : !den=omega(io)**2-(ppm%omegatw(iqibz)%vals(ig1,ig2)-delta)**2
863 : : !em1w(io)=em1w(io)+eigvec(ig1,idm,iqibz)*conjg(eigvec(ig2,idm,iqibz))*bigomegatwsq(ig1,ig2,iqibz)/den
864 : 0 : ug1 =ppm%eigpot(iqibz)%vals(ig1,idm)
865 : 0 : ug2 =ppm%eigpot(iqibz)%vals(ig2,idm)
866 : 0 : otw =ppm%bigomegatwsq(iqibz)%vals(idm,1)*ppm%omegatw(iqibz)%vals(idm,1)
867 : 0 : zzpq=ppm%bigomegatwsq(iqibz)%vals(idm,1)
868 : 0 : den=half*REAL(zzpq*otw*( one/(omega(io)-otw+delta) - one/(omega(io)+otw-delta) ))
869 : 0 : em1ggp=em1ggp+ug1*CONJG(ug2)*den
870 : : !eigenvalues(idm,io)=one + half*REAL(zzpq*otw*( one/(omega(io)-otw+delta) - one/(omega(io)+otw-delta) ))
871 : : end do
872 : :
873 [ # # ]: 0 : if (ig2==ig1) em1ggp=em1ggp+one
874 : 0 : em1q(io)=em1ggp
875 : : end do !iomega
876 : :
877 : : case (PPM_ENGEL_FARID)
878 : : ! Make e^-1
879 [ # # ]: 0 : do io=1,nomega
880 : 0 : qpg2=one/Vcp%vc_sqrt(ig2,iqibz)
881 : 0 : qpg1=one/Vcp%vc_sqrt(ig1,iqibz)
882 : :
883 : 0 : chig1g2=czero
884 [ # # ]: 0 : do idm=1,ppm%npwc
885 : 0 : otw =ppm%omegatw(iqibz)%vals(idm,1)
886 : 0 : bot1=ppm%bigomegatwsq(iqibz)%vals(ig1,idm)
887 : 0 : bot2=ppm%bigomegatwsq(iqibz)%vals(ig2,idm)
888 : 0 : yg1=SQRT(otw/four_pi)*qpg1*bot1
889 : 0 : yg2=SQRT(otw/four_pi)*qpg2*bot2
890 : 0 : chig1g2=chig1g2 + yg1*CONJG(yg2)/(omega(io)**2-(otw-delta)**2)
891 : : end do
892 : :
893 : 0 : em1ggp=four_pi*chig1g2/(qpg1*qpg2)
894 [ # # ]: 0 : if (ig1==ig2) em1ggp=em1ggp+one
895 : 0 : em1q(io)=em1ggp !*Vcp%vc_sqrt(ig1,iqibz)*Vcp%vc_sqrt(ig2,iqibz)
896 : : end do ! io
897 : :
898 : : case default
899 [ # # # # ]: 0 : ABI_BUG(sjoin('Wrong ppm%model:',itoa(ppm%model)))
900 : : end select
901 : :
902 : 0 : end subroutine ppm_getem1_one_ggp
903 : : !!***
904 : :
905 : : !----------------------------------------------------------------------
906 : :
907 : : !!****f* m_ppmodel/ppm_get_eigenvalues
908 : : !! NAME
909 : : !! ppm_get_eigenvalues
910 : : !!
911 : : !! FUNCTION
912 : : !! Constructs the inverse dielectri matrixc starting from the plasmon-pole
913 : : !! parameters and calculates the frequency-dependent eigenvalues for each
914 : : !! of the nomega frequencies specified in the array omega.
915 : : !!
916 : : !! INPUTS
917 : : !!
918 : : !! OUTPUT
919 : : !!
920 : : !! SOURCE
921 : :
922 : 0 : subroutine ppm_get_eigenvalues(ppm, iqibz, zcut, nomega, omega, Vcp, eigenvalues)
923 : :
924 : : !Arguments ------------------------------------
925 : : !scalars
926 : : class(ppmodel_t),intent(in) :: ppm
927 : : integer,intent(in) :: iqibz,nomega
928 : : type(vcoul_t),intent(in) :: Vcp
929 : : real(dp),intent(in) :: zcut
930 : : !arrays
931 : : complex(dp),intent(in) :: omega(nomega)
932 : : complex(dp),intent(out) :: eigenvalues(ppm%npwc,nomega)
933 : :
934 : : !Local variables-------------------------------
935 : : !scalars
936 : : integer :: info,lwork,negw,ig1,ig2,idx,sdim,iomega,ierr
937 : : character(len=500) :: msg
938 : : !arrays
939 : 0 : real(dp),allocatable :: ww(:),rwork(:)
940 : 0 : complex(dp),allocatable :: work(:),Adpp(:),eigvec(:,:),wwc(:),vs(:,:),Afull(:,:)
941 : 0 : complex(dp),allocatable :: em1q(:,:,:)
942 : 0 : logical,allocatable :: bwork(:)
943 : : logical :: sortcplx !BUG in abilint
944 : : ! *************************************************************************
945 : :
946 [ # # ]: 0 : ABI_CHECK(ppm%mqmem/=0,'mqmem==0 not implemented')
947 : :
948 [ # # # # : 0 : ABI_MALLOC(em1q, (ppm%npwc,ppm%npwc,nomega))
# # # # #
# # # # #
# # # # #
# ]
949 : :
950 : 0 : call ppm%getem1(ppm%npwc,iqibz,zcut,nomega,omega,Vcp,em1q)
951 : :
952 [ # # ]: 0 : do iomega=1,nomega
953 [ # # ]: 0 : if (ABS(REAL(omega(iomega)))>0.00001) then
954 : : ! Eigenvalues for a generic complex matrix.
955 : :
956 : 0 : lwork=4*2*ppm%npwc
957 [ # # # # : 0 : ABI_MALLOC(wwc,(ppm%npwc))
# # # # #
# # # ]
958 [ # # # # : 0 : ABI_MALLOC(work,(lwork))
# # # # #
# # # ]
959 [ # # # # : 0 : ABI_MALLOC(rwork,(ppm%npwc))
# # # # #
# # # ]
960 [ # # # # : 0 : ABI_MALLOC(bwork,(ppm%npwc))
# # # # #
# # # ]
961 [ # # # # : 0 : ABI_MALLOC(vs,(ppm%npwc,ppm%npwc))
# # # # #
# # # # #
# # ]
962 [ # # # # : 0 : ABI_MALLOC(Afull,(ppm%npwc,ppm%npwc))
# # # # #
# # # ]
963 [ # # # # : 0 : Afull=em1q(:,:,iomega)
# # # # #
# ]
964 : :
965 : : !for the time being, no sorting. Maybe here I should sort using the real part?
966 : 0 : call ZGEES('V','N',sortcplx,ppm%npwc,Afull,ppm%npwc,sdim,wwc,vs,ppm%npwc,work,lwork,rwork,bwork,info)
967 [ # # ]: 0 : if (info/=0) then
968 : 0 : write(msg,'(2a,i10)')' ppm_get_eigenvalues: Error in ZGEES, diagonalizing complex matrix, info = ',info
969 : 0 : call wrtout(std_out,msg)
970 : : end if
971 : :
972 [ # # ]: 0 : eigenvalues(:,iomega)=wwc(:)
973 : :
974 [ # # ]: 0 : ABI_FREE(wwc)
975 [ # # ]: 0 : ABI_FREE(work)
976 [ # # ]: 0 : ABI_FREE(rwork)
977 [ # # ]: 0 : ABI_FREE(bwork)
978 [ # # ]: 0 : ABI_FREE(vs)
979 [ # # ]: 0 : ABI_FREE(Afull)
980 : :
981 : : else
982 : : ! === Hermitian Case ===
983 : 0 : lwork=2*ppm%npwc-1
984 [ # # # # : 0 : ABI_MALLOC(ww,(ppm%npwc))
# # # # #
# # # ]
985 [ # # # # : 0 : ABI_MALLOC(work,(lwork))
# # # # #
# # # ]
986 [ # # # # : 0 : ABI_MALLOC(rwork,(3*ppm%npwc-2))
# # # # #
# # # #
# ]
987 [ # # # # : 0 : ABI_MALLOC(eigvec,(ppm%npwc,ppm%npwc))
# # # # #
# # # # #
# # ]
988 : :
989 [ # # # # : 0 : ABI_MALLOC_OR_DIE(Adpp,(ppm%npwc*(ppm%npwc+1)/2), ierr)
# # # # #
# # # # #
# # ]
990 : : !write(std_out,*) 'in hermitian'
991 : :
992 : 0 : idx=0
993 [ # # ]: 0 : do ig2=1,ppm%npwc
994 [ # # ]: 0 : do ig1=1,ig2
995 : 0 : idx=idx+1
996 : 0 : Adpp(idx)=em1q(ig1,ig2,iomega)
997 : : end do
998 : : end do
999 : :
1000 : : ! Require eigenvectors as well
1001 : 0 : call ZHPEV('V','U',ppm%npwc,Adpp,ww,eigvec,ppm%npwc,work,rwork,info)
1002 : :
1003 [ # # ]: 0 : ABI_CHECK(info == 0, sjoin('Error diagonalizing matrix, info: ', itoa(info)))
1004 : :
1005 [ # # # # ]: 0 : negw = (COUNT((REAL(ww)<tol6)))
1006 [ # # ]: 0 : if (negw /= 0) then
1007 [ # # # # : 0 : write(msg,'(a,i0,a,i0,a,f8.4)')'Found negative eigenvalues. No. ',negw,' at iqibz= ',iqibz,' minval= ',MINVAL(REAL(ww))
# # # # ]
1008 : 0 : ABI_WARNING(msg)
1009 : : end if
1010 : :
1011 [ # # ]: 0 : eigenvalues(:,iomega)=ww(:)
1012 : :
1013 [ # # ]: 0 : ABI_FREE(ww)
1014 [ # # ]: 0 : ABI_FREE(work)
1015 [ # # ]: 0 : ABI_FREE(rwork)
1016 [ # # ]: 0 : ABI_FREE(eigvec)
1017 [ # # ]: 0 : ABI_FREE(Adpp)
1018 : : end if
1019 : : end do !iomega
1020 : :
1021 [ # # ]: 0 : ABI_FREE(em1q)
1022 : :
1023 : 0 : end subroutine ppm_get_eigenvalues
1024 : : !!***
1025 : :
1026 : : !----------------------------------------------------------------------
1027 : :
1028 : : !!****f* m_ppmodel/cppm1par
1029 : : !! NAME
1030 : : !! cppm1par
1031 : : !!
1032 : : !! FUNCTION
1033 : : !! Calculate the plasmon-pole parameters big-omega-twiddle-squared and omega-twiddle from
1034 : : !! epsilon-twiddle^-1 calculated for nomega (usually 2) frequencies omega=0 and omega=iE0.
1035 : : !!
1036 : : !! INPUTS
1037 : : !! epsm1(npwc,npwc,nomega)=dielectric matrix at nomega frequencies.
1038 : : !! npwc=number of plane waves
1039 : : !! nomega=number of frequencies (usually 2)
1040 : : !! omega(nomega)=frequencies
1041 : : !! omegaplasma=input variable or Drude plasma frequency
1042 : : !!
1043 : : !! OUTPUT
1044 : : !! bigomegatwsq(npwc,npwc)=parameter of the plasmon-pole model (see gwa.pdf file)
1045 : : !! omegatw(npwc,npwc)=parameter of the plasmon-pole model (see gwa.pdf file)
1046 : : !!
1047 : : !! TODO
1048 : : !! Calculation can be done in place.
1049 : : !!
1050 : : !! SOURCE
1051 : :
1052 : 1085 : subroutine cppm1par(npwc, nomega, omega, omegaplasma, epsm1, omegatw, bigomegatwsq)
1053 : :
1054 : : !Arguments ------------------------------------
1055 : : !scalars
1056 : : integer,intent(in) :: nomega,npwc
1057 : : real(dp),intent(in) :: omegaplasma
1058 : : !arrays
1059 : : complex(dp),intent(in) :: omega(nomega)
1060 : : complex(gwp),intent(in) :: epsm1(npwc,npwc,nomega)
1061 : : complex(gwp),intent(out) :: omegatw(npwc,npwc), bigomegatwsq(npwc,npwc)
1062 : :
1063 : : !Local variables-------------------------------
1064 : : !scalars
1065 : : integer :: ig,igp,io,io0,ioe0
1066 : : real(dp) :: e0,minomega
1067 : : character(len=500) :: msg
1068 : : complex(gwp) :: AA,omegatwsq,diff,ratio,epsm1_io0,epsm1_ioe0
1069 : : ! *************************************************************************
1070 : :
1071 : : ! Find omega=0 and omega=imag (closest to omegaplasma) to fit the ppm parameters
1072 : 1085 : minomega=1.0d-3; io0=0
1073 [ + + ]: 3406 : do io=1,nomega
1074 [ + + ]: 3406 : if (ABS(omega(io))<minomega) then
1075 : 1085 : io0=io; minomega=ABS(omega(io))
1076 : : end if
1077 : : end do
1078 [ - + ]: 1085 : ABI_CHECK(io0 /= 0, "omega=0 not found")
1079 : :
1080 : 1085 : minomega=1.0d-3; e0=200.0; ioe0=0
1081 [ + + ]: 3406 : do io=1,nomega
1082 [ + + + + ]: 3406 : if (REAL(omega(io))<minomega.and.AIMAG(omega(io))>minomega) then
1083 [ + + ]: 1145 : if (ABS(AIMAG(omega(io))-omegaplasma)<ABS(e0-omegaplasma)) then
1084 : 2321 : ioe0=io; e0=AIMAG(omega(io))
1085 : : end if
1086 : : end if
1087 : : end do
1088 : :
1089 : 1085 : write(msg,'(a,f9.4,a)')' Imaginary frequency for fit located at: ',e0*Ha_eV,' [eV] '
1090 : 1085 : call wrtout(std_out, msg)
1091 [ - + ]: 1085 : ABI_CHECK(ioe0 /= 0,"Imaginary omega not found")
1092 : :
1093 : : ! ================================================================
1094 : : ! === Calculate plasmon-pole A parameter A=epsilon^-1(0)-delta ===
1095 : : ! ================================================================
1096 [ + + ]: 53890 : do ig=1,npwc
1097 [ + + ]: 4140407 : do igp=1,npwc
1098 : 4086517 : epsm1_io0 = epsm1(ig,igp,io0)
1099 : 4086517 : epsm1_ioe0 = epsm1(ig,igp,ioe0)
1100 : :
1101 : 4086517 : AA=epsm1_io0
1102 [ + + ]: 4086517 : if (ig==igp) AA=AA-one
1103 : :
1104 : : ! === Calculate plasmon-pole omega-twiddle-square parameter ===
1105 : : ! XG201009 Strangely, the next formula does not work with gcc43-debug
1106 : : ! omegatwsq=(AA/(epsm1_io0-epsm1_ioe0)-one)*e0**2
1107 : : ! This seems to be due to precision issue at the level of division by a complex whose norm squared
1108 : : ! is below the smallest representable number.
1109 : : ! After many trials, I have decided to shift the difference by a small number ... well, not so small ...
1110 : : ! for numerical issues
1111 : 4086517 : diff=epsm1_io0-epsm1_ioe0
1112 : 4086517 : diff=diff+cmplx(tol10,tol10)
1113 : 4086517 : ratio=AA/diff
1114 : 4086517 : omegatwsq=(ratio-cone)*e0**2
1115 : : !
1116 : : ! If omega-twiddle-squared is negative,set omega-twiddle-squared to 1.0 (a reasonable way of treating
1117 : : ! such terms, in which epsilon**-1 was originally increasing along this part of the imaginary axis)
1118 : : ! (note: originally these terms were ignored in Sigma; this was changed on 6 March 1990.)
1119 : :
1120 [ + + ]: 4086517 : if (REAL(omegatwsq)<=zero) omegatwsq=one
1121 : : !
1122 : : ! Get omega-twiddle. Neglect the imag part (if any) in omega-twiddle-squared
1123 : 4086517 : omegatw(ig,igp)=SQRT(REAL(omegatwsq))
1124 : :
1125 : : ! Get big-omega-twiddle-squared=-omega-twiddle-squared AA
1126 : 4139322 : bigomegatwsq(ig,igp)=-AA*omegatw(ig,igp)**2
1127 : :
1128 : : end do !igp
1129 : : end do !ig
1130 : :
1131 : 1085 : write(msg,'(2a,f15.12,2a,2i5,a)')ch10,&
1132 [ + + + + : 4140407 : ' cppm1par : omega twiddle minval [eV] = ',MINVAL(ABS(omegatw))*Ha_eV,ch10,&
+ + + - -
+ - - ]
1133 [ + - + + : 4143662 : ' omega twiddle min location = ',MINLOC(ABS(omegatw)),ch10
+ + + + ]
1134 : 1085 : call wrtout(std_out,msg)
1135 : :
1136 : 1085 : end subroutine cppm1par
1137 : : !!***
1138 : :
1139 : : !----------------------------------------------------------------------
1140 : :
1141 : : !!****f* m_ppmodel/cppm2par
1142 : : !! NAME
1143 : : !! cppm2par
1144 : : !!
1145 : : !! FUNCTION
1146 : : !! Calculate plasmon-pole parameters of the Hybertsen and Louie model (PRB 34, 5390 (1986) [[cite:Hybertsen1986]])
1147 : : !!
1148 : : !! INPUTS
1149 : : !! qpt(3)=The coordinates of the q-point in the IBZ.
1150 : : !! epsm1(npwc,npwc)=symmetrized inverse dielectric (static limit is used)
1151 : : !! gmet(3,3)=metric in reciprocal space
1152 : : !! ngfftf(18)=contain all needed information about the 3D fine FFT, see ~abinit/doc/variables/vargs.htm#ngfft
1153 : : !! npwc=number of plane waves in epsm1
1154 : : !! rhor(nfftf)=charge density on the real space FFT grid
1155 : : !! nfftf= total number of points in the fine FFT mesh (for this processor)
1156 : : !! invalid_freq: what to do when PPM omega is found negative or imaginary
1157 : : !! 0) drop it (default as specified in Hybersen-Louie original GW paper)
1158 : : !! 1) set to 1 hartree
1159 : : !! 2) set to infinity
1160 : : !!
1161 : : !! OUTPUT
1162 : : !! bigomegatwsq(npwc,npwc)= squared bare plasma frequencies
1163 : : !! \Omega^2_{G1 G2}(q) = 4\pi \frac {(q+G1).(q+G2)}/{|q+G1|^2} n(G1-G2)
1164 : : !! omegatw(npwc,npwc)= plasmon frequencies \tilde\omega_{G1 G2}(q) where:
1165 : : !! \tilde\omega^2_{G1 G2}(q) =
1166 : : !! \frac {\Omega^2_{G1 G2}(q)} {\delta_{G1 G2}-\tilde\epsilon^{-1}_{G1 G2} (q, \omega=0)}
1167 : : !!
1168 : : !! SOURCE
1169 : :
1170 : 24 : subroutine cppm2par(qpt, npwc, epsm1, ngfftf, gvec, gprimd, rhor, nfftf, gmet, bigomegatwsq, omegatw, invalid_freq)
1171 : :
1172 : : !Arguments ------------------------------------
1173 : : !scalars
1174 : : integer,intent(in) :: npwc,nfftf,invalid_freq
1175 : : !arrays
1176 : : integer,intent(in) :: gvec(3,npwc), ngfftf(18)
1177 : : real(dp),intent(in) :: qpt(3),gmet(3,3),gprimd(3,3), rhor(nfftf)
1178 : : complex(gwp),intent(in) :: epsm1(npwc,npwc)
1179 : : complex(gwp),intent(out) :: bigomegatwsq(npwc,npwc), omegatw(npwc,npwc)
1180 : :
1181 : : !Local variables-------------------------------
1182 : : !scalars
1183 : : integer :: ig,igp,nimwp,ngfft1,ngfft2,ngfft3,gmgp_idx,ierr
1184 : : real(dp) :: lambda,phi,AA
1185 : : logical,parameter :: use_symmetrized=.TRUE., check_imppf=.FALSE.
1186 : : character(len=500) :: msg
1187 : 24 : type(MPI_type) :: MPI_enreg_seq
1188 : : !arrays
1189 : : real(dp) :: qlist(3,1)
1190 : 24 : real(dp),allocatable :: tmp_rhor(:),qratio(:,:),qplusg(:),rhog_dp(:,:)
1191 : 24 : complex(gwp),allocatable :: omegatwsq(:,:)
1192 : 24 : complex(gwp),allocatable :: rhog(:),rhogg(:,:),temp(:,:) !MG these should be double precision TODO
1193 : : !*************************************************************************
1194 : :
1195 : 24 : call initmpi_seq(MPI_enreg_seq)
1196 : 24 : call MPI_enreg_seq%distribfft%init_seq('c',ngfftf(2),ngfftf(3),'all')
1197 : :
1198 : : ! Calculate qratio(npwec,npvec) = (q+G).(q+Gp)/|q+G|^2 ===
1199 [ - + + - : 96 : ABI_MALLOC_OR_DIE(qratio,(npwc,npwc), ierr)
- + + - +
- + - + -
- + - + ]
1200 : :
1201 : 24 : call cqratio(npwc,gvec,qpt,gmet,gprimd,qratio)
1202 : : !
1203 : : ! Compute the density in G space rhor(R)--> rhog(G)
1204 [ - + + - : 72 : ABI_MALLOC(rhog_dp,(2,nfftf))
+ - + - -
+ - + ]
1205 [ - + + - : 72 : ABI_MALLOC(rhog,(nfftf))
+ - + - -
+ - + ]
1206 : 24 : ngfft1=ngfftf(1); ngfft2=ngfftf(2); ngfft3=ngfftf(3)
1207 : :
1208 [ - + + - : 72 : ABI_MALLOC(tmp_rhor,(nfftf))
+ - + - -
+ - + ]
1209 [ - + + + ]: 183222 : tmp_rhor = rhor ! To avoid having to use intent(inout).
1210 : 24 : call fourdp(1,rhog_dp,tmp_rhor,-1,MPI_enreg_seq,nfftf,1,ngfftf,0)
1211 [ - + ]: 24 : ABI_FREE(tmp_rhor)
1212 : :
1213 [ + + ]: 183198 : rhog(1:nfftf)=CMPLX(rhog_dp(1,1:nfftf),rhog_dp(2,1:nfftf))
1214 : :
1215 : : ! Calculate the FFT index of each (G-Gp) vector and assign
1216 : : ! the value of the correspondent density simultaneously
1217 [ - + + - : 96 : ABI_MALLOC_OR_DIE(rhogg,(npwc, npwc), ierr)
- + + - +
- + - + -
- + - + ]
1218 : :
1219 : 24 : ierr=0
1220 [ + + ]: 708 : do ig=1,npwc
1221 [ + + ]: 20364 : do igp=1,npwc
1222 [ + + ]: 78624 : gmgp_idx = g2ifft(gvec(:,ig)-gvec(:,igp),ngfftf)
1223 [ + - ]: 20340 : if (gmgp_idx/=0) then
1224 : 19656 : rhogg(ig,igp)=rhog(gmgp_idx)
1225 : : else
1226 : 0 : ierr=ierr+1
1227 : 0 : rhogg(ig,igp)=czero
1228 : : end if
1229 : : end do
1230 : : end do
1231 : :
1232 [ - + ]: 24 : if (ierr /= 0) then
1233 : : write(msg,'(a,i0,1x,3a)')&
1234 : 0 : 'Found ',ierr,' G1-G2 vectors falling outside the FFT box. ',ch10,&
1235 : 0 : 'Enlarge the FFT mesh to get rid of this problem. '
1236 : 0 : ABI_WARNING(msg)
1237 : : end if
1238 : :
1239 [ + + + + ]: 20364 : rhogg=four_pi*rhogg
1240 [ - + ]: 24 : ABI_FREE(rhog_dp)
1241 [ - + ]: 24 : ABI_FREE(rhog)
1242 : :
1243 : : ! Calculate GPP parameters
1244 : : ! unsymmetrized epsm1 -> epsm1=|q+Gp|/|q+G|*epsm1
1245 [ - + + - : 72 : ABI_MALLOC(qplusg,(npwc))
+ - + - -
+ - + ]
1246 [ - + + - : 72 : ABI_MALLOC(temp,(npwc,npwc))
- + + - -
+ - + ]
1247 [ - + + - : 72 : ABI_MALLOC_OR_DIE(omegatwsq,(npwc,npwc), ierr)
- + + - +
- - + -
+ ]
1248 : :
1249 [ - + + + : 20388 : temp = -epsm1(:,:)
+ + ]
1250 : : !
1251 : : ! RS still not obvious for me whether one shall use the symmetrized inverse DM or the unsymmetrized one
1252 : : ! the default here is to use the symmetrized one, I must discuss this with XG
1253 : : !
1254 : : ! MG it turns out that using the symmetrized inverse DM in the plasmon-pole
1255 : : ! equations give the same results for the squared plasmon frequencies omegatwsq while the
1256 : : ! squared bare plasma frequencies bigomegatwsq related to the symmetrized dielectric matrix
1257 : : ! are obtained multiplying by |q+G1|/|q+G2|.
1258 : : !
1259 : : if (.not.use_symmetrized) then
1260 : : qlist(:,1) = qpt
1261 : : call cmod_qpg(1,1,qlist,npwc,gvec,gprimd,qplusg) !MG TODO here take care of small q
1262 : : do ig=1,npwc
1263 : : do igp=1,npwc
1264 : : temp(ig,igp)=qplusg(igp)/qplusg(ig)*temp(ig,igp)
1265 : : end do
1266 : : end do
1267 : : end if
1268 : :
1269 : 24 : nimwp=0
1270 [ + + ]: 708 : do ig=1,npwc
1271 : 684 : temp(ig,ig)=temp(ig,ig)+one
1272 [ + + ]: 20364 : do igp=1,npwc
1273 : 19656 : bigomegatwsq(ig,igp) = rhogg(ig,igp)*qratio(ig,igp)
1274 : 19656 : omegatwsq(ig,igp)=bigomegatwsq(ig,igp)/temp(ig,igp)
1275 : : !
1276 : : ! Set to an arbitrary value the omegawsq which become negative or imaginary
1277 : : ! in principle these correspond to cases where the imaginary part of epsm1 does not have
1278 : : ! a well defined peak. The imaginary part of epsm1 in these cases oscillates with a small amplitude
1279 : : ! since the amplitude A_GGpr=-pi/2*bigomegatwsq/omegatw,
1280 : : ! it follows that bigomegatwsq shall be set to zero for these cases
1281 [ + + - + ]: 20340 : if ( REAL(omegatwsq(ig,igp))<= tol12 .or. AIMAG(omegatwsq(ig,igp))**2*tol12> REAL(omegatwsq(ig,igp))**2) then
1282 : 10280 : nimwp=nimwp+1
1283 : :
1284 [ - + ]: 10280 : if ( invalid_freq == 1 ) then
1285 : : ! set omegatwsq to 1 hartree
1286 : 0 : omegatwsq(ig,igp)=cone
1287 : 0 : AA = epsm1(ig,igp)
1288 [ # # ]: 0 : if ( ig == igp ) AA = AA - one
1289 : 0 : omegatw(ig,igp)=SQRT(REAL(omegatwsq(ig,igp)))
1290 : 0 : bigomegatwsq(ig,igp)=-AA*omegatw(ig,igp)**2
1291 [ + + ]: 10280 : elseif ( invalid_freq == 2 ) then
1292 : : ! set omegatwsq to infinity
1293 : 3911 : omegatwsq(ig,igp)=cone/tol6
1294 : 3911 : AA = epsm1(ig,igp)
1295 [ - + ]: 3911 : if ( ig == igp ) AA = AA - one
1296 : 3911 : omegatw(ig,igp)=SQRT(REAL(omegatwsq(ig,igp)))
1297 : 3911 : bigomegatwsq(ig,igp)=-AA*omegatw(ig,igp)**2
1298 : : else
1299 : : ! simply ignore all cases of omegatw with imaginary values
1300 : 6369 : bigomegatwsq(ig,igp)=(0.,0.)
1301 : 6369 : omegatw(ig,igp)=(ten,0.)
1302 : : end if
1303 : : if (check_imppf) then
1304 : : write(msg,'(a,2(i0,1x))')' Imaginary plasmon frequency at : ',ig,igp
1305 : : call wrtout(std_out,msg)
1306 : : end if
1307 : : else
1308 : : ! this part has been added to deal with systems without inversion symmetry
1309 : : ! this new implementation gives the same results as the previous one if
1310 : : ! omegatwsq is a pure real number and has the advantage of being an improved
1311 : : ! approach for systems without an inversion center.
1312 : 9376 : lambda=ABS(omegatwsq(ig,igp))
1313 : 9376 : phi=ATAN(AIMAG(omegatwsq(ig,igp))/REAL(omegatwsq(ig,igp)))
1314 : 9376 : omegatw(ig,igp)=SQRT(lambda/COS(phi))
1315 : 9376 : bigomegatwsq(ig,igp)=bigomegatwsq(ig,igp)*(1.-(0.,1.)*TAN(phi))
1316 : : ! Uncomment the following line and comment the previous to restore the old version.
1317 : : !omegatw(ig,igp)=sqrt(real(omegatwsq(ig,igp)))
1318 : : end if
1319 : : end do
1320 : : end do
1321 : :
1322 [ + - ]: 24 : write(msg,'(3a,i0,a,i0)')' At q-point : ',trim(ktoa(qpt)), ' # imaginary plasmonpole frequencies: ',nimwp,' / ',npwc**2
1323 : 24 : call wrtout(std_out, msg)
1324 : : write(msg,'(a,f12.8,a,3(i0,1x))') &
1325 [ + + + + : 40752 : " omega twiddle minval: ", MINVAL(ABS(omegatw))*Ha_eV, "[eV], min location: ",MINLOC(ABS(omegatw))
+ + + - -
+ - - + -
+ + + + +
+ ]
1326 : 24 : call wrtout(std_out, msg)
1327 : :
1328 : 24 : call destroy_mpi_enreg(MPI_enreg_seq)
1329 : :
1330 [ - + ]: 24 : ABI_FREE(omegatwsq)
1331 [ - + ]: 24 : ABI_FREE(rhogg)
1332 [ - + ]: 24 : ABI_FREE(temp)
1333 [ - + ]: 24 : ABI_FREE(qplusg)
1334 [ - + ]: 24 : ABI_FREE(qratio)
1335 : :
1336 [ - + - + : 24 : end subroutine cppm2par
- + - + -
+ - + - +
- + - + -
+ - + ]
1337 : : !!***
1338 : :
1339 : : !----------------------------------------------------------------------
1340 : :
1341 : : !!****f* m_ppmodel/cppm3par
1342 : : !! NAME
1343 : : !! cppm3par
1344 : : !!
1345 : : !! FUNCTION
1346 : : !! Calculate the plasmon-pole parameters using the von Linden-Horsh model (PRB 37, 8351, 1988) [[cite:vonderLinden1988]]
1347 : : !! (see also Pag 22 of Quasiparticle Calculations in Solids [[cite:Aulbur2001]].
1348 : : !!
1349 : : !! INPUTS
1350 : : !! epsm1(npwc,npwc))= symmetrized inverse dielectric
1351 : : !! ngfftf(18)=contain all needed information about 3D fine FFT, see ~abinit/doc/variables/vargs.htm#ngfft
1352 : : !! npwc=number of plane waves in epsm1
1353 : : !! qratio=(q+G1).(q+G2)/(|q+G1|.|q+G2|)
1354 : : !! rhor(nfftf)=charge density on the real space FFT grid
1355 : : !! nfftf=number of points in the FFT grid (for this processor)
1356 : : !! gvec(3,npwc)= G vectors in reduced coordinates
1357 : : !!
1358 : : !! OUTPUT
1359 : : !! omegatw(npwc,npwc)= plasmon pole positions
1360 : : !! bigomegatwsq(npwc,npwc)=(E_{q,ii}^{-1}-1)*omegatw
1361 : : !! where E^{-1} is the eigenvalue of the inverse dielectric matrix
1362 : : !! eigtot(npwc,npwc)=the eigvectors of the symmetrized inverse dielectric matrix
1363 : : !! (first index for G, second index for bands)
1364 : : !!
1365 : : !! SOURCE
1366 : :
1367 : 18 : subroutine cppm3par(qpt,npwc,epsm1,ngfftf,gvec,gprimd,rhor,nfftf,bigomegatwsq,omegatw,eigtot)
1368 : :
1369 : : !Arguments ------------------------------------
1370 : : !scalars
1371 : : integer,intent(in) :: nfftf,npwc
1372 : : !arrays
1373 : : integer,intent(in) :: gvec(3,npwc),ngfftf(18)
1374 : : real(dp),intent(in) :: qpt(3),gprimd(3,3),rhor(nfftf)
1375 : : complex(gwp),intent(in) :: epsm1(npwc,npwc)
1376 : : complex(gwp),intent(out) :: bigomegatwsq(npwc,1),omegatw(npwc) ,eigtot(npwc,npwc)
1377 : :
1378 : : !Local variables-------------------------------
1379 : : !TODO these should be dp
1380 : : !scalars
1381 : : integer :: idx,ierr,ig,igp,ii,jj,ngfft1,ngfft2,ngfft3,gmgp_idx
1382 : : real(dp) :: num,qpg_dot_qpgp
1383 : : complex(dp) :: conjg_eig
1384 : : logical :: qiszero
1385 : : character(len=500) :: msg
1386 : 18 : type(MPI_type) :: MPI_enreg_seq
1387 : : !arrays
1388 : : real(dp) :: b1(3),b2(3),b3(3),gppq(3),gpq(3),qlist(3,1)
1389 : 18 : real(dp),allocatable :: eigval(:),qplusg(:),rhog_dp(:,:),zhpev2(:),tmp_rhor(:)
1390 : 18 : complex(dp),allocatable :: eigvec(:,:),matr(:),mm(:,:),rhog(:),rhogg(:,:), zhpev1(:),zz(:)
1391 : : !*************************************************************************
1392 : :
1393 : : ! Fake MPI_type for the sequential part.
1394 : 18 : call initmpi_seq(MPI_enreg_seq)
1395 : 18 : call MPI_enreg_seq%distribfft%init_seq('c',ngfftf(2),ngfftf(3),'all')
1396 : :
1397 [ + + + + ]: 30 : qiszero = (ALL(ABS(qpt)<1.0e-3))
1398 : :
1399 [ + + + + : 180 : b1 = two_pi*gprimd(:,1); b2 = two_pi*gprimd(:,2); b3 = two_pi*gprimd(:,3)
+ + ]
1400 : :
1401 : 18 : ngfft1=ngfftf(1); ngfft2=ngfftf(2); ngfft3=ngfftf(3)
1402 : :
1403 [ - + + - : 54 : ABI_MALLOC(rhog_dp,(2,nfftf))
+ - + - -
+ - + ]
1404 [ - + + - : 54 : ABI_MALLOC(rhog,(nfftf))
+ - + - -
+ - + ]
1405 [ - + + - : 72 : ABI_MALLOC_OR_DIE(rhogg,(npwc, npwc), ierr)
- + + - +
- + - + -
- + - + ]
1406 : : !
1407 : : ! === Compute the density in G space rhog(r)--> rho(G) ===
1408 : : ! FIXME this has to be fixed, rho(G) should be passed instead of doing FFT for each q
1409 : :
1410 [ - + + - : 54 : ABI_MALLOC(tmp_rhor,(nfftf))
+ - + - -
+ - + ]
1411 [ - + + + ]: 65112 : tmp_rhor=rhor ! To avoid having to use intent(inout).
1412 : 18 : call fourdp(1,rhog_dp,tmp_rhor,-1,MPI_enreg_seq,nfftf,1,ngfftf,0)
1413 [ - + ]: 18 : ABI_FREE(tmp_rhor)
1414 : :
1415 [ + + ]: 65094 : rhog(1:nfftf)=CMPLX(rhog_dp(1,1:nfftf),rhog_dp(2,1:nfftf))
1416 : : !
1417 : : ! Calculate the FFT index of each (G-Gp) vector and assign the value
1418 : : ! of the correspondent density simultaneously
1419 : 18 : ierr=0
1420 [ + + ]: 504 : do ig=1,npwc
1421 [ + + ]: 13626 : do igp=1,npwc
1422 [ + + ]: 52488 : gmgp_idx = g2ifft(gvec(:,ig)-gvec(:,igp),ngfftf)
1423 [ + - ]: 13608 : if (gmgp_idx/=0) then
1424 : 13122 : rhogg(ig,igp)=rhog(gmgp_idx)
1425 : : else
1426 : 0 : ierr=ierr+1
1427 : 0 : rhogg(ig,igp)=czero
1428 : : end if
1429 : : end do
1430 : : end do
1431 : :
1432 [ - + ]: 18 : if (ierr /= 0) then
1433 : : write(msg,'(a,i0,3a)')&
1434 : 0 : 'Found ',ierr,' G1-G2 vectors falling outside the FFT box. ',ch10,&
1435 : 0 : 'Enlarge the FFT mesh to get rid of this problem. '
1436 : 0 : ABI_WARNING(msg)
1437 : : end if
1438 : :
1439 : : ! mm(G,Gp) = (q+G) \cdot (q+Gp) n(G-Gp)
1440 [ - + + - : 54 : ABI_MALLOC_OR_DIE(mm, (npwc,npwc), ierr)
- + + - +
- - + -
+ ]
1441 : :
1442 [ + + ]: 504 : do ig=1,npwc
1443 [ + + ]: 486 : if (qiszero) then
1444 : : ! To be discussed with Riad, here we should use the small q
1445 : : ! to be consistent and consider the limit q-->0
1446 [ + + ]: 324 : gpq(:)=gvec(:,ig)
1447 : : else
1448 [ + + ]: 1620 : gpq(:)=gvec(:,ig)+qpt
1449 : : end if
1450 [ + + ]: 13626 : do igp=1,npwc
1451 [ + + ]: 13122 : if (qiszero) then
1452 [ + + ]: 8748 : gppq(:)=gvec(:,igp)
1453 : : else
1454 [ + + ]: 43740 : gppq(:)=gvec(:,igp)+qpt
1455 : : end if
1456 : 13122 : qpg_dot_qpgp=zero
1457 [ + + ]: 52488 : do ii=1,3
1458 : : qpg_dot_qpgp=qpg_dot_qpgp+&
1459 : : ( gpq(1)*b1(ii) +gpq(2)*b2(ii) +gpq(3)*b3(ii))*&
1460 : 52488 : (gppq(1)*b1(ii)+gppq(2)*b2(ii)+gppq(3)*b3(ii))
1461 : : end do
1462 : 13608 : mm(ig,igp)=rhogg(ig,igp)*qpg_dot_qpgp
1463 : : end do !igp
1464 : : end do !ig
1465 : :
1466 [ - + ]: 18 : ABI_FREE(rhog_dp)
1467 [ - + ]: 18 : ABI_FREE(rhog)
1468 : : ! === Now we have rhogg,rho0 ===
1469 : : !
1470 : : ! Calculate the dielectric matrix eigenvalues and vectors
1471 : : ! Use only the static epsm1 i.e., only the w=0 part (eps(:,:,1,:))
1472 [ - + + - : 54 : ABI_MALLOC(eigval,(npwc))
+ - + - -
+ - + ]
1473 [ - + + - : 54 : ABI_MALLOC_OR_DIE(eigvec, (npwc, npwc), ierr)
- + + - +
- - + -
+ ]
1474 : :
1475 [ - + + - : 54 : ABI_MALLOC(zz,(npwc))
+ - + - -
+ - + ]
1476 [ + + ]: 504 : zz=czero
1477 : :
1478 [ - + + - : 36 : ABI_MALLOC(qplusg,(npwc))
- + - + ]
1479 : :
1480 : : ! Store the susceptibility matrix in upper mode before calling zhpev.
1481 [ + - - + : 54 : ABI_MALLOC_OR_DIE(matr,(npwc*(npwc+1)/2), ierr)
+ - + - +
- + - - +
- + ]
1482 : :
1483 : : idx=1
1484 [ + + ]: 504 : do ii=1,npwc
1485 [ + + ]: 7308 : do jj=1,ii
1486 : 7290 : matr(idx)=epsm1(jj,ii); idx=idx+1
1487 : : end do
1488 : : end do
1489 : :
1490 [ + - - + : 54 : ABI_MALLOC(zhpev2,(3*npwc-2))
+ - + - +
- - + -
+ ]
1491 [ + - - + : 54 : ABI_MALLOC(zhpev1,(2*npwc-1))
+ - + - +
- - + -
+ ]
1492 : :
1493 : 18 : call ZHPEV('V','U',npwc,matr,eigval,eigvec,npwc,zhpev1,zhpev2,ierr)
1494 [ - + ]: 18 : ABI_FREE(matr)
1495 [ - + ]: 18 : ABI_FREE(zhpev2)
1496 [ - + ]: 18 : ABI_FREE(zhpev1)
1497 : :
1498 [ - + ]: 18 : if (ierr < 0) then
1499 : : write (msg,'(2a,i0,a)')&
1500 : 0 : ' Failed to calculate the eigenvalues and eigenvectors of the dielectric matrix ',ch10,&
1501 : 0 : ierr*(-1),'-th argument in the matrix has an illegal value. '
1502 : 0 : ABI_ERROR(msg)
1503 : : end if
1504 : :
1505 [ - + ]: 18 : if (ierr > 0) then
1506 : : write(msg,'(3a,i0,2a)')&
1507 : 0 : ' Failed to calculate the eigenvalues and eigenvectors of the dielectric matrix ',ch10,&
1508 : 0 : ' the algorithm failed to converge; ierr = ', ierr,ch10,&
1509 : 0 : ' off-diagonal elements of an intermediate tridiagonal form did not converge to zero. '
1510 : 0 : ABI_ERROR(msg)
1511 : : end if
1512 : :
1513 : : ! Calculate the PPM parameters and the eigenpotentials needed for
1514 : : ! the calculation of the generalized overlap matrix
1515 : : ! Note: the eigenpotentials has to be calculated on the FFT (G-Gp) index
1516 : : !
1517 : : ! Save eigenvectors of \tilde\epsilon^{-1}
1518 : : ! MG well it is better to save \Theta otherwise
1519 : : ! we have to calculare \Theta for each band, spin, k-point but oh well
1520 [ + + + + ]: 13626 : eigtot=eigvec
1521 : :
1522 : 18 : qlist(:,1) = qpt
1523 : 18 : call cmod_qpg(1,1,qlist,npwc,gvec,gprimd,qplusg) !MG TODO here take care of small q
1524 : : !
1525 : : ! Basic Equation:
1526 : : !
1527 : : ! \Theta_{q,ii}(G)=\Psi_{q,ii}(G)/|q+G|
1528 : : ! where \Psi_{q,ii}(G) is the eigenvector of \tilde\epsilon^{-1}
1529 : :
1530 : : ! \tilde\omega_{ii,q}^2= 4\pi (1-eigenval(ii,q)))
1531 : : ! \sum_{G,Gp} \Theta^*_{q,ii}(G) (q+G)\cdot(q+Gp) n(G-Gp) \Theta_{q,ii}(Gp)
1532 : :
1533 [ + + ]: 504 : do ii=1,npwc !DM band
1534 : : ! Calculate \Theta_{q,ii}(G)
1535 : : ! why the first element is not modified? if the problem is the small value of qplusg(1)
1536 : : ! we could multiply by sqrt(mod((q+G)(q+G'))) and then add the sing at the end
1537 [ + + ]: 486 : if (qiszero)then
1538 [ + + ]: 2187 : eigvec(2:,ii)=eigvec(2:,ii)/qplusg(2:)
1539 : : else
1540 [ + + ]: 11340 : eigvec(:,ii)=eigvec(:,ii)/qplusg(:)
1541 : : end if
1542 [ + + ]: 13608 : do ig=1,npwc
1543 : 13122 : conjg_eig=CONJG(eigvec(ig,ii))
1544 [ + + ]: 367902 : do igp=1,npwc
1545 [ + + + + ]: 367416 : if(qiszero .and. ig==1 .and. igp==1)then
1546 : 81 : zz(ii)=zz(ii)+conjg_eig*rhogg(ig,igp)*eigvec(igp,ii)
1547 : : else
1548 : 354213 : zz(ii)=zz(ii)+conjg_eig*mm(ig,igp)*eigvec(igp,ii)
1549 : : end if
1550 : : end do
1551 : : end do
1552 : :
1553 : 486 : num=one-eigval(ii)
1554 [ + - ]: 486 : if (num<=zero) then
1555 : : ! here I think we should set bigomegatwsq=0 and omegatw to an arbitrary value
1556 : : ! maybe we can output a warning TO BE discussed with Riad
1557 [ # # ]: 0 : if (ABS(num)<1.0d-4) then
1558 : : num=1.0d-5
1559 : : else
1560 : 0 : ABI_ERROR("One or more imaginary plasmon pole energies")
1561 : : end if
1562 : : end if
1563 : :
1564 : 486 : omegatw(ii)=SQRT(4*pi*REAL(zz(ii))/num)
1565 : : ! this should be \alpha = 2\pi omegatw * (1-eigenval)
1566 : : ! MG check this, in the review I found a factor 2\pi, maybe it is reintroduced later
1567 : 504 : bigomegatwsq(ii,1)=num*omegatw(ii)
1568 : : end do
1569 : :
1570 [ - + ]: 18 : ABI_FREE(rhogg)
1571 [ - + ]: 18 : ABI_FREE(mm)
1572 [ - + ]: 18 : ABI_FREE(eigval)
1573 [ - + ]: 18 : ABI_FREE(zz)
1574 [ - + ]: 18 : ABI_FREE(eigvec)
1575 [ - + ]: 18 : ABI_FREE(qplusg)
1576 : :
1577 : 18 : call destroy_mpi_enreg(MPI_enreg_seq)
1578 : :
1579 : 18 : write(msg,'(2a,f12.8,2a,3i5)')ch10,&
1580 [ + - - + : 522 : ' cppm3par : omega twiddle minval [eV] = ',MINVAL(ABS(omegatw))*Ha_eV,ch10,&
- - + + ]
1581 [ + - - + : 540 : ' omega twiddle min location = ',MINLOC(ABS(omegatw))
- - + + +
+ + + ]
1582 : 18 : call wrtout(std_out,msg)
1583 : :
1584 [ - + - + : 18 : end subroutine cppm3par
- + - + -
+ - + - +
- + - + -
+ - + ]
1585 : : !!***
1586 : :
1587 : : !----------------------------------------------------------------------
1588 : :
1589 : : !!****f* m_ppmodel/cppm4par
1590 : : !! NAME
1591 : : !! cppm4par
1592 : : !!
1593 : : !! FUNCTION
1594 : : !! Calculate the plasmon-pole parameters using Engel-Farid model (PRB47,15931,1993) [[cite:Engel1993]].
1595 : : !! See also Quasiparticle Calculations in Solids [[cite:Aulbur2001]] page. 23.
1596 : : !!
1597 : : !! INPUTS
1598 : : !! qpt(3)=Reduced coordinates of the q-point.
1599 : : !! npwc=number of plane waves in epsm1
1600 : : !! epsm1(npwc,npwc)=symmetrized inverse dielectric matrix.
1601 : : !! ngfftf(18)=contain all needed information about 3D fine FFT, see ~abinit/doc/variables/vargs.htm#ngfft
1602 : : !! gvec(3,npwc)=G vectors in reduced coordinated
1603 : : !! gprimd(3,3)=dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$)
1604 : : !! rhor(nfftf)=charge density on the real space FFT grid
1605 : : !! nfftf=Number of FFT points.
1606 : : !!
1607 : : !! OUTPUT
1608 : : !! bigomegatwsq(npwc,npwc)=plasmon-pole strength.
1609 : : !! omegatw(npwc)=plasmon-pole frequencies.
1610 : : !!
1611 : : !! SOURCE
1612 : :
1613 : 18 : subroutine cppm4par(qpt, npwc, epsm1, ngfftf, gvec, gprimd, rhor, nfftf, bigomegatwsq, omegatw)
1614 : :
1615 : : !Arguments ------------------------------------
1616 : : !scalars
1617 : : integer,intent(in) :: nfftf,npwc
1618 : : !arrays
1619 : : integer,intent(in) :: gvec(3,npwc),ngfftf(18)
1620 : : real(dp),intent(in) :: gprimd(3,3),qpt(3), rhor(nfftf)
1621 : : complex(gwp),intent(in) :: epsm1(npwc,npwc)
1622 : : complex(gwp),intent(out) :: bigomegatwsq(npwc,npwc),omegatw(npwc)
1623 : :
1624 : : !Local variables-------------------------------
1625 : : !scalars
1626 : : integer :: ierr,ig,igp,ii,ngfft1,ngfft2,ngfft3,gmgp_idx
1627 : : real(dp) :: qpg_dot_qpgp
1628 : : character(len=500) :: msg
1629 : : character(len=80) :: bar
1630 : 18 : type(MPI_type) :: MPI_enreg_seq
1631 : : !arrays
1632 : : real(dp) :: b1(3),b2(3),b3(3),gppq(3),gpq(3),qlist(3,1)
1633 : 18 : real(dp),allocatable :: eigval(:),qplusg(:),rhog_dp(:,:),tmp_rhor(:)
1634 : 18 : complex(dp),allocatable :: chi(:,:), mm(:,:),mtemp(:,:),rhog(:), tmp1(:),zz2(:,:)
1635 : : !*************************************************************************
1636 : :
1637 : : ! Calculate density in G space rhog(G)
1638 : : ! FIXME this has to be fixed, rho(G) should be passed instead of doing FFT for each q
1639 : 18 : call initmpi_seq(MPI_enreg_seq)
1640 : 18 : call MPI_enreg_seq%distribfft%init_seq('c',ngfftf(2),ngfftf(3),'all')
1641 : :
1642 [ - + + - : 54 : ABI_MALLOC(rhog_dp, (2,nfftf))
+ - + - -
+ - + ]
1643 : :
1644 : : ! Conduct FFT tho(r)-->rhog(G)
1645 [ - + + - : 54 : ABI_MALLOC(tmp_rhor,(nfftf))
+ - + - -
+ - + ]
1646 [ - + + + ]: 65112 : tmp_rhor = rhor ! To avoid having to use intent(inout).
1647 : 18 : call fourdp(1,rhog_dp,tmp_rhor,-1,MPI_enreg_seq,nfftf,1,ngfftf,0)
1648 : :
1649 [ - + ]: 18 : ABI_FREE(tmp_rhor)
1650 : 18 : call destroy_mpi_enreg(MPI_enreg_seq)
1651 : :
1652 [ - + + - : 54 : ABI_MALLOC(rhog, (nfftf))
+ - + - -
+ - + ]
1653 [ + + ]: 65094 : rhog(1:nfftf)=CMPLX(rhog_dp(1,1:nfftf),rhog_dp(2,1:nfftf))
1654 [ - + ]: 18 : ABI_FREE(rhog_dp)
1655 : :
1656 : : ! Calculate the FFT index of each (G-Gp) vector and assign the value
1657 : : ! of the correspondent density simultaneously
1658 : 18 : ngfft1=ngfftf(1)
1659 : 18 : ngfft2=ngfftf(2)
1660 : 18 : ngfft3=ngfftf(3)
1661 : :
1662 [ - + + - : 72 : ABI_MALLOC_OR_DIE(mm, (npwc,npwc), ierr)
- + + - +
- + - + -
- + - + ]
1663 : :
1664 : 18 : ierr = 0
1665 [ + + ]: 504 : do ig=1,npwc
1666 [ + + ]: 13626 : do igp=1,npwc
1667 [ + + ]: 52488 : gmgp_idx = g2ifft(gvec(:,ig)-gvec(:,igp),ngfftf)
1668 [ + - ]: 13608 : if (gmgp_idx /= 0) then
1669 : 13122 : mm(ig,igp) = rhog(gmgp_idx)
1670 : : else
1671 : 0 : ierr = ierr + 1
1672 : 0 : mm(ig,igp) = czero
1673 : : end if
1674 : : end do
1675 : : end do
1676 : :
1677 [ - + ]: 18 : if (ierr /= 0) then
1678 : : write(msg,'(a,i0,3a)')&
1679 : 0 : 'Found ',ierr,' G1-G2 vectors falling outside the FFT box. ',ch10,&
1680 : 0 : 'Enlarge the FFT mesh to get rid of this problem. '
1681 : 0 : ABI_WARNING(msg)
1682 : : end if
1683 : :
1684 [ - + ]: 18 : ABI_FREE(rhog)
1685 : :
1686 : : ! Now we have rhogg, calculate the M matrix (q+G1).(q+G2) n(G1-G2)
1687 [ + + + + : 180 : b1=two_pi*gprimd(:,1); b2=two_pi*gprimd(:,2); b3=two_pi*gprimd(:,3)
+ + ]
1688 [ + + ]: 504 : do ig=1,npwc
1689 [ + + ]: 1944 : gpq(:)=gvec(:,ig)+qpt
1690 [ + + ]: 13626 : do igp=1,npwc
1691 [ + + ]: 52488 : gppq(:)=gvec(:,igp)+qpt
1692 : : qpg_dot_qpgp=zero
1693 [ + + ]: 52488 : do ii=1,3
1694 : : qpg_dot_qpgp = qpg_dot_qpgp + &
1695 : : ( gpq(1)*b1(ii) +gpq(2)*b2(ii) +gpq(3)*b3(ii))*&
1696 : 52488 : (gppq(1)*b1(ii)+gppq(2)*b2(ii)+gppq(3)*b3(ii))
1697 : : end do
1698 : 13608 : mm(ig,igp) = mm(ig,igp)*qpg_dot_qpgp
1699 : : end do ! igp
1700 : : end do ! ig
1701 : :
1702 : : ! Extract the reducible polarizability chi: e^{-1} = 1 + v chi
1703 : : ! \tilde\epsilon^{-1}_{G1 G2} = \delta_{G1 G2} + 4\pi \frac{\chi_{G1 G2}}{|q+G1| |q+G2|}
1704 : : !MG TODO too much memory in chi, we can do all this stuff inside a loop
1705 [ - + + - : 54 : ABI_MALLOC_OR_DIE(chi, (npwc,npwc), ierr)
- + + - +
- - + -
+ ]
1706 [ - + + - : 54 : ABI_MALLOC(qplusg, (npwc))
+ - + - -
+ - + ]
1707 : :
1708 [ + + + + ]: 13626 : chi(:,:)=epsm1(:,:)
1709 : 18 : qlist(:,1) = qpt
1710 : 18 : call cmod_qpg(1,1,qlist,npwc,gvec,gprimd,qplusg) !MG TODO here take care of small q
1711 : :
1712 [ + + ]: 504 : do ig=1,npwc
1713 : 504 : chi(ig,ig)=chi(ig,ig) - one
1714 : : end do
1715 : :
1716 [ + + ]: 504 : do ig=1,npwc
1717 [ + + ]: 13626 : do igp=1,npwc
1718 : 13608 : chi(ig,igp) = chi(ig,igp) * qplusg(ig) * qplusg(igp) / four_pi
1719 : : end do
1720 : : end do
1721 : :
1722 : : ! Solve chi(w=)*X = Lambda M*X where Lambda=-1/em(q)**2
1723 [ - + + - : 36 : ABI_MALLOC(eigval, (npwc))
- + - + ]
1724 [ - + + - : 54 : ABI_MALLOC_OR_DIE(mtemp, (npwc,npwc), ierr)
- + + - +
- - + -
+ ]
1725 : :
1726 : : ! Copy mm into working array as xhegv changes input matrices
1727 [ + + + + ]: 13626 : mtemp(:,:) = mm(:,:)
1728 : :
1729 : 18 : call xhegv(1,"Vectors","Upper",npwc,chi,mtemp,eigval)
1730 [ - + ]: 18 : ABI_FREE(mtemp)
1731 : :
1732 : : ! Now chi contains the eigenvectors.
1733 : : ! Eigenvectors are normalized as: X_i^* M X_j = \delta_{ij}
1734 : :
1735 : : ! Calculate the plasmon pole parameters
1736 : : ! good check: the lowest plasmon energy on gamma should be
1737 : : ! close to experimental plasma energy within an error of 10%
1738 : : ! this error can be reduced further if one includes the non local
1739 : : ! commutators in the calculation of the polarizability at q==0
1740 : :
1741 [ - + + - : 54 : ABI_MALLOC(tmp1,(npwc))
+ - + - -
+ - + ]
1742 [ - + + - : 54 : ABI_MALLOC_OR_DIE(zz2, (npwc, npwc), ierr)
- + + - +
- - + -
+ ]
1743 [ + + + + ]: 13626 : zz2(:,:)= zero
1744 : :
1745 : : ! Caller is responsible for handing small q case
1746 : 18 : qlist(:,1) = qpt
1747 : 18 : call cmod_qpg(1,1,qlist,npwc,gvec,gprimd,qplusg)
1748 : :
1749 [ + + ]: 504 : do ii=1,npwc
1750 : : ! keeping in mind that the above matrix is negative definite
1751 : : ! we might have a small problem with the eigvals corresponding to large G vectors
1752 : : ! i.e. DM band index, where the eigevalues become very small with
1753 : : ! possibility of being small positive numbers (due to numerical problems)
1754 : : ! thus as a caution one can use the following condition
1755 : : ! this will not affect the result since such a huge plasmon energy give almost zero
1756 : : ! contribution to the self-energy correlation energy.
1757 : :
1758 [ - + ]: 486 : if (eigval(ii)>=zero) then
1759 : : !write(msg,'(a,i0,a,es16.6)')' Imaginary plasmon pole eigenenergy, eigenvector number ',ii,' with eigval',eigval(ii),ch10
1760 : : !ABI_ERROR(msg)
1761 : 0 : eigval(ii) = -1.0d-4
1762 : : end if
1763 : :
1764 : : ! Save plasmon energies omega_p(q)
1765 : 486 : omegatw(ii) = SQRT(-one/eigval(ii))
1766 : :
1767 : : ! Calculate and save scaled plasmon-pole eigenvectors
1768 : : ! defined as \sqrt{4\pi} \frac{Mx}{\sqrt{\tilde\omega} |q+G|}
1769 [ + + ]: 13608 : tmp1(:)=chi(:,ii)
1770 : :
1771 [ + + ]: 13626 : do ig=1,npwc
1772 [ + + ]: 367416 : do igp=1,npwc
1773 : 367416 : zz2(ig,ii)=zz2(ig,ii)+mm(ig,igp)*tmp1(igp) ! z --> y
1774 : : end do
1775 : 13608 : bigomegatwsq(ig,ii)= SQRT(four_pi) * zz2(ig,ii) / SQRT(omegatw(ii)) / qplusg(ig)
1776 : : end do
1777 : :
1778 : : end do ! ii
1779 : :
1780 [ - + ]: 18 : ABI_FREE(tmp1)
1781 [ - + ]: 18 : ABI_FREE(eigval)
1782 [ - + ]: 18 : ABI_FREE(zz2)
1783 [ - + ]: 18 : ABI_FREE(qplusg)
1784 [ - + ]: 18 : ABI_FREE(chi)
1785 [ - + ]: 18 : ABI_FREE(mm)
1786 : :
1787 : 18 : bar = repeat('-', 80)
1788 : 18 : write(msg,'(3a)')bar,ch10,' plasmon energies in eV vs q vector shown for the lowest 10 bands'
1789 : 18 : call wrtout(std_out,msg)
1790 [ + - + + ]: 198 : write(msg,'(2x,5x,10f7.3)')(REAL(omegatw(ig))*Ha_eV, ig=1,min(10, npwc))
1791 : 18 : call wrtout(std_out,msg)
1792 : 18 : write(msg,'(a)')bar
1793 : 18 : call wrtout(std_out,msg)
1794 : :
1795 : 18 : write(msg,'(2a,f12.8,2a,3i5)')ch10,&
1796 [ + - - + : 522 : ' cppm4par: omega twiddle minval [eV] = ',MINVAL(ABS(omegatw))*Ha_eV,ch10,&
- - + + ]
1797 [ + - - + : 540 : ' omega twiddle min location = ',MINLOC(ABS(omegatw))
- - + + -
+ + + ]
1798 : 18 : call wrtout(std_out,msg)
1799 : :
1800 [ - + - + : 18 : end subroutine cppm4par
- + - + -
+ - + - +
- + - + -
+ - + ]
1801 : : !!***
1802 : :
1803 : : !----------------------------------------------------------------------
1804 : :
1805 : : !!****f* m_ppmodel/cqratio
1806 : : !! NAME
1807 : : !! cqratio
1808 : : !!
1809 : : !! FUNCTION
1810 : : !! Calculate qratio(G,Gp,q)= (q+G)\cdot(q+Gp) / |q+G|^2 needed for Hybertsen-Louie and Plasmonpole model
1811 : : !!
1812 : : !! INPUTS
1813 : : !! npwc=number of planewaves considered (used for the correlation part)
1814 : : !! gvec(3,npwc)=reduced coordinates of the plane waves
1815 : : !! q(3)=coordinates of q points
1816 : : !! gmet(3,3)=metric in reciprocal space
1817 : : !! gprimd(3,3)=reciprocal lattice vectors
1818 : : !!
1819 : : !! OUTPUT
1820 : : !! qratio(npwc,npwc)=(q+G).(q+Gp)
1821 : : !!
1822 : : !! SOURCE
1823 : :
1824 : 24 : subroutine cqratio(npwc, gvec, q, gmet, gprimd, qratio)
1825 : :
1826 : : !Arguments ------------------------------------
1827 : : !scalars
1828 : : integer,intent(in) :: npwc
1829 : : !arrays
1830 : : integer,intent(in) :: gvec(3,npwc)
1831 : : real(dp),intent(in) :: gmet(3,3),gprimd(3,3),q(3)
1832 : : real(dp),intent(out) :: qratio(npwc,npwc)
1833 : :
1834 : : !Local variables ------------------------------
1835 : : !scalars
1836 : : integer :: ig,igp,ii
1837 : : real(dp),parameter :: tol = 0.001_dp
1838 : : real(dp) :: qpg_dot_qpgp
1839 : : !arrays
1840 : 24 : real(dp) :: b1(3),b2(3),b3(3),gppq(3),gpq(3),norm(npwc)
1841 : : !************************************************************************
1842 : :
1843 [ + + + + : 240 : b1=two_pi*gprimd(:,1); b2=two_pi*gprimd(:,2); b3=two_pi*gprimd(:,3)
+ + ]
1844 : :
1845 [ + + + + : 21048 : norm(:)=zero; qratio=zero
+ + ]
1846 : :
1847 : : !FIXME this loops have to be rewritten!!!!
1848 [ + + ]: 708 : do ig=1,npwc
1849 [ + + ]: 2736 : gpq(:)=gvec(:,ig)+q
1850 [ + + + + : 10968 : norm(ig)=two_pi*SQRT(DOT_PRODUCT(gpq,MATMUL(gmet,gpq)))
+ + ]
1851 : : !norm(ig)=normv(gpq,gmet,'g')
1852 : : end do
1853 : :
1854 [ + + ]: 708 : do ig=1,npwc
1855 [ + + ]: 2736 : gpq(:)=gvec(:,ig)+q
1856 [ + + ]: 20364 : do igp=1,npwc
1857 [ + + ]: 78624 : gppq(:)=gvec(:,igp)+q
1858 : : qpg_dot_qpgp=zero
1859 : : !qpg_dot_qpgp=vdotw(gpq,gppq,gmet,'g')
1860 [ + + ]: 78624 : do ii=1,3
1861 : : qpg_dot_qpgp=qpg_dot_qpgp+&
1862 : : ( gpq(1)*b1(ii) + gpq(2)*b2(ii) + gpq(3)*b3(ii))*&
1863 : 78624 : (gppq(1)*b1(ii) + gppq(2)*b2(ii) +gppq(3)*b3(ii))
1864 : : end do
1865 : :
1866 : : ! Now calculate qratio = (q+G).(q+Gp)/|q+G|^2
1867 : : ! when |q+G|^2 and (q+G).(q+Gp) are both zero set (q+G).(q+Gp)/|q+G|^2 = 1
1868 : : ! when |q+G|^2 is zero and |q+Gp| is not zero set (q+G).(q+Gp)/|q+G|^2 = 0
1869 [ + + ]: 20340 : if (norm(ig) < tol) then
1870 [ + + ]: 114 : if (norm(igp) < tol) then ! Case q=0 and G=Gp=0
1871 : 4 : qratio(ig,igp) = one
1872 : : else ! Case q=0 and G=0 and Gp !=0
1873 : 110 : qratio(ig,igp) = zero
1874 : : end if
1875 [ + + ]: 19542 : else if (norm(igp) < tol) then ! Case q=0 and G= !0 and Gp=0
1876 : 110 : qratio(ig,igp)=zero
1877 : : else
1878 : 19432 : qratio(ig,igp)=qpg_dot_qpgp / norm(ig)**2
1879 : : end if
1880 : :
1881 : : end do
1882 : : end do
1883 : :
1884 : 24 : end subroutine cqratio
1885 : : !!***
1886 : :
1887 : : !----------------------------------------------------------------------
1888 : :
1889 : : !!****f* m_ppmodel/ppm_calc_sigc
1890 : : !!
1891 : : !! NAME
1892 : : !! ppm_calc_sigc
1893 : : !!
1894 : : !! FUNCTION
1895 : : !! Calculate the contribution to self-energy operator for a single band s in the band sum
1896 : : !! using a plasmon-pole model.
1897 : : !!
1898 : : !! INPUTS
1899 : : !! nspinor=Number of spinor components.
1900 : : !! npwc=Number of G vectors in the plasmon pole (correlation part)
1901 : : !! nomega=Number of frequencies.
1902 : : !! rhotwgp(npwx)=oscillator matrix elements divided by |q+G| i.e. $\frac{\langle b1 k-q s | e^{-i(q+G)r | b2 k s \rangle}{|q+G|}$.
1903 : : !! botsq(npwc,dm2_botsq)=Plasmon pole parameters for this q-point.
1904 : : !! otq(npwc,dm2_otq)=Plasmon pole parameters for this q-point.
1905 : : !! omegame0i(nomega)=($\omega$ - $\epsilon_i)$
1906 : : !! zcut=Small imaginary part to avoid the divergence. (see related input variable)
1907 : : !! theta_mu_minus_e0i= $\theta(\mu-\epsilon_{k-q,b1,s}), defines if the state is occupied or not.
1908 : : !! eig(dm_eig,dm_eig)=The eigvectors of the symmetrized inverse dielectric matrix for this q point
1909 : : !! (first index for G, second index for bands).
1910 : : !! npwx=number of G vectors in rhotwgp.
1911 : : !!
1912 : : !! OUTPUT
1913 : : !! ket(npwc,nomega):
1914 : : !!
1915 : : !! === model==1,2 ====
1916 : : !!
1917 : : !! ket(G,omega) += Sum_G2 Omega(G,G2) * rhotw(G2)
1918 : : !! ---------------------------------------------------
1919 : : !! 2 omegatw(G,G2) (omega - E_i + omegatw(G,G2)(2f-1))
1920 : : !!
1921 : : !! sigcme(nomega) (to be described), only relevant if ppm3 or ppm4
1922 : : !!
1923 : : !! NOTES
1924 : : !! The i/two_pi factor in the convolution between G and W is included in this routine.
1925 : : !!
1926 : : !! TODO:
1927 : : !! Use BLAS for better efficiency
1928 : : !!
1929 : : !! SOURCE
1930 : :
1931 : 1558650 : subroutine ppm_calc_sigc(ppm, nspinor, npwc, nomega, rhotwgp, botsq, otq, &
1932 : 1558650 : omegame0i, zcut, theta_mu_minus_e0i, eig, npwx, ket, sigcme)
1933 : :
1934 : : !Arguments ------------------------------------
1935 : : !scalars
1936 : : class(ppmodel_t),intent(in) :: ppm
1937 : : integer,intent(in) :: nomega, npwc, npwx, nspinor
1938 : : real(dp),intent(in) :: theta_mu_minus_e0i, zcut
1939 : : !arrays
1940 : : real(dp),intent(in) :: omegame0i(nomega)
1941 : : complex(gwp),intent(in) :: botsq(npwc, ppm%dm2_botsq), eig(ppm%dm_eig, ppm%dm_eig), otq(npwc, ppm%dm2_otq)
1942 : : complex(gwp),intent(in) :: rhotwgp(npwx, nspinor)
1943 : : complex(gwp),intent(inout) :: ket(npwc, nspinor, nomega)
1944 : : complex(gwp),intent(out) :: sigcme(nomega)
1945 : :
1946 : : !Local variables-------------------------------
1947 : : !scalars
1948 : : integer :: ig,igp,ii,iw,ispinor
1949 : : real(dp),parameter :: tol_occ = tol3, tol_omega = tol6
1950 : : real(dp) :: den, den2, ff, inv_den, omegame0i_io, otw, twofm1, twofm1_zcut, twofm1_zcut2, zcut2
1951 : : complex(gwp) :: ct, num, numf, rhotwgdp_igp
1952 : : logical :: fully_occupied, totally_empty
1953 : : !character(len=500) :: msg
1954 : : !arrays
1955 : 1558650 : complex(gwp),allocatable :: rhotwgdpcc(:)
1956 : : !*************************************************************************
1957 : :
1958 : 1558650 : zcut2 = zcut ** 2
1959 : :
1960 : 3113460 : select case (ppm%model)
1961 : :
1962 : : case (PPM_GODBY_NEEDS, PPM_HYBERTSEN_LOUIE)
1963 : 1554810 : fully_occupied = (abs(theta_mu_minus_e0i-one) < tol_occ)
1964 : 1554810 : totally_empty = (abs(theta_mu_minus_e0i ) < tol_occ)
1965 : :
1966 [ + + ]: 3116340 : do ispinor=1,nspinor
1967 : :
1968 [ + + ]: 1561530 : if (.not. totally_empty) then
1969 : : ! \Bomega^2_{G1G2}/\omegat_{G1G2} M_{G1,G2}. \theta(\mu-e_s) / (\omega+\omegat_{G1G2}-e_s-i\delta)
1970 : 4436630 : twofm1_zcut = zcut
1971 : 4436630 : twofm1_zcut2 = zcut2
1972 : : !$omp parallel do private(omegame0i_io, rhotwgdp_igp, otw, num, den, den2)
1973 [ + + ]: 4436630 : do iw=1,nomega
1974 : 3915487 : omegame0i_io = omegame0i(iw)
1975 : : !if (iw > 1 .and. abs(omegame0i(iw) - omegame0i(iw-1)) < tol_omega) then
1976 : : ! ket(:,:,iw) = ket(:,:,iw-1); cycle
1977 : : !end if
1978 : :
1979 [ + + ]: 96965859 : do igp=1,npwc
1980 : 92529229 : rhotwgdp_igp = rhotwgp(igp, ispinor)
1981 [ + + ]: 4833905347 : do ig=1,npwc
1982 : 4737460631 : otw = DBLE(otq(ig,igp)) !in principle otw -> otw - ieta
1983 : 4737460631 : num = botsq(ig,igp) * rhotwgdp_igp
1984 : 4737460631 : den = omegame0i_io + otw
1985 : 4737460631 : den2 = den ** 2
1986 [ + + ]: 4829989860 : if (den2 > zcut2) then
1987 : 4734382654 : ket(ig,ispinor, iw) = ket(ig,ispinor,iw) + num/(den*otw) * theta_mu_minus_e0i
1988 : : else
1989 : : ket(ig,ispinor,iw) = ket(ig,ispinor,iw) + &
1990 : 3077977 : num * CMPLX(den,twofm1_zcut) / ((den2 + twofm1_zcut2) * otw) * theta_mu_minus_e0i
1991 : : end if
1992 : : end do ! ig
1993 : : end do ! igp
1994 : : end do ! iw
1995 : : end if ! not totally empty
1996 : :
1997 [ + + ]: 3116340 : if (.not. fully_occupied) then
1998 : : ! \Bomega^2_{G1G2}/\omegat_{G1G2} M_{G1,G2}. \theta(e_s-\mu) / (\omega-\omegat_{G1G2}-e_s+i\delta)
1999 : 1041183 : twofm1_zcut = -zcut
2000 : 1041183 : twofm1_zcut2 = twofm1_zcut**2
2001 : : !$omp parallel do private(omegame0i_io, rhotwgdp_igp, otw, num, den, den2)
2002 [ + + ]: 18225046 : do iw=1,nomega
2003 : 17183863 : omegame0i_io = omegame0i(iw)
2004 : : !if (iw > 1 .and. abs(omegame0i(iw) - omegame0i(iw-1)) < tol_omega) then
2005 : : ! ket(:,:,iw) = ket(:,:,iw-1); cycle
2006 : : !end if
2007 : :
2008 [ + + ]: 464817479 : do igp=1,npwc
2009 : 446592433 : rhotwgdp_igp = rhotwgp(igp, ispinor)
2010 [ + + ]: 31046850919 : do ig=1,npwc
2011 : 30583074623 : otw = DBLE(otq(ig,igp)) !in principle otw -> otw + ieta
2012 : 30583074623 : num = botsq(ig,igp) * rhotwgdp_igp
2013 : 30583074623 : den = omegame0i_io - otw
2014 : 30583074623 : den2 = den ** 2
2015 [ + + ]: 31029667056 : if (den2 > zcut2) then
2016 : 30582157413 : ket(ig,ispinor,iw) = ket(ig,ispinor, iw) + num / (den*otw) * (one-theta_mu_minus_e0i)
2017 : : else
2018 : : ket(ig,ispinor, iw) = ket(ig,ispinor,iw) + &
2019 : 917210 : num * CMPLX(den,twofm1_zcut) / ((den2 + twofm1_zcut2) * otw) * (one-theta_mu_minus_e0i)
2020 : : end if
2021 : : end do ! ig
2022 : : end do ! igp
2023 : : end do ! iw
2024 : : end if ! not fully occupied
2025 : :
2026 : : end do ! ispinor
2027 : :
2028 [ + + + + : 582613416 : ket=ket*half
+ + ]
2029 : :
2030 : : case (PPM_LINDEN_HORSH, PPM_ENGEL_FARID)
2031 [ - + ]: 3840 : ABI_CHECK(nspinor == 1, "nspinor/=1 not allowed")
2032 : :
2033 : : ! rho-twiddle(G) is formed, introduce rhotwgdpcc, for speed reason
2034 [ - + + - : 11520 : ABI_MALLOC(rhotwgdpcc, (npwx))
+ - + - -
+ - + ]
2035 : :
2036 : 3840 : ff = theta_mu_minus_e0i ! occupation number f (include poles if ...)
2037 : 3840 : twofm1 = two*ff-one ! 2f-1
2038 : 3840 : twofm1_zcut = twofm1*zcut
2039 [ + + ]: 107520 : rhotwgdpcc(:) = CONJG(rhotwgp(:, 1))
2040 : :
2041 [ + + ]: 33280 : do iw=1,nomega
2042 : 29440 : omegame0i_io = omegame0i(iw)
2043 : 29440 : ct = czero_gw
2044 [ + + ]: 824320 : do ii=1,npwc ! Loop over the DM bands
2045 : 794880 : num = czero_gw
2046 : :
2047 : 794880 : select case (ppm%model)
2048 : : case (PPM_LINDEN_HORSH)
2049 : : ! Calculate \beta (eq. 106 pag 47)
2050 [ + + ]: 11128320 : do ig=1,npwc
2051 : 11128320 : num = num + rhotwgdpcc(ig)*eig(ig,ii)
2052 : : end do
2053 : 397440 : numf=num*CONJG(num) !MG this means that we cannot do SCGW
2054 : 397440 : numf=numf*botsq(ii,1)
2055 : :
2056 : : case (PPM_ENGEL_FARID)
2057 [ + + ]: 11128320 : do ig=1,npwc
2058 : 11128320 : num = num + rhotwgdpcc(ig)*botsq(ig,ii)
2059 : : end do
2060 : 397440 : numf = num*CONJG(num) !MG this means that we cannot do SCGW
2061 : :
2062 : : case default
2063 [ + + - ]: 794880 : ABI_ERROR("Wrong ppm%model")
2064 : : end select
2065 : :
2066 : 794880 : otw=DBLE(otq(ii,1)) ! in principle otw -> otw - ieta
2067 : 794880 : den=omegame0i_io+otw*twofm1
2068 : :
2069 [ + - ]: 824320 : if (den**2 > zcut**2) then
2070 : 794880 : inv_den=one/den
2071 : 794880 : ct=ct+numf*inv_den
2072 : : else
2073 : 0 : inv_den = one/((den**2+twofm1_zcut**2))
2074 : 0 : ct = ct + numf*CMPLX(den,twofm1_zcut)*inv_den
2075 : : end if
2076 : :
2077 : : end do ! ii DM bands
2078 : 33280 : sigcme(iw) = ct*half
2079 : :
2080 : : !if (ppm%model == PPM_ENGEL_FARID) then
2081 : : !ct = dot_product(ket(:, iw), ket(:, iw))
2082 : : !if (abs(sigcme(iw) - ct) > tol12) then
2083 : : ! ABI_ERROR("foo bar")
2084 : : !end if
2085 : : !end if
2086 : : end do ! iw
2087 : :
2088 [ - + ]: 3840 : ABI_FREE(rhotwgdpcc)
2089 : :
2090 : : case default
2091 [ + + - ]: 1558650 : ABI_BUG(sjoin('Wrong ppm%model:',itoa(ppm%model)))
2092 : : end select
2093 : :
2094 : 1558650 : end subroutine ppm_calc_sigc
2095 : : !!***
2096 : :
2097 : : !----------------------------------------------------------------------
2098 : :
2099 : : !!****f* m_ppmodel/ppm_calc_dsigc
2100 : : !! NAME
2101 : : !! ppm_calc_dsigc
2102 : : !!
2103 : : !! FUNCTION
2104 : : !! Compute the frequency derivative of the plasmon-pole convolution performed in ppm_calc_sigc,
2105 : : !! i.e. the same contraction but with the simple pole 1/(omega - E_i +- omegatw) replaced by
2106 : : !! its omega-derivative -1/(omega - E_i +- omegatw)**2.
2107 : : !!
2108 : : !! This kernel is needed for the additional GWPT term at q=Gamma arising from the first-order
2109 : : !! derivative of the KS eigenvalue in Delta G0 (double pole, see Appendix A of the GWPT paper):
2110 : : !!
2111 : : !! Delta_{kappa,alpha,Gamma} G0_{nk}(e) = - (d eps_nk / d tau_{kappa,alpha}) * dG0_{nk}(e)/de
2112 : : !!
2113 : : !! The caller (gwpt_run) uses the output ket(G,omega) as
2114 : : !!
2115 : : !! g^Sigma_{mn}(omega) += - Sum_G <m,k| e^{i(p+G)r} |chi_{n'k-p}> vc_sqrt(p+G) * ket(G,omega)
2116 : : !!
2117 : : !!
2118 : : !! INPUTS
2119 : : !! Same conventions as ppm_calc_sigc (GN and HL plasmon-pole models only).
2120 : : !!
2121 : : !! OUTPUT
2122 : : !! ket(npwc,nspinor,nomega): incremented with
2123 : : !!
2124 : : !! ket(G,omega) += - Sum_G2 Omega(G,G2) * rhotw(G2) * [2f-1 pole structure]
2125 : : !! --------------------------------------------------
2126 : : !! 2 omegatw(G,G2) (omega - E_i + omegatw(G,G2)(2f-1))**2
2127 : : !!
2128 : : !! NOTES
2129 : : !! The i/two_pi factor in the convolution between G and W is included, as in ppm_calc_sigc.
2130 : : !! The zcut regularization is applied to the complex simple pole before squaring so that
2131 : : !! the derivative kernel is consistent with the regularized parent kernel.
2132 : : !!
2133 : : !! SOURCE
2134 : :
2135 : 0 : subroutine ppm_calc_dsigc(ppm, nspinor, npwc, nomega, rhotwgp, botsq, otq, &
2136 : 0 : omegame0i, zcut, theta_mu_minus_e0i, npwx, ket)
2137 : :
2138 : : !Arguments ------------------------------------
2139 : : !scalars
2140 : : class(ppmodel_t),intent(in) :: ppm
2141 : : integer,intent(in) :: nomega, npwc, npwx, nspinor
2142 : : real(dp),intent(in) :: theta_mu_minus_e0i, zcut
2143 : : !arrays
2144 : : real(dp),intent(in) :: omegame0i(nomega)
2145 : : complex(gwp),intent(in) :: botsq(npwc, ppm%dm2_botsq), otq(npwc, ppm%dm2_otq)
2146 : : complex(gwp),intent(in) :: rhotwgp(npwx, nspinor)
2147 : : complex(gwp),intent(inout) :: ket(npwc, nspinor, nomega)
2148 : :
2149 : : !Local variables-------------------------------
2150 : : !scalars
2151 : : integer :: ig,igp,iw,ispinor
2152 : : real(dp),parameter :: tol_occ = tol3
2153 : : real(dp) :: den, den2, omegame0i_io, otw, twofm1_zcut, twofm1_zcut2, zcut2
2154 : : complex(gwp) :: num, inv_cden, rhotwgdp_igp
2155 : : logical :: fully_occupied, totally_empty
2156 : : !*************************************************************************
2157 : :
2158 : 0 : zcut2 = zcut ** 2
2159 : :
2160 : 0 : select case (ppm%model)
2161 : :
2162 : : case (PPM_GODBY_NEEDS, PPM_HYBERTSEN_LOUIE)
2163 : 0 : fully_occupied = (abs(theta_mu_minus_e0i-one) < tol_occ)
2164 : 0 : totally_empty = (abs(theta_mu_minus_e0i ) < tol_occ)
2165 : :
2166 [ # # ]: 0 : do ispinor=1,nspinor
2167 : :
2168 [ # # ]: 0 : if (.not. totally_empty) then
2169 : : ! d/domega of \Bomega^2_{G1G2}/\omegat_{G1G2} M_{G1,G2} \theta(\mu-e_s) / (\omega+\omegat_{G1G2}-e_s-i\delta)
2170 : 0 : twofm1_zcut = zcut
2171 : 0 : twofm1_zcut2 = zcut2
2172 : : !$omp parallel do private(omegame0i_io, rhotwgdp_igp, otw, num, den, den2, inv_cden)
2173 [ # # ]: 0 : do iw=1,nomega
2174 : 0 : omegame0i_io = omegame0i(iw)
2175 [ # # ]: 0 : do igp=1,npwc
2176 : 0 : rhotwgdp_igp = rhotwgp(igp, ispinor)
2177 [ # # ]: 0 : do ig=1,npwc
2178 : 0 : otw = DBLE(otq(ig,igp)) !in principle otw -> otw - ieta
2179 : 0 : num = botsq(ig,igp) * rhotwgdp_igp
2180 : 0 : den = omegame0i_io + otw
2181 : 0 : den2 = den ** 2
2182 [ # # ]: 0 : if (den2 > zcut2) then
2183 : 0 : ket(ig,ispinor,iw) = ket(ig,ispinor,iw) - num/(den2*otw) * theta_mu_minus_e0i
2184 : : else
2185 : 0 : inv_cden = CMPLX(den,twofm1_zcut) / (den2 + twofm1_zcut2)
2186 : : ket(ig,ispinor,iw) = ket(ig,ispinor,iw) - &
2187 : 0 : num * inv_cden**2 / otw * theta_mu_minus_e0i
2188 : : end if
2189 : : end do ! ig
2190 : : end do ! igp
2191 : : end do ! iw
2192 : : end if ! not totally empty
2193 : :
2194 [ # # ]: 0 : if (.not. fully_occupied) then
2195 : : ! d/domega of \Bomega^2_{G1G2}/\omegat_{G1G2} M_{G1,G2} \theta(e_s-\mu) / (\omega-\omegat_{G1G2}-e_s+i\delta)
2196 : 0 : twofm1_zcut = -zcut
2197 : 0 : twofm1_zcut2 = twofm1_zcut**2
2198 : : !$omp parallel do private(omegame0i_io, rhotwgdp_igp, otw, num, den, den2, inv_cden)
2199 [ # # ]: 0 : do iw=1,nomega
2200 : 0 : omegame0i_io = omegame0i(iw)
2201 [ # # ]: 0 : do igp=1,npwc
2202 : 0 : rhotwgdp_igp = rhotwgp(igp, ispinor)
2203 [ # # ]: 0 : do ig=1,npwc
2204 : 0 : otw = DBLE(otq(ig,igp)) !in principle otw -> otw + ieta
2205 : 0 : num = botsq(ig,igp) * rhotwgdp_igp
2206 : 0 : den = omegame0i_io - otw
2207 : 0 : den2 = den ** 2
2208 [ # # ]: 0 : if (den2 > zcut2) then
2209 : 0 : ket(ig,ispinor,iw) = ket(ig,ispinor,iw) - num/(den2*otw) * (one-theta_mu_minus_e0i)
2210 : : else
2211 : 0 : inv_cden = CMPLX(den,twofm1_zcut) / (den2 + twofm1_zcut2)
2212 : : ket(ig,ispinor,iw) = ket(ig,ispinor,iw) - &
2213 : 0 : num * inv_cden**2 / otw * (one-theta_mu_minus_e0i)
2214 : : end if
2215 : : end do ! ig
2216 : : end do ! igp
2217 : : end do ! iw
2218 : : end if ! not fully occupied
2219 : :
2220 : : end do ! ispinor
2221 : :
2222 [ # # # # : 0 : ket = ket * half
# # ]
2223 : :
2224 : : case default
2225 [ # # ]: 0 : ABI_ERROR(sjoin('ppm_calc_dsigc not implemented for ppm%model:', itoa(ppm%model)))
2226 : : end select
2227 : :
2228 : 0 : end subroutine ppm_calc_dsigc
2229 : : !!***
2230 : :
2231 : : !----------------------------------------------------------------------
2232 : :
2233 : : !!****f* m_ppmodel/ppm_rotate_iqbz
2234 : : !! NAME
2235 : : !! ppm_rotate_iqbz
2236 : : !!
2237 : : !! FUNCTION
2238 : : !! Symmetrize the plasmonpole parameters in the full BZ.
2239 : : !!
2240 : : !! INPUTS
2241 : : !! iq_bz=Index of the q-point in the BZ where the ppmodel parameters are wanted.
2242 : : !! Gsph<gsphere_t>=data related to the G-sphere.
2243 : : !! Cryst<crystal_t>=Info on the unit cell and crystal symmetries.
2244 : : !! Qmesh<kmesh_t>=the q-mesh used for the inverse dielectric matrix
2245 : : !! iq_ibz=Index of the q-point in the BZ.
2246 : : !! npwe=number of G vectors for the correlation part
2247 : : !! nomega=number of frequencies in $\epsilon^{-1}$
2248 : : !! omega=frequencies in epsm1_ggw
2249 : : !! epsm1_ggw(npwe,npwe,nomega)=the inverse dielctric matrix
2250 : : !! ngfftf(18)=contain all needed information about the 3D fine FFT mesh, see ~abinit/doc/variables/vargs.htm#ngfft
2251 : : !! gprimd(3,3)=dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$)
2252 : : !! nfftf=the number of points in the FFT mesh (for this processor)
2253 : : !! rhor_tot(nfftf)=the total charge in real space
2254 : : !!
2255 : : !! SIDE EFFECTS
2256 : : !! ppm<ppmodel_t>=data type containing information on the plasmonpole technique.
2257 : : !! Internal tables are modified so that they (point|store) the plasmon-pole parameters
2258 : : !! for the specified q-point in the BZ.
2259 : : !!
2260 : : !! SOURCE
2261 : :
2262 : 0 : subroutine ppm_rotate_iqbz(ppm, iq_bz, Cryst, Qmesh, Gsph, npwe, nomega, omega, epsm1_ggw, &
2263 : 0 : nfftf, ngfftf, rhor_tot)
2264 : :
2265 : : !Arguments ------------------------------------
2266 : : !scalars
2267 : : class(ppmodel_t),target,intent(inout) :: ppm
2268 : : integer,intent(in) :: nfftf,npwe,nomega,iq_bz
2269 : : type(crystal_t),intent(in) :: Cryst
2270 : : type(gsphere_t),intent(in) :: Gsph
2271 : : type(kmesh_t),intent(in) :: Qmesh
2272 : : !arrays
2273 : : integer,intent(in) :: ngfftf(18)
2274 : : real(dp),intent(in) :: rhor_tot(nfftf)
2275 : : complex(dp),intent(in) :: omega(nomega)
2276 : : complex(gwp),intent(in) :: epsm1_ggw(npwe,npwe,nomega)
2277 : :
2278 : : !Local variables-------------------------------
2279 : : !scalars
2280 : : integer :: iq_ibz,itim_q,isym_q,iq_curr
2281 : : logical :: q_isirred
2282 : : !character(len=500) :: msg
2283 : : !arrays
2284 : : real(dp) :: qbz(3)
2285 : : ! *********************************************************************
2286 : :
2287 : : ! Save the index of the q-point in the BZ for checking purpose.
2288 : 0 : ppm%iq_bz = iq_bz
2289 : :
2290 : 0 : call qmesh%get_bz_item(iq_bz, qbz, iq_ibz, isym_q, itim_q, isirred=q_isirred)
2291 : 0 : iq_curr = iq_ibz; if (ppm%mqmem == 0) iq_curr = 1
2292 : :
2293 : : ! =======================================================
2294 : : ! ==== Branching for in-core or out-of-core solution ====
2295 : : ! =======================================================
2296 : :
2297 : : ! Allocate the tables for this q_ibz
2298 : : !print *, "ppm%has_qibz(iq_ibz)", ppm%has_qibz(iq_ibz), "q_isirred:", q_isirred
2299 [ # # ]: 0 : if (ppm%has_qibz(iq_ibz) == PPM_NOTAB) call ppm%malloc_iqibz(iq_ibz)
2300 : :
2301 [ # # ]: 0 : if (ppm%has_qibz(iq_ibz) == PPM_TAB_ALLOCATED) then
2302 : : ! Calculate the ppmodel tables for this q_ibz
2303 : 0 : call ppm%new_setup(iq_ibz, Cryst, Qmesh, npwe, nomega, omega, epsm1_ggw, nfftf, Gsph%gvec, ngfftf, rhor_tot)
2304 : : end if
2305 : :
2306 : : ! Allocate memory if not done yet.
2307 : : #ifdef FC_LLVM
2308 : : !FIXME I don't understand why LLVM fails here...
2309 : : !I put preproc so others know extra spaces are on purpose
2310 : : ABI_REMALLOC(ppm%bigomegatwsq_qbz_vals, (ppm%npwc, ppm%dm2_botsq) )
2311 : : ABI_REMALLOC(ppm%omegatw_qbz_vals, (ppm%npwc, ppm%dm2_otq) )
2312 : : ABI_REMALLOC(ppm%eigpot_qbz_vals, (ppm%dm_eig, ppm%dm_eig) )
2313 : : #else
2314 [ # # # # : 0 : ABI_REMALLOC(ppm%bigomegatwsq_qbz_vals, (ppm%npwc, ppm%dm2_botsq))
# # # # #
# # # # #
# # # # #
# ]
2315 [ # # # # : 0 : ABI_REMALLOC(ppm%omegatw_qbz_vals, (ppm%npwc, ppm%dm2_otq))
# # # # #
# # # # #
# # # # #
# ]
2316 [ # # # # : 0 : ABI_REMALLOC(ppm%eigpot_qbz_vals, (ppm%dm_eig, ppm%dm_eig))
# # # # #
# # # # #
# # # # #
# ]
2317 : : #endif
2318 : :
2319 [ # # ]: 0 : if (q_isirred) then
2320 : : ! Symmetrization is not needed. Copy the data in memory and change the status.
2321 [ # # # # : 0 : ppm%bigomegatwsq_qbz_vals = ppm%bigomegatwsq(iq_ibz)%vals
# # # # #
# # # # #
# # # # #
# ]
2322 [ # # # # : 0 : ppm%omegatw_qbz_vals = ppm%omegatw(iq_ibz)%vals
# # # # #
# # # # #
# # # # #
# ]
2323 [ # # # # : 0 : ppm%eigpot_qbz_vals = ppm%eigpot(iq_ibz)%vals
# # # # #
# # # # #
# # # # #
# ]
2324 : :
2325 : : else
2326 : : ! q-point in the BZ. Calculate new table for this q-point in the BZ. Beware: Dimensions should not change.
2327 [ # # # # : 0 : call ppm%get_qbz(Gsph, Qmesh, iq_bz, ppm%bigomegatwsq_qbz_vals, ppm%omegatw_qbz_vals, ppm%eigpot_qbz_vals)
# # ]
2328 : :
2329 : : ! Release the table in the IBZ if required.
2330 [ # # ]: 0 : if (.not. ppm%keep_qibz(iq_ibz)) call ppm%table_free_iqibz(iq_ibz)
2331 : : end if
2332 : :
2333 : 0 : end subroutine ppm_rotate_iqbz
2334 : : !!***
2335 : :
2336 : : !----------------------------------------------------------------------
2337 : :
2338 : : !!****f* m_ppmodel/ppm_new_setup
2339 : : !! NAME
2340 : : !! ppm_new_setup
2341 : : !!
2342 : : !! FUNCTION
2343 : : !! Initialize some values of several arrays of the ppm datastructure
2344 : : !! that are used in case of plasmonpole calculations
2345 : : !! Just a wrapper around different plasmonpole routines.
2346 : : !!
2347 : : !! INPUTS
2348 : : !! iq_ibz=Index of the q-point in the BZ.
2349 : : !! Cryst<crystal_t>=Info on the unit cell and crystal symmetries.
2350 : : !! Qmesh<kmesh_t>=the q-mesh used for the inverse dielectric matrix
2351 : : !! npwe=number of G vectors for the correlation part
2352 : : !! nomega=number of frequencies in $\epsilon^{-1}$
2353 : : !! omega=frequencies in epsm1_ggw
2354 : : !! epsm1_ggw(npwe,npwe,nomega)=the inverse dielctric matrix
2355 : : !! nfftf=the number of points in the FFT mesh (for this processor)
2356 : : !! ngfftf(18)=contain all needed information about the 3D fine FFT mesh, see ~abinit/doc/variables/vargs.htm#ngfft
2357 : : !! rhor(nfftf)=the total charge in real space.
2358 : : !!
2359 : : !! SIDE EFFECTS
2360 : : !! == if ppmodel 1 or 2 ==
2361 : : !! %omegatw and %bigomegatwsq
2362 : : !! == if ppmodel 3 ==
2363 : : !! %omegatw, %bigomegatwsq and %eigpot
2364 : : !! == if ppmodel 4 ==
2365 : : !! %omegatw and %bigomegatwsq
2366 : : !!
2367 : : !! NOTES
2368 : : !! * FFT parallelism not implemented.
2369 : : !! * TODO: rhor_tot should be replaced by rhog_tot to avoid nq_ibz FFTs.
2370 : : !!
2371 : : !! SOURCE
2372 : :
2373 : 19 : subroutine ppm_new_setup(ppm, iq_ibz, Cryst, Qmesh, npwe, nomega, omega, epsm1_ggw, nfftf, gvec, ngfftf, rhor_tot)
2374 : :
2375 : : !Arguments ------------------------------------
2376 : : !scalars
2377 : : class(ppmodel_t),intent(inout) :: ppm
2378 : : integer,intent(in) :: nfftf,npwe,nomega,iq_ibz
2379 : : type(kmesh_t),intent(in) :: Qmesh
2380 : : type(crystal_t),intent(in) :: Cryst
2381 : : !arrays
2382 : : integer,intent(in) :: gvec(3,npwe),ngfftf(18)
2383 : : real(dp),intent(in) :: rhor_tot(nfftf)
2384 : : complex(dp),intent(in) :: omega(nomega)
2385 : : complex(gwp),intent(in) :: epsm1_ggw(npwe,npwe,nomega)
2386 : :
2387 : : !Local variables-------------------------------
2388 : : !scalars
2389 : : real(dp) :: n_at_G_zero
2390 : : character(len=500) :: msg
2391 : : !scalars
2392 : : real(dp) :: qpt(3)
2393 : : ! *************************************************************************
2394 : :
2395 [ - + ]: 19 : if (ppm%has_qibz(iq_ibz) /= PPM_TAB_ALLOCATED) then
2396 : 0 : ABI_ERROR(sjoin("ppmodel tables for iq_ibz:", itoa(iq_ibz), "are not allocated! has_qibz=", itoa(ppm%has_qibz(iq_ibz))))
2397 : : end if
2398 : :
2399 [ + + ]: 76 : qpt = Qmesh%ibz(:,iq_ibz)
2400 : 19 : ppm%has_qibz(iq_ibz) = PPM_TAB_STORED
2401 : :
2402 : : ! Calculate plasmonpole parameters
2403 : 19 : select case (ppm%model)
2404 : :
2405 : : case (PPM_NONE)
2406 : 0 : ABI_COMMENT('Skipping plasmonpole model calculation')
2407 : :
2408 : : case (PPM_GODBY_NEEDS)
2409 : : ! Note: the q-dependence enters only through epsilon^-1.
2410 : 19 : call cppm1par(npwe, nomega, omega, ppm%drude_plsmf, epsm1_ggw, ppm%omegatw(iq_ibz)%vals, ppm%bigomegatwsq(iq_ibz)%vals)
2411 : :
2412 : : case (PPM_HYBERTSEN_LOUIE)
2413 : : call cppm2par(qpt, npwe, epsm1_ggw(:,:,1), ngfftf, gvec, Cryst%gprimd, rhor_tot, nfftf, Cryst%gmet, &
2414 : 0 : ppm%bigomegatwsq(iq_ibz)%vals, ppm%omegatw(iq_ibz)%vals, ppm%invalid_freq)
2415 : :
2416 : : ! Quick-and-dirty change of the plasmon frequency. Never executed in standard runs.
2417 [ # # ]: 0 : if (ppm%force_plsmf > tol6) then
2418 : : ! Integrate the real-space density
2419 [ # # ]: 0 : n_at_G_zero = SUM(rhor_tot(:))/nfftf
2420 : : ! Change the prefactor
2421 : 0 : write(msg,'(2(a,es16.8))') 'Forced ppmfreq: ',ppm%force_plsmf*Ha_eV,' nelect/ucvol: ',n_at_G_zero
2422 : 0 : ABI_WARNING(msg)
2423 : :
2424 : 0 : ppm%force_plsmf = (ppm%force_plsmf**2)/(four_pi*n_at_G_zero)
2425 [ # # # # ]: 0 : ppm%bigomegatwsq(iq_ibz)%vals = ppm%force_plsmf * ppm%bigomegatwsq(iq_ibz)%vals
2426 [ # # # # ]: 0 : ppm%omegatw(iq_ibz)%vals = ppm%force_plsmf * ppm%omegatw(iq_ibz)%vals
2427 : 0 : write(msg,'(a,es16.8)') 'Plasma frequency forced in HL ppmodel, new prefactor is: ',ppm%force_plsmf
2428 : 0 : ABI_WARNING(msg)
2429 : : end if
2430 : :
2431 : : case (PPM_LINDEN_HORSH)
2432 : : call cppm3par(qpt, npwe,epsm1_ggw(:,:,1), ngfftf,gvec, Cryst%gprimd, rhor_tot, nfftf, &
2433 : 0 : ppm%bigomegatwsq(iq_ibz)%vals, ppm%omegatw(iq_ibz)%vals(:,1), ppm%eigpot(iq_ibz)%vals)
2434 : :
2435 : : case (PPM_ENGEL_FARID)
2436 [ # # # # : 0 : if ((ALL(ABS(qpt)<1.0e-3))) qpt = GW_Q0_DEFAULT ! FIXME
# # ]
2437 : :
2438 : : call cppm4par(qpt, npwe,epsm1_ggw(:,:,1), ngfftf, gvec, Cryst%gprimd, rhor_tot, nfftf, &
2439 : 0 : ppm%bigomegatwsq(iq_ibz)%vals, ppm%omegatw(iq_ibz)%vals(:,1))
2440 : :
2441 : : case default
2442 [ - + - - : 19 : ABI_BUG(sjoin('Wrong ppm%model:', itoa(ppm%model)))
- - ]
2443 : : end select
2444 : :
2445 : 19 : end subroutine ppm_new_setup
2446 : : !!***
2447 : :
2448 : : !!****f* m_ppmodel/ppm_print
2449 : : !! NAME
2450 : : !! ppm_print
2451 : : !!
2452 : : !! FUNCTION
2453 : : !! Print info on object
2454 : : !!
2455 : : !! SOURCE
2456 : :
2457 : 6 : subroutine ppm_print(ppm, units, header)
2458 : :
2459 : : !Arguments ------------------------------------
2460 : : class(ppmodel_t),intent(in) :: ppm
2461 : : integer,intent(in) :: units(:)
2462 : : character(len=*),optional,intent(in) :: header
2463 : :
2464 : : !Local variables-------------------------------
2465 : : character(len=500) :: msg
2466 : : type(yamldoc_t) :: ydoc
2467 : : !*************************************************************************
2468 : :
2469 [ # # # # ]: 0 : msg = ' ==== Info on the ppm_t object ==== '; if (present(header)) msg=' ==== '//trim(adjustl(header))//' ==== '
2470 : 3 : call wrtout(units, msg)
2471 : :
2472 : 3 : ydoc = yamldoc_open('Plasmonpole_params') !, width=11, real_fmt='(3f8.3)')
2473 : : !call ydoc%add_string("gwr_task", )
2474 : 3 : call ydoc%add_int("dm2_botsq", ppm%dm2_botsq)
2475 : 3 : call ydoc%add_int("dm_eig", ppm%dm_eig)
2476 : 3 : call ydoc%add_int("dm2_otq", ppm%dm2_otq)
2477 : 3 : call ydoc%add_int("invalid_freq", ppm%invalid_freq)
2478 : 3 : call ydoc%add_int("model", ppm%model)
2479 : 3 : call ydoc%add_int("mqmem", ppm%mqmem)
2480 : 3 : call ydoc%add_int("nqibz", ppm%nqibz)
2481 : 3 : call ydoc%add_int("npwc", ppm%npwc)
2482 : 3 : call ydoc%add_int("userho", ppm%userho)
2483 : 3 : call ydoc%add_int("iq_bz", ppm%iq_bz)
2484 : 3 : call ydoc%add_real("drude_plsmf", ppm%drude_plsmf)
2485 : 3 : call ydoc%add_real("force_plsmf", ppm%force_plsmf)
2486 : : !call ydoc%add_int1d("keep_qibz", ppm%keep_qibz)
2487 : : !call ydoc%add_int1d("has_qibz", ppm%has_qibz)
2488 : :
2489 : 3 : call ydoc%write_units_and_free(units)
2490 : :
2491 [ - + - + ]: 3 : end subroutine ppm_print
2492 : : !!***
2493 : :
2494 : : !----------------------------------------------------------------------
2495 : :
2496 [ - + - - : 414 : end module m_ppmodel
+ + - + +
- - + - +
+ - - + -
+ - + - +
- - - - -
+ - + - -
- - - + -
+ - - - -
- + ]
2497 : : !!***
|
