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