Branch data Line data Source code
1 : : !!****m* ABINIT/m_sigma
2 : : !! NAME
3 : : !! m_sigma
4 : : !!
5 : : !! FUNCTION
6 : : !! This module provides the definition of the sigma_t data type
7 : : !! used to store results of the GW calculation.
8 : : !!
9 : : !! COPYRIGHT
10 : : !! Copyright (C) 2008-2026 ABINIT group (MG, FB, GMR, VO, LR, RWG)
11 : : !! This file is distributed under the terms of the
12 : : !! GNU General Public License, see ~abinit/COPYING
13 : : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : : !!
15 : : !! SOURCE
16 : :
17 : : #if defined HAVE_CONFIG_H
18 : : #include "config.h"
19 : : #endif
20 : :
21 : : #include "abi_common.h"
22 : :
23 : : module m_sigma
24 : :
25 : : use, intrinsic :: iso_c_binding
26 : : use defs_basis
27 : : use m_gwdefs
28 : : use m_xmpi
29 : : use m_abicore
30 : : use m_errors
31 : : use m_nctk
32 : : use m_yaml
33 : : use m_melemts
34 : : use netcdf
35 : : use m_wfd
36 : :
37 : : use defs_abitypes, only : MPI_type
38 : : !use m_gwdefs, only : unt_gw, unt_sig, unt_sgr, unt_sgm, unt_gwdiag, sigparams_t, unt_sigc
39 : : use m_fstrings, only : itoa, sjoin
40 : : use m_numeric_tools, only : c2r
41 : : use m_crystal, only : crystal_t
42 : : use m_ebands, only : ebands_t
43 : : use m_bz_mesh, only : kmesh_t, littlegroup_t, findqg0
44 : : use m_screening, only : epsm1_t
45 : : use m_dtset, only : dataset_type
46 : :
47 : : implicit none
48 : :
49 : : private
50 : : !!***
51 : :
52 : : !----------------------------------------------------------------------
53 : :
54 : : !!****t* m_sigma/sigma_t
55 : : !! NAME
56 : : !! sigma_t
57 : : !!
58 : : !! FUNCTION
59 : : !! The sigma_t structured datatype gathers the results of a GW calculation.
60 : : !!
61 : : !! TODO
62 : : !! ragged arrays (nk,nsppol) --> values ?
63 : : !!
64 : : !! SOURCE
65 : :
66 : : type,public :: sigma_t
67 : :
68 : : integer :: b1gw, b2gw ! min and Max gw band indices over spin and k-points (used to dimension arrays)
69 : : integer :: gwcalctyp ! Flag defining the calculation type.
70 : : integer :: nkptgw ! No. of points calculated
71 : : integer :: nkibz ! No. of irreducible k-points.
72 : : integer :: nbnds ! Total number of bands
73 : : integer :: nomega_r ! No. of real frequencies for the spectral function.
74 : : integer :: nomega_i ! No. of frequencies along the imaginary axis.
75 : : integer :: nomega4sd ! No. of real frequencies to evaluate the derivative of $\Sigma(E)$.
76 : : integer :: nsig_ab ! 1 if nspinor=1,4 for noncollinear case.
77 : : integer :: nsppol ! No. of spin polarizations.
78 : : integer :: usepawu ! 1 if we are using DFT+U as starting point (only for PAW)
79 : :
80 : : real(dp) :: deltae ! Frequency step for the calculation of d\Sigma/dE
81 : : real(dp) :: maxomega4sd ! Max frequency around E_ks for d\Sigma/dE.
82 : : real(dp) :: maxomega_r ! Max frequency for spectral function.
83 : : real(dp) :: scissor_ene ! Scissor energy value. zero for None.
84 : :
85 : : integer,allocatable :: maxbnd(:,:)
86 : : ! (nkptgw, nsppol)
87 : : ! Max band index considered in GW for this k-point.
88 : :
89 : : integer,allocatable :: minbnd(:,:)
90 : : ! (nkptgw, nsppol)
91 : : ! Min band index considered in GW for this k-point.
92 : :
93 : : real(dp),allocatable :: degwgap(:,:)
94 : : ! (nkibz, nsppol)
95 : : ! Difference btw the QP and the KS direct gap.
96 : :
97 : : real(dp),allocatable :: egwgap(:,:)
98 : : ! (nkibz, nsppol))
99 : : ! QP direct gap at each k-point and spin.
100 : :
101 : : real(dp),allocatable :: en_qp_diago(:,:,:)
102 : : ! (nbnds, nkibz, nsppol))
103 : : ! QP energies obtained from the diagonalization of the Hermitian approximation to Sigma (QPSCGW)
104 : :
105 : : real(dp),allocatable :: e0(:,:,:)
106 : : ! (nbnds, nkibz, nsppol)
107 : : ! KS eigenvalues for each band, k-point and spin. In case of self-consistent?
108 : :
109 : : real(dp),allocatable :: e0gap(:,:)
110 : : ! (nkibz, nsppol),
111 : : ! KS gap at each k-point, for each spin.
112 : :
113 : : real(dp),allocatable :: omega_r(:)
114 : : ! (nomega_r)
115 : : ! real frequencies used for the self energy.
116 : :
117 : : real(dp),allocatable :: kptgw(:,:)
118 : : ! (3, nkptgw)
119 : : ! ! TODO there is a similar array in sigparams_t
120 : : ! List of calculated k-points.
121 : :
122 : : real(dp),allocatable :: sigxme(:,:,:)
123 : : ! (b1gw:b2gw, nkibz, nsppol*nsig_ab))
124 : : ! Diagonal matrix elements $\<nks|\Sigma_x|nks\>$
125 : :
126 : : real(dp),allocatable :: sigxcnofme(:,:,:)
127 : : ! (b1gw:b2gw, nkibz, nsppol*nsig_ab))
128 : : ! Diagonal matrix elements $\<nks|\Sigma_xc|nks\>$ taking sqrt(occs) in \Sigma_x, occs in [0,1]
129 : :
130 : : complex(dp),allocatable :: x_mat(:,:,:,:)
131 : : ! (b1gw:b2gw, b1gw:b2gw, nkibz, nsppol*nsig_ab)
132 : : ! Matrix elements of $\<nks|\Sigma_x|mks\>$
133 : :
134 : : real(dp),allocatable :: vxcme(:,:,:)
135 : : ! (b1gw:b2gw, nkibz, nsppol*nsig_ab))
136 : : ! $\<nks|v_{xc}[n_val]|nks\>$ matrix elements of vxc
137 : : ! NB: valence-only contribution i.e. computed without model core charge
138 : :
139 : : real(dp),allocatable :: vUme(:,:,:)
140 : : ! (b1gw:b2gw, nkibz, nsppol*nsig_ab))
141 : : ! $\<nks|v_{U}|nks\>$ for DFT+U.
142 : :
143 : : complex(dp),allocatable :: degw(:,:,:)
144 : : ! (b1gw:b2gw, nkibz, nsppol))
145 : : ! Difference between the QP and the KS energies.
146 : :
147 : : complex(dp),allocatable :: dsigmee0(:,:,:)
148 : : ! (b1gw:b2gw, nkibz, nsppol*nsig_ab))
149 : : ! Derivative of $\Sigma_c(E)$ calculated at the KS eigenvalue.
150 : :
151 : : complex(dp),allocatable :: egw(:,:,:)
152 : : ! (nbnds, nkibz, nsppol))
153 : : ! QP energies, $\epsilon_{nks}^{QP}$.
154 : :
155 : : logical :: needs_eigvec_qp = .True.
156 : :
157 : : ! FIXME: These arrays are huge and should be allocated only if self-consistent
158 : : complex(dp),allocatable :: eigvec_qp(:,:,:,:)
159 : : ! (nbnds, nbnds, nkibz, nsppol))
160 : : ! Expansion of the QP amplitudes in the QP basis set of the previous iteration.
161 : :
162 : : complex(dp),allocatable :: m_ks_to_qp(:,:,:,:)
163 : : ! (nbnds, nbnds, nkibz, nsppol))
164 : : ! m_ks_to_qp(ib,jb,k,s) := <\psi_{ib,k,s}^{KS}|\psi_{jb,k,s}^{QP}>
165 : :
166 : : complex(dp),allocatable :: hhartree(:,:,:,:)
167 : : ! (b1gw:b2gw, b1gw:b2gw, nkibz, nsppol*nsig_ab)
168 : : ! $\<nks|T+v_H+v_{loc}+v_{nl}|mks\>$
169 : : ! Note that v_{loc} does not include the contribution to vxc(r) given by the model core charge.
170 : :
171 : : complex(dp),allocatable :: sigcme(:,:,:,:)
172 : : ! (b1gw:b2gw, nkibz, nomega_r, nsppol*nsig_ab))
173 : : ! $\<nks|\Sigma_{c}(E)|nks\>$ at each nomega_r frequency
174 : :
175 : : complex(dp),allocatable :: sigmee(:,:,:)
176 : : ! (b1gw:b2gw, nkibz, nsppol*nsig_ab))
177 : : ! $\Sigma_{xc}E_{KS} + (E_{QP}- E_{KS})*dSigma/dE_KS
178 : :
179 : : complex(dp),allocatable :: sigcmee0(:,:,:)
180 : : ! (b1gw:b2gw, nkibz, nsppol*nsig_ab))
181 : : ! Diagonal matrix elements of $\Sigma_c(E)$ calculated at the KS energy $E_{KS}$
182 : :
183 : : complex(dp),allocatable :: sigcmesi(:,:,:,:)
184 : : ! (b1gw:b2gw, nkibz, nomega_i, nsppol*nsig_ab))
185 : : ! Matrix elements of $\Sigma_c$ along the imaginary axis.
186 : : ! Only used in case of analytical continuation.
187 : :
188 : : complex(dp),allocatable :: sigcme4sd(:,:,:,:)
189 : : ! (b1gw:b2gw, nkibz, nomega4sd, nsppol*nsig_ab))
190 : : ! Diagonal matrix elements of \Sigma_c around the zeroth order eigenvalue (usually KS).
191 : :
192 : : complex(dp),allocatable :: sigxcme(:,:,:,:)
193 : : ! (b1gw:b2gw, nkibz, nomega_r, nsppol*nsig_ab))
194 : : ! $\<nks|\Sigma_{xc}(E)|nks\>$ at each real frequency frequency.
195 : :
196 : : complex(dp),allocatable :: sigxcmesi(:,:,:,:)
197 : : ! (b1gw:b2gw, nkibz, nomega_i, nsppol*nsig_ab))
198 : : ! Matrix elements of $\Sigma_{xc}$ along the imaginary axis.
199 : : ! Only used in case of analytical continuation.
200 : :
201 : : complex(dp),allocatable :: sigxcme4sd(:,:,:,:)
202 : : ! (b1gw:b2gw, nkibz, nomega4sd, nsppol*nsig_ab))
203 : : ! Diagonal matrix elements of \Sigma_xc for frequencies around the zeroth order eigenvalues.
204 : :
205 : : complex(dp),allocatable :: ze0(:,:,:)
206 : : ! (b1gw:b2gw, nkibz, nsppol))
207 : : ! renormalization factor. $(1-\dfrac{\partial\Sigma_c} {\partial E_{KS}})^{-1}$
208 : :
209 : : complex(dp),allocatable :: omega_i(:)
210 : : ! (nomega_i)
211 : : ! Frequencies along the imaginary axis used for the analytical continuation.
212 : :
213 : : complex(dp),allocatable :: omega4sd(:,:,:,:)
214 : : ! (b1gw:b2gw, nkibz, nomega4sd, nsppol).
215 : : ! Frequencies used to evaluate the Derivative of Sigma.
216 : :
217 : : contains
218 : : procedure :: init => sigma_init
219 : : ! Initialize the object.
220 : :
221 : : procedure :: free => sigma_free
222 : : ! Deallocate memory.
223 : :
224 : : procedure :: get_exene => sigma_get_exene
225 : : ! Compute exchange energy.
226 : :
227 : : procedure :: get_excene => sigma_get_excene
228 : : ! Compute exchange-correlation MBB (Nat. Orb. Funct. Approx.) energy.
229 : :
230 : : procedure :: get_haene => sigma_get_haene
231 : : ! Compute hartree energy.
232 : :
233 : : procedure :: get_kiene => sigma_get_kiene
234 : : ! Compute kinetic energy.
235 : :
236 : : procedure :: ncwrite => sigma_ncwrite
237 : : ! Write data in netcdf format.
238 : :
239 : : procedure :: write_results => sigma_write_results
240 : : procedure :: print_perturbative => sigma_print_pertubative
241 : : procedure :: print_qpsc => sigma_print_qpsc
242 : : end type sigma_t
243 : :
244 : : public :: sigma_distribute_bks
245 : : public :: write_sigma_header
246 : : !!***
247 : :
248 : : contains !========================================================================================
249 : : !!***
250 : :
251 : : !!****f* m_sigma/write_sigma_header
252 : : !! NAME
253 : : !! write_sigma_header
254 : : !!
255 : : !! FUNCTION
256 : : !! Write basic info and dimensions used during the calculation
257 : : !! of the QP correctoions (optdriver==4).
258 : : !!
259 : : !! INPUTS
260 : : !! Sigp=sigparams_t
261 : : !! Cryst<crystal_t>= Info on the Crystal structure
262 : : !! Kmesh<kmesh_t>= Description of the BZ sampling.
263 : : !!
264 : : !! OUTPUT
265 : : !! (for writing routines, no output) otherwise, should be described
266 : : !!
267 : : !! NOTES
268 : : !!
269 : : !! SOURCE
270 : :
271 : 173 : subroutine write_sigma_header(Sigp, epsm1, Cryst, Kmesh, Qmesh)
272 : :
273 : : !Arguments ------------------------------------
274 : : !scalars
275 : : class(sigparams_t),intent(in) :: Sigp
276 : : type(kmesh_t),intent(in) :: Kmesh,Qmesh
277 : : type(crystal_t),intent(in) :: Cryst
278 : : type(epsm1_t),intent(in) :: epsm1
279 : :
280 : : !Local variables-------------------------------
281 : : !scalars
282 : : integer :: gwcalctyp, mod10
283 : : character(len=500) :: msg
284 : : integer :: units(2)
285 : : ! *************************************************************************
286 : :
287 [ + + ]: 519 : units = [std_out, ab_out]
288 : 173 : call wrtout(units, ' SIGMA fundamental parameters:')
289 : :
290 : 173 : gwcalctyp = Sigp%gwcalctyp
291 : 173 : mod10 = MOD(Sigp%gwcalctyp,10)
292 : :
293 : 87 : SELECT CASE (mod10)
294 : : CASE (SIG_GW_PPM)
295 : 87 : write(msg,'(a,i2)')' PLASMON POLE MODEL ',Sigp%ppmodel
296 : : CASE (SIG_GW_AC)
297 : 9 : write(msg,'(a)')' ANALYTIC CONTINUATION'
298 : : CASE (SIG_GW_CD)
299 : 16 : write(msg,'(a)')' CONTOUR DEFORMATION'
300 : : CASE (SIG_HF)
301 : 35 : write(msg,'(a)')' Hartree-Fock'
302 : : CASE (SIG_SEX)
303 : 1 : write(msg,'(a)')' Screened Exchange'
304 : : CASE (SIG_COHSEX)
305 : 2 : write(msg,'(a)')' COHSEX'
306 : : CASE (SIG_QPGW_PPM)
307 : 19 : write(msg,'(a,i2)')' MODEL GW with PLASMON POLE MODEL ',Sigp%ppmodel
308 : : CASE (SIG_QPGW_CD)
309 : 4 : write(msg,'(a)')' MODEL GW without PLASMON POLE MODEL'
310 : : CASE DEFAULT
311 : 0 : write(msg,'(a,i0)')' Wrong value for Sigp%gwcalctyp = ',Sigp%gwcalctyp
312 [ + + + + : 173 : ABI_BUG(msg)
+ + + +
- ]
313 : : END SELECT
314 : 173 : call wrtout(units, msg)
315 : :
316 : 173 : write(msg,'(a,i12)')' number of plane-waves for SigmaX ',Sigp%npwx
317 : 173 : call wrtout(units, msg)
318 : 173 : write(msg,'(a,i12)')' number of plane-waves for SigmaC and W ',Sigp%npwc
319 : 173 : call wrtout(units, msg)
320 : 173 : write(msg,'(a,i12)')' number of plane-waves for wavefunctions ',Sigp%npwwfn
321 : 173 : call wrtout(units, msg)
322 : 173 : write(msg,'(a,i12)')' number of bands ',Sigp%nbnds
323 : 173 : call wrtout(units, msg)
324 : 173 : write(msg,'(a,i12)')' number of independent spin polarizations ',Sigp%nsppol
325 : 173 : call wrtout(units, msg)
326 : 173 : write(msg,'(a,i12)')' number of spinorial components ',Sigp%nspinor
327 : 173 : call wrtout(units, msg)
328 : 173 : write(msg,'(a,i12)')' number of k-points in IBZ ',Kmesh%nibz
329 : 173 : call wrtout(units, msg)
330 : 173 : write(msg,'(a,i12)')' number of q-points in IBZ ',Qmesh%nibz
331 : 173 : call wrtout(units, msg)
332 : 173 : write(msg,'(a,i12)')' number of symmetry operations ',Cryst%nsym
333 : 173 : call wrtout(units, msg)
334 : 173 : write(msg,'(a,i12)')' number of k-points in BZ ',Kmesh%nbz
335 : 173 : call wrtout(units, msg)
336 : 173 : write(msg,'(a,i12)')' number of q-points in BZ ',Qmesh%nbz
337 : 173 : call wrtout(units, msg)
338 : 173 : write(msg,'(a,i12)')' number of frequencies for dSigma/dE ',Sigp%nomegasrd
339 : 173 : call wrtout(units, msg)
340 : 173 : write(msg,'(a,f12.2)')' frequency step for dSigma/dE [eV] ',Sigp%deltae*Ha_eV
341 : 173 : call wrtout(units, msg)
342 : 173 : write(msg,'(a,i12)')' number of omega for Sigma on real axis ',Sigp%nomegasr
343 : 173 : call wrtout(units, msg)
344 : 173 : write(msg,'(a,f12.2)')' max omega for Sigma on real axis [eV] ',Sigp%maxomega_r*Ha_eV
345 : 173 : call wrtout(units, msg)
346 : 173 : write(msg,'(a,f12.2)')' zcut for avoiding poles [eV] ',Sigp%zcut*Ha_eV
347 : 173 : call wrtout(units, msg)
348 : :
349 [ + + ]: 173 : if (Sigp%mbpt_sciss>0.1d-4) then
350 : 2 : write(msg,'(a,f12.2)')' scissor energy [eV] ',Sigp%mbpt_sciss*Ha_eV
351 : 2 : call wrtout(units, msg)
352 : : end if
353 : :
354 [ + + ]: 173 : if (mod10 == SIG_GW_AC) then
355 : 9 : write(msg,'(a,i12)')' number of imaginary frequencies for Sigma',Sigp%nomegasi
356 : 9 : call wrtout(units, msg)
357 : : ! MRM not needed for GW 1RDM
358 [ + + ]: 9 : if (gwcalctyp/=21) then
359 : 4 : write(msg,'(a,f12.2)')' max omega for Sigma on imag axis [eV] ',Sigp%omegasimax*Ha_eV
360 : 4 : call wrtout(units, msg)
361 : : endif
362 : : end if
363 : :
364 [ + + ]: 173 : if (Sigp%needs_w()) then
365 : 138 : write(msg,'(2a)')ch10,' EPSILON^-1 parameters (SCR file):'
366 : 138 : call wrtout(units, msg)
367 : 138 : write(msg,'(a,i12)')' dimension of the eps^-1 matrix on file ',epsm1%Hscr%npwe
368 : 138 : call wrtout(units, msg)
369 : 138 : write(msg,'(a,i12)')' dimension of the eps^-1 matrix used ',epsm1%npwe
370 : 138 : call wrtout(units, msg)
371 : 138 : write(msg,'(a,i12)')' number of plane-waves for wavefunctions ',epsm1%Hscr%npwwfn_used
372 : 138 : call wrtout(units, msg)
373 : 138 : write(msg,'(a,i12)')' number of bands ',epsm1%Hscr%nbnds_used
374 : 138 : call wrtout(units, msg)
375 : 138 : write(msg,'(a,i12)')' number of q-points in IBZ ',Qmesh%nibz
376 : 138 : call wrtout(units, msg)
377 : 138 : write(msg,'(a,i12)')' number of frequencies ',epsm1%nomega
378 : 138 : call wrtout(units, msg)
379 : 138 : write(msg,'(a,i12)')' number of real frequencies ',epsm1%nomega_r
380 : 138 : call wrtout(units, msg)
381 : 138 : write(msg,'(a,i12)')' number of imag frequencies ',epsm1%nomega_i
382 : 138 : call wrtout(units, msg)
383 : : end if
384 : :
385 : : ! MRM not needed for GW 1RDM
386 [ + + ]: 173 : if (gwcalctyp /= 21) then
387 : 168 : write(msg,'(3a)')ch10,' matrix elements of self-energy operator (all in [eV])',ch10
388 : 168 : call wrtout(units, msg)
389 : : !call wrtout(units, "(a)")" Notations:"
390 : : !call wrtout(units, "(a)")"E0: KS eigenvalue.")
391 : : !call wrtout(units, "(a)")"VxcDFT: KS exchange-correlation potential expectation value.")
392 : : !call wrtout(units, "(a)")"SigX: exchange part of the self-energy.")
393 : : !call wrtout(units, "(a)")"SigC(E0) correlation part of the self-energy, evaluated at the KS eigenenergy.")
394 : : !call wrtout(units, "(a)")"Z: renormalization factor.")
395 : : !call wrtout(units, "(a)")"dSigC/dE: energy derivative of SigC with respect to the energy.")
396 : : !call wrtout(units, "(a)")"SigC(E): correlation part of the self-energy, evaluated at the QP energy.")
397 : : !call wrtout(units, "(a)")"E-E0: difference between QP energy and KS eigenenergy.")
398 : : !call wrtout(units, "(a)")"E: quasiparticle energy.")
399 : : !if (mod10 == SIG_GW_AC) then
400 : : ! call wrtout(units, "For AC calculations, the KS Fermi level has been set to zero.")
401 : : ! call wrtout(units, "KS and QP energies are shifted accordingly.")
402 : : ! call wrtout(units, "IMPORTANT: In AC calculations, the QP energies are obtained by solving the non-linear QP equation along the real-axis")
403 : : !end if
404 : : end if
405 : :
406 [ + + ]: 168 : if (gwcalctyp < 10) then
407 : 113 : write(msg,'(a)')' Perturbative Calculation'
408 [ + + ]: 60 : else if (gwcalctyp < 20) then
409 : 7 : write(msg,'(a)')' Self-Consistent on Energies only'
410 : : else
411 : 53 : write(msg,'(a)')' Self-Consistent on Energies and Wavefunctions'
412 : : end if
413 : 173 : call wrtout(units, msg)
414 : :
415 : 173 : end subroutine write_sigma_header
416 : : !!***
417 : :
418 : : !----------------------------------------------------------------------
419 : :
420 : : !!****f* m_sigma/sigma_write_results
421 : : !! NAME
422 : : !! sigma_write_results
423 : : !!
424 : : !! FUNCTION
425 : : !! Write the final results of the GW calculation.
426 : : !!
427 : : !! INPUTS
428 : : !! Sigp=sigparams_t datatype
429 : : !! ik_ibz= index of the k-point in the array kibz, where GW corrections are calculated
430 : : !! ikcalc= index of the k-point in the array Sigp%kptgw2bz
431 : : !! ks_ebands<ebands_t>=Info on the KS band structure energies.
432 : : !!
433 : : !! SOURCE
434 : :
435 : 497 : subroutine sigma_write_results(sigma, ikcalc, ik_ibz, Sigp, ks_ebands)
436 : :
437 : : !Arguments ------------------------------------
438 : : !scalars
439 : : class(sigma_t),intent(in) :: sigma
440 : : integer,intent(in) :: ikcalc,ik_ibz
441 : : type(ebands_t),intent(in) :: ks_ebands
442 : : type(sigparams_t),intent(in) :: Sigp
443 : :
444 : : !Local variables-------------------------------
445 : : !scalars
446 : : integer :: ib,io,is,gwcalctyp,mod10
447 : : character(len=500) :: msg
448 : : type(yamldoc_t) :: ydoc
449 : : !arrays
450 : : character(len=12) :: tag_spin(2)
451 : : ! *************************************************************************
452 : :
453 : 497 : gwcalctyp = Sigp%gwcalctyp
454 : 497 : mod10 = MOD(Sigp%gwcalctyp,10)
455 : :
456 : : ! unt_gw: File with GW corrections.
457 : : ! unt_sig: Self-energy as a function of frequency.
458 : : ! unt_sgr: Derivative wrt omega of the Self-energy.
459 : : ! unt_sigc: Sigma_c(eik) MRM
460 : : ! unt_sgm: Sigma on the Matsubara axis (imag axis)
461 : :
462 [ + + + + : 1507 : tag_spin=(/' ',' '/); if (sigma%nsppol==2) tag_spin=(/', SPIN UP ',', SPIN DOWN'/)
+ + ]
463 : :
464 [ + + ]: 1002 : do is=1,sigma%nsppol
465 : 505 : write(msg,'(2a,3f8.3,a)')ch10,' k = ',Sigp%kptgw(:,ikcalc),tag_spin(is)
466 : 505 : call wrtout(std_out,msg)
467 : : !call wrtout(ab_out,msg)
468 : :
469 : 505 : msg = ' Band E0 <VxcDFT> SigX SigC(E0) Z dSigC/dE Sig(E) E-E0 E'
470 [ - + ]: 505 : if (sigma%usepawu/=0) then
471 : 0 : msg = ' Band E0 <VxcDFT> <H_U> SigX SigC(E0) Z dSigC/dE Sig(E) E-E0 E'
472 : : end if
473 : :
474 [ + + ]: 505 : if (gwcalctyp>=10) then
475 : : write(msg,'(2a)')&
476 : 297 : ' Band E_DFT <VxcDFT> E(N-1) <Hhartree> SigX SigC[E(N-1)]',&
477 : 594 : ' Z dSigC/dE Sig[E(N)] DeltaE E(N)_pert E(N)_diago'
478 : : end if
479 : 505 : call wrtout(std_out,msg)
480 : :
481 : 505 : ydoc = yamldoc_open('SelfEnergy_ee', width=11, real_fmt='(3f8.3)')
482 : 505 : call ydoc%add_real1d('kpoint', Sigp%kptgw(:,ikcalc))
483 : 505 : call ydoc%add_int('spin', is, int_fmt="(i1)")
484 : 505 : call ydoc%add_real('KS_gap', sigma%e0gap(ik_ibz,is)*Ha_eV)
485 : 505 : call ydoc%add_real('QP_gap', sigma%egwgap(ik_ibz,is)*Ha_eV)
486 : 505 : call ydoc%add_real('Delta_QP_KS', sigma%degwgap(ik_ibz,is)*Ha_eV)
487 : 505 : call ydoc%open_tabular('data', tag='SigmaeeData')
488 : 505 : call ydoc%add_tabular_line(msg)
489 : :
490 : 505 : write(unt_gw,'(3f10.6)')Sigp%kptgw(:,ikcalc)
491 : 505 : write(unt_gw,'(i4)')Sigp%maxbnd(ikcalc,is)-Sigp%minbnd(ikcalc,is)+1
492 : :
493 : 505 : write(unt_gwdiag,'(3f10.6)')Sigp%kptgw(:,ikcalc)
494 : 505 : write(unt_gwdiag,'(i4)')Sigp%maxbnd(ikcalc,is)-Sigp%minbnd(ikcalc,is)+1
495 : :
496 : 505 : write(unt_sig,'("# k = ",3f10.6)')Sigp%kptgw(:,ikcalc)
497 : 505 : write(unt_sig,'("# b = ",2i10)')Sigp%minbnd(ikcalc,is),Sigp%maxbnd(ikcalc,is)
498 : :
499 : 505 : write(unt_sgr,'("# k = ",3f10.6)')Sigp%kptgw(:,ikcalc)
500 : 505 : write(unt_sgr,'("# b = ",2i10)')Sigp%minbnd(ikcalc,is),Sigp%maxbnd(ikcalc,is)
501 : :
502 : 505 : write(unt_sigc,'("# k = ",3f10.6)')Sigp%kptgw(:,ikcalc)
503 : 505 : write(unt_sigc,'("# b = ",2i10)')Sigp%minbnd(ikcalc,is),Sigp%maxbnd(ikcalc,is)
504 : :
505 [ + + ]: 4840 : do ib=Sigp%minbnd(ikcalc,is),Sigp%maxbnd(ikcalc,is)
506 [ + + ]: 4335 : if (gwcalctyp >= 10) then
507 : 3182 : call sigma%print_QPSC(ik_ibz, ib, is, ks_ebands, units=[dev_null], ydoc=ydoc)
508 [ + + ]: 6364 : call sigma%print_QPSC(ik_ibz, ib, is, ks_ebands, units=[std_out], prtvol=1)
509 : :
510 : : write(unt_gwdiag,'(i6,3f9.4)') &
511 : 3182 : ib, &
512 : 3182 : sigma%en_qp_diago(ib,ik_ibz,is)*Ha_eV, &
513 : 3182 : (sigma%en_qp_diago(ib,ik_ibz,is) - ks_ebands%eig(ib,ik_ibz,is))*Ha_eV,&
514 : 6364 : zero
515 : :
516 : : else
517 : : ! If not ppmodel, write out also the imaginary part in ab_out
518 : 378 : select case(mod10)
519 : : case (SIG_GW_AC, SIG_GW_CD)
520 : 378 : call sigma%print_perturbative(ik_ibz, ib, is, units=[dev_null], ydoc=ydoc, prtvol=1)
521 : : case default
522 [ + + ]: 1153 : call sigma%print_perturbative(ik_ibz, ib, is, units=[dev_null], ydoc=ydoc)
523 : : end select
524 [ + + ]: 2306 : call sigma%print_perturbative(ik_ibz, ib, is, units=[std_out], prtvol=1)
525 : : end if
526 : :
527 : : write(unt_gw,'(i6,3f9.4)') &
528 : 4335 : ib, &
529 : 4335 : REAL (sigma%egw (ib,ik_ibz,is))*Ha_eV,&
530 : 4335 : REAL (sigma%degw(ib,ik_ibz,is))*Ha_eV,&
531 : 9175 : AIMAG(sigma%egw (ib,ik_ibz,is))*Ha_eV
532 : : end do !ib
533 : :
534 [ + + ]: 505 : if (sigma%e0gap(ik_ibz,is)**2+sigma%egwgap(ik_ibz,is)**2+sigma%degwgap(ik_ibz,is)**2 > tol10) then
535 : : ! Output the direct gap for each spin
536 : : ! If all the gaps are zero, this means that they could not be computed in the calling routine
537 : 485 : write(msg,'(2a,f8.3)')ch10,' E^0_gap ',sigma%e0gap(ik_ibz,is)*Ha_eV
538 : 485 : call wrtout(std_out,msg)
539 : 485 : write(msg,'(a,f8.3)') ' E^GW_gap ',sigma%egwgap(ik_ibz,is)*Ha_eV
540 : 485 : call wrtout(std_out,msg)
541 : 485 : write(msg,'(a,f8.3,a)') ' DeltaE^GW_gap ',sigma%degwgap(ik_ibz,is)*Ha_eV,ch10
542 : 485 : call wrtout(std_out,msg)
543 : : end if
544 : :
545 : 505 : call ydoc%write_and_free(ab_out)
546 : :
547 : : ! Output of the spectral function.
548 [ + + ]: 1363 : do io=1,sigma%nomega_r
549 : : write(unt_sig,'(100(e12.5,2x))')&
550 : 858 : REAL(sigma%omega_r(io))*Ha_eV,&
551 : 3740 : (REAL(sigma%sigxcme(ib,ik_ibz,io,is))*Ha_eV,&
552 : 3740 : AIMAG(sigma%sigxcme(ib,ik_ibz,io,is))*Ha_eV,&
553 : 4598 : gw_spectral_function(sigma,io,ib,ik_ibz,is),&
554 [ + - + + ]: 6819 : ib=Sigp%minbnd(ikcalc,is),Sigp%maxbnd(ikcalc,is))
555 : : end do
556 : :
557 [ + + ]: 4840 : do ib=Sigp%minbnd(ikcalc,is),Sigp%maxbnd(ikcalc,is)
558 : 4335 : write(unt_sgr,'("# ik, ib",2i5)')ik_ibz,ib
559 [ + + ]: 19279 : do io=1,sigma%nomega4sd
560 : : write(unt_sgr,'(100(e12.5,2x))') &
561 : 14439 : REAL (sigma%omega4sd (ib,ik_ibz,io,is)) * Ha_eV,&
562 : 14439 : REAL (sigma%sigxcme4sd(ib,ik_ibz,io,is)) * Ha_eV,&
563 : 33213 : AIMAG(sigma%sigxcme4sd(ib,ik_ibz,io,is)) * Ha_eV
564 : : end do
565 : : end do
566 : :
567 : : !MRM
568 [ + + ]: 4840 : do ib=Sigp%minbnd(ikcalc,is),Sigp%maxbnd(ikcalc,is)
569 : 4335 : write(unt_sigc,'("# ik, ib",2i5)')ik_ibz,ib
570 [ + + ]: 19279 : do io=1,sigma%nomega4sd
571 : : write(unt_sigc,'(100(e12.5,2x))') &
572 : 14439 : REAL (sigma%omega4sd (ib,ik_ibz,io,is)) * Ha_eV,&
573 : 14439 : REAL (sigma%sigcme4sd(ib,ik_ibz,io,is)) * Ha_eV,&
574 : 33213 : AIMAG(sigma%sigcme4sd(ib,ik_ibz,io,is)) * Ha_eV
575 : : end do
576 : : end do
577 : :
578 [ + + ]: 1002 : if (mod10 == SIG_GW_AC) then
579 : : ! For AC, write sigma matrix elements along the imaginary axis
580 [ + + ]: 60 : do ib=Sigp%minbnd(ikcalc,is),Sigp%maxbnd(ikcalc,is)
581 : 53 : write(unt_sgm,'("# ik, ib",2i5)')ik_ibz,ib
582 [ + + ]: 670 : do io=1,sigma%nomega_i
583 : : write(unt_sgm,'(3(e12.5,2x))') &
584 : 610 : AIMAG(sigma%omega_i(io)) * Ha_eV,&
585 : 610 : REAL (sigma%sigxcmesi(ib,ik_ibz,io,is)) * Ha_eV,&
586 : 1273 : AIMAG(sigma%sigxcmesi(ib,ik_ibz,io,is)) * Ha_eV
587 : : end do
588 : : end do
589 : : end if
590 : :
591 : : end do !is
592 : :
593 : 497 : end subroutine sigma_write_results
594 : : !!***
595 : :
596 : : !----------------------------------------------------------------------
597 : :
598 : : !!****f* m_sigma/gw_spectral_function
599 : : !! NAME
600 : : !! gw_spectral_function
601 : : !!
602 : : !! FUNCTION
603 : : !! Compute the spectral function
604 : : !!
605 : : !! INPUTS
606 : : !! io,ib,ik_ibz,spin=Frequency, band, k-point, spin index
607 : : !!
608 : : !! SOURCE
609 : :
610 : 3740 : real(dp) pure function gw_spectral_function(sigma, io, ib, ik_ibz, spin) result(aw)
611 : :
612 : : !Arguments ------------------------------------
613 : : class(sigma_t),intent(in) :: sigma
614 : : integer,intent(in) :: io,ib,ik_ibz,spin
615 : : ! *********************************************************************
616 : :
617 : : aw = one / pi * abs(aimag(sigma%sigcme(ib,ik_ibz,io,spin))) &
618 : : /( (real(sigma%omega_r(io) - sigma%hhartree(ib,ib,ik_ibz,spin) - sigma%sigxcme(ib,ik_ibz,io,spin)))**2 &
619 : 3740 : +(aimag(sigma%sigcme(ib,ik_ibz,io,spin))) ** 2) / Ha_eV
620 : :
621 : 3740 : end function gw_spectral_function
622 : : !!***
623 : :
624 : : !----------------------------------------------------------------------
625 : :
626 : : !!****f* m_sigma/sigma_print_pertubative
627 : : !! NAME
628 : : !! sigma_print_pertubative
629 : : !!
630 : : !! FUNCTION
631 : : !! Write the results of the GW calculation done with the perturbative approach
632 : : !!
633 : : !! INPUTS
634 : : !!
635 : : !! OUTPUT
636 : : !!
637 : : !! SOURCE
638 : :
639 [ + - + - : 2306 : subroutine sigma_print_pertubative(sigma, ik_ibz, band, spin, units, &
- + ]
640 : : prtvol, with_header, ydoc) ! Optional
641 : :
642 : : !Arguments ------------------------------------
643 : : !scalars
644 : : class(sigma_t),intent(in) :: sigma
645 : : integer,intent(in) :: band,ik_ibz,spin
646 : : integer,optional,intent(in) :: prtvol,units(:)
647 : : logical,optional,intent(in) :: with_header
648 : : type(yamldoc_t),intent(inout),optional :: ydoc
649 : :
650 : : !Local variables-------------------------------
651 : : !scalars
652 : : integer :: verbose
653 : : character(len=500) :: msg
654 : : ! *********************************************************************
655 : :
656 [ + + ]: 2306 : verbose=0 ; if (PRESENT(prtvol)) verbose=prtvol
657 : :
658 [ - + ]: 2306 : if (present(with_header)) then
659 [ # # ]: 0 : if (with_header) then
660 : 0 : call wrtout(units,' Band E0 <VxcDFT> SigX SigC(E0) Z dSigC/dE Sig(E) E-E0 E ')
661 : : end if
662 : : end if
663 : :
664 [ + - ]: 2306 : if (sigma%usepawu == 0) then
665 : :
666 [ + + ]: 2306 : if (sigma%nsig_ab /= 1) then
667 : : write(msg,'(i5,9f8.3)') &
668 : 288 : band, &
669 : 288 : sigma%e0 (band,ik_ibz,1)*Ha_eV, &
670 [ + + ]: 1440 : SUM(sigma%vxcme (band,ik_ibz,:))*Ha_eV, &
671 [ + + ]: 1440 : SUM(sigma%sigxme (band,ik_ibz,:))*Ha_eV, &
672 [ + + ]: 1440 : REAL(SUM(sigma%sigcmee0(band,ik_ibz,:)))*Ha_eV,&
673 : 288 : REAL(sigma%ze0 (band,ik_ibz,1)), &
674 [ + + ]: 1440 : REAL(SUM(sigma%dsigmee0(band,ik_ibz,:))), &
675 [ + + ]: 1440 : REAL(SUM(sigma%sigmee (band,ik_ibz,:)))*Ha_eV,&
676 : 288 : REAL(sigma%degw (band,ik_ibz,1))*Ha_eV, &
677 : 576 : REAL(sigma%egw (band,ik_ibz,1))*Ha_eV
678 : 288 : call wrtout(units, msg)
679 [ + + ]: 288 : if (present(ydoc)) call ydoc%add_tabular_line(msg)
680 [ + + ]: 288 : if (verbose /= 0) then
681 : : write(msg,'(i5,9f8.3)') &
682 : 198 : band, &
683 : 198 : zero, &
684 : 198 : zero, &
685 : 198 : zero, &
686 [ + + ]: 990 : AIMAG(SUM(sigma%sigcmee0(band,ik_ibz,:)))*Ha_eV,&
687 : 198 : AIMAG(sigma%ze0 (band,ik_ibz,1)), &
688 [ + + ]: 990 : AIMAG(SUM(sigma%dsigmee0(band,ik_ibz,:))), &
689 [ + + ]: 990 : AIMAG(SUM(sigma%sigmee (band,ik_ibz,:)))*Ha_eV,&
690 : 198 : AIMAG(sigma%degw (band,ik_ibz,1))*Ha_eV, &
691 : 396 : AIMAG(sigma%egw (band,ik_ibz,1))*Ha_eV
692 : 198 : call wrtout(units, msg)
693 [ + + ]: 198 : if(present(ydoc)) call ydoc%add_tabular_line(msg)
694 : : end if
695 : : else
696 : : write(msg,'(i5,9f8.3)') &
697 : 2018 : band, &
698 : 2018 : sigma%e0 (band,ik_ibz,spin)*Ha_eV, &
699 : 2018 : sigma%vxcme (band,ik_ibz,spin)*Ha_eV, &
700 : 2018 : sigma%sigxme (band,ik_ibz,spin)*Ha_eV, &
701 : 2018 : REAL(sigma%sigcmee0(band,ik_ibz,spin))*Ha_eV,&
702 : 2018 : REAL(sigma%ze0 (band,ik_ibz,spin)), &
703 : 2018 : REAL(sigma%dsigmee0(band,ik_ibz,spin)), &
704 : 2018 : REAL(sigma%sigmee (band,ik_ibz,spin))*Ha_eV,&
705 : 2018 : REAL(sigma%degw (band,ik_ibz,spin))*Ha_eV,&
706 : 4036 : REAL(sigma%egw (band,ik_ibz,spin))*Ha_eV
707 : 2018 : call wrtout(units, msg)
708 [ + + ]: 2018 : if (present(ydoc)) call ydoc%add_tabular_line(msg)
709 : :
710 [ + + ]: 2018 : if (verbose /= 0) then
711 : : write(msg,'(i5,9f8.3)') &
712 : 1333 : band, &
713 : 1333 : zero, &
714 : 1333 : zero, &
715 : 1333 : zero, &
716 : 1333 : AIMAG(sigma%sigcmee0(band,ik_ibz,spin))*Ha_eV,&
717 : 1333 : AIMAG(sigma%ze0 (band,ik_ibz,spin)), &
718 : 1333 : AIMAG(sigma%dsigmee0(band,ik_ibz,spin)), &
719 : 1333 : AIMAG(sigma%sigmee (band,ik_ibz,spin))*Ha_eV,&
720 : 1333 : AIMAG(sigma%degw (band,ik_ibz,spin))*Ha_eV,&
721 : 2666 : AIMAG(sigma%egw (band,ik_ibz,spin))*Ha_eV
722 : 1333 : call wrtout(units,msg)
723 [ + + ]: 1333 : if (present(ydoc)) call ydoc%add_tabular_line(msg)
724 : : end if
725 : : end if
726 : :
727 : : else
728 : : ! PAW+U+GW calculation.
729 [ # # ]: 0 : ABI_CHECK(sigma%nsig_ab==1, 'DFT+U with spinor not implemented')
730 : : write(msg,'(i5,10f8.3)') &
731 : 0 : band, &
732 : 0 : sigma%e0 (band,ik_ibz,spin)*Ha_eV, &
733 : 0 : sigma%vxcme (band,ik_ibz,spin)*Ha_eV, &
734 : 0 : sigma%vUme (band,ik_ibz,spin)*Ha_eV, &
735 : 0 : sigma%sigxme (band,ik_ibz,spin)*Ha_eV, &
736 : 0 : REAL(sigma%sigcmee0(band,ik_ibz,spin))*Ha_eV,&
737 : 0 : REAL(sigma%ze0 (band,ik_ibz,spin)), &
738 : 0 : REAL(sigma%dsigmee0(band,ik_ibz,spin)), &
739 : 0 : REAL(sigma%sigmee (band,ik_ibz,spin))*Ha_eV,&
740 : 0 : REAL(sigma%degw (band,ik_ibz,spin))*Ha_eV,&
741 : 0 : REAL(sigma%egw (band,ik_ibz,spin))*Ha_eV
742 : 0 : call wrtout(units,msg)
743 [ # # ]: 0 : if(present(ydoc)) call ydoc%add_tabular_line(msg)
744 : :
745 [ # # ]: 0 : if (verbose/=0) then
746 : : write(msg,'(i5,10f8.3)') &
747 : 0 : band, &
748 : 0 : zero, &
749 : 0 : zero, &
750 : 0 : zero, &
751 : 0 : zero, &
752 : 0 : AIMAG(sigma%sigcmee0(band,ik_ibz,spin))*Ha_eV,&
753 : 0 : AIMAG(sigma%ze0 (band,ik_ibz,spin)), &
754 : 0 : AIMAG(sigma%dsigmee0(band,ik_ibz,spin)), &
755 : 0 : AIMAG(sigma%sigmee (band,ik_ibz,spin))*Ha_eV,&
756 : 0 : AIMAG(sigma%degw (band,ik_ibz,spin))*Ha_eV,&
757 : 0 : AIMAG(sigma%egw (band,ik_ibz,spin))*Ha_eV
758 : 0 : call wrtout(units, msg)
759 [ # # ]: 0 : if(present(ydoc)) call ydoc%add_tabular_line(msg)
760 : : end if
761 : : end if
762 : :
763 : 2306 : end subroutine sigma_print_pertubative
764 : : !!***
765 : :
766 : : !----------------------------------------------------------------------
767 : :
768 : : !!****f* m_sigma/sigma_print_qpsc
769 : : !! NAME
770 : : !! sigma_print_qpsc
771 : : !!
772 : : !! FUNCTION
773 : : !! Write the results of the GW calculation in case of self-consistency
774 : : !!
775 : : !! SOURCE
776 : :
777 [ + - + - : 6364 : subroutine sigma_print_qpsc(sigma, ik_ibz, band, spin, ks_ebands, units, &
- + ]
778 : : prtvol, ydoc) ! Optional
779 : :
780 : : !Arguments ------------------------------------
781 : : !scalars
782 : : class(sigma_t),intent(in) :: sigma
783 : : integer,intent(in) :: band,ik_ibz,spin
784 : : integer,intent(in),optional :: prtvol, units(:)
785 : : type(ebands_t),intent(in) :: ks_ebands
786 : : type(yamldoc_t),intent(inout),optional :: ydoc
787 : :
788 : : !Local variables-------------------------------
789 : : !scalars
790 : : integer :: verbose
791 : : character(len=500) :: msg
792 : : ! *********************************************************************
793 : :
794 [ + + ]: 6364 : verbose=0 ; if (PRESENT(prtvol )) verbose=prtvol
795 : :
796 : : ! write(msg,'(a)')&
797 : : !& ' Band E_DFT <VxcDFT> E(N-1) <Hhartree> SigX SigC[E(N-1)]',&
798 : : !& ' Z dSigC/dE Sig[E(N)] DeltaE E(N)_pert E(N)_diago'
799 : :
800 : : if (sigma%usepawu==0 .or. .TRUE.) then
801 [ - + ]: 6364 : if (sigma%nsig_ab/=1) then
802 : : write(msg,'(i5,12(2x,f8.3))') &
803 : 0 : band, &
804 : 0 : ks_ebands%eig (band,ik_ibz,1)*Ha_eV, &
805 [ # # ]: 0 : SUM(sigma%vxcme (band,ik_ibz,:))*Ha_eV, &
806 : 0 : sigma%e0 (band,ik_ibz,1)*Ha_eV, &
807 [ # # ]: 0 : REAL(SUM(sigma%hhartree(band,band,ik_ibz,:)))*Ha_eV,&
808 [ # # ]: 0 : SUM(sigma%sigxme (band,ik_ibz,:))*Ha_eV, &
809 [ # # ]: 0 : REAL(SUM(sigma%sigcmee0(band,ik_ibz,:)))*Ha_eV, &
810 : 0 : REAL(sigma%ze0 (band,ik_ibz,1)), &
811 [ # # ]: 0 : REAL(SUM(sigma%dsigmee0(band,ik_ibz,:))), &
812 [ # # ]: 0 : REAL(SUM(sigma%sigmee (band,ik_ibz,:)))*Ha_eV, &
813 : 0 : REAL(sigma%degw (band,ik_ibz,1))*Ha_eV, &
814 : 0 : REAL(sigma%egw (band,ik_ibz,1))*Ha_eV, &
815 : 0 : sigma%en_qp_diago (band,ik_ibz,1)*Ha_eV
816 : 0 : call wrtout(units, msg)
817 [ # # ]: 0 : if (present(ydoc)) call ydoc%add_tabular_line(msg)
818 : :
819 : : write(msg,'(i5,12(2x,f8.3))') &
820 : 0 : band, &
821 : 0 : zero, &
822 : 0 : zero, &
823 : 0 : zero, &
824 [ # # ]: 0 : AIMAG(SUM(sigma%hhartree(band,band,ik_ibz,:)))*Ha_eV,&
825 : 0 : zero, &
826 [ # # ]: 0 : AIMAG(SUM(sigma%sigcmee0(band,ik_ibz,:)))*Ha_eV, &
827 : 0 : AIMAG(sigma%ze0 (band,ik_ibz,1)), &
828 [ # # ]: 0 : AIMAG(SUM(sigma%dsigmee0(band,ik_ibz,:))), &
829 [ # # ]: 0 : AIMAG(SUM(sigma%sigmee (band,ik_ibz,:)))*Ha_eV, &
830 : 0 : AIMAG(sigma%degw (band,ik_ibz,1))*Ha_eV, &
831 : 0 : AIMAG(sigma%egw (band,ik_ibz,1))*Ha_eV, &
832 : 0 : zero
833 [ # # ]: 0 : if (verbose/=0) then
834 : 0 : call wrtout(units, msg)
835 [ # # ]: 0 : if (present(ydoc)) call ydoc%add_tabular_line(msg)
836 : : end if
837 : : else
838 : : write(msg,'(i5,12(2x,f8.3))') &
839 : 6364 : band, &
840 : 6364 : ks_ebands%eig (band,ik_ibz,spin)*Ha_eV, &
841 : 6364 : sigma%vxcme (band,ik_ibz,spin)*Ha_eV, &
842 : 6364 : sigma%e0 (band,ik_ibz,spin)*Ha_eV, &
843 : 6364 : REAL(sigma%hhartree (band,band,ik_ibz,spin))*Ha_eV, &
844 : 6364 : sigma%sigxme (band,ik_ibz,spin)*Ha_eV, &
845 : 6364 : REAL(sigma%sigcmee0 (band,ik_ibz,spin))*Ha_eV, &
846 : 6364 : REAL(sigma%ze0 (band,ik_ibz,spin)), &
847 : 6364 : REAL(sigma%dsigmee0 (band,ik_ibz,spin)), &
848 : 6364 : REAL(sigma%sigmee (band,ik_ibz,spin))*Ha_eV, &
849 : 6364 : REAL(sigma%degw (band,ik_ibz,spin))*Ha_eV, &
850 : 6364 : REAL(sigma%egw (band,ik_ibz,spin))*Ha_eV, &
851 : 12728 : sigma%en_qp_diago(band,ik_ibz,spin)*Ha_eV
852 : 6364 : call wrtout(units, msg)
853 [ + + ]: 6364 : if (present(ydoc)) call ydoc%add_tabular_line(msg)
854 : :
855 : : write(msg,'(i5,12(2x,f8.3))') &
856 : 6364 : band, &
857 : 6364 : zero, &
858 : 6364 : zero, &
859 : 6364 : zero, &
860 : 6364 : AIMAG(sigma%hhartree (band,band,ik_ibz,spin))*Ha_eV,&
861 : 6364 : zero, &
862 : 6364 : AIMAG(sigma%sigcmee0 (band,ik_ibz,spin))*Ha_eV, &
863 : 6364 : AIMAG(sigma%ze0 (band,ik_ibz,spin)), &
864 : 6364 : AIMAG(sigma%dsigmee0 (band,ik_ibz,spin)), &
865 : 6364 : AIMAG(sigma%sigmee (band,ik_ibz,spin))*Ha_eV, &
866 : 6364 : AIMAG(sigma%degw (band,ik_ibz,spin))*Ha_eV, &
867 : 6364 : AIMAG(sigma%egw (band,ik_ibz,spin))*Ha_eV, &
868 : 12728 : zero
869 [ + + ]: 6364 : if (verbose/=0) then
870 : 3182 : call wrtout(units, msg)
871 [ - + ]: 3182 : if (present(ydoc)) call ydoc%add_tabular_line(msg)
872 : : end if
873 : : end if
874 : :
875 : : else
876 : : ! PAW+U+GW calculation.
877 : : ABI_ERROR("PAW+U+GW not yet implemented")
878 : : end if
879 : :
880 : 6364 : end subroutine sigma_print_qpsc
881 : : !!***
882 : :
883 : : !----------------------------------------------------------------------
884 : :
885 : : !!****f* m_sigma/sigma_init
886 : : !! NAME
887 : : !! sigma_init
888 : : !!
889 : : !! FUNCTION
890 : : !! Main creation method for the sigma_t data type.
891 : : !!
892 : : !! INPUTS
893 : : !! usepawu= /=0 if we used DFT+U as starting point (only for PAW)
894 : : !!
895 : : !! SOURCE
896 : :
897 : 201 : subroutine sigma_init(sigma, Sigp, nkibz, usepawu)
898 : :
899 : : !Arguments ------------------------------------
900 : : class(sigma_t),intent(inout) :: sigma
901 : : type(sigparams_t),intent(in) :: Sigp
902 : : integer,intent(in) :: nkibz, usepawu
903 : :
904 : : !Local variables-------------------------------
905 : : integer :: b1gw,b2gw,mod10
906 : : ! *************************************************************************
907 : :
908 : 201 : mod10 = MOD(Sigp%gwcalctyp,10)
909 : :
910 : : ! Copy important dimensions
911 : 201 : sigma%nkptgw =Sigp%nkptgw
912 : 201 : sigma%gwcalctyp =Sigp%gwcalctyp
913 : 201 : sigma%deltae =Sigp%deltae
914 : 201 : sigma%maxomega4sd=Sigp%maxomega4sd
915 : 201 : sigma%maxomega_r =Sigp%maxomega_r
916 : 201 : sigma%scissor_ene=Sigp%mbpt_sciss
917 : :
918 : : !FIXME this should be done in sigma_allocate
919 [ - + + - : 804 : ABI_MALLOC(sigma%minbnd, (sigma%nkptgw,Sigp%nsppol))
- + + - +
- + - - +
- + - + ]
920 [ - + + - : 603 : ABI_MALLOC(sigma%maxbnd, (sigma%nkptgw,Sigp%nsppol))
- + + - -
+ - + -
+ ]
921 [ + - + - : 2307 : sigma%minbnd=Sigp%minbnd; sigma%maxbnd=Sigp%maxbnd
- + - - -
- - - + +
+ + + - +
- - + - -
- - - - +
+ + + ]
922 [ - + + - : 603 : ABI_MALLOC(sigma%kptgw, (3,sigma%nkptgw))
+ - + - -
+ - + -
+ ]
923 [ + - + - : 2958 : sigma%kptgw=Sigp%kptgw
- + - - -
- - - + +
+ + ]
924 : :
925 : 201 : sigma%b1gw =Sigp%minbdgw ! min and Max GW band index over k and spin.
926 : 201 : sigma%b2gw =Sigp%maxbdgw ! Used to dimension arrays.
927 : 201 : sigma%nbnds =Sigp%nbnds
928 : 201 : sigma%nkibz =nkibz
929 : 201 : sigma%nsppol =Sigp%nsppol
930 : 201 : sigma%nsig_ab =Sigp%nsig_ab
931 : 201 : sigma%nomega_r =Sigp%nomegasr !FIXME change name
932 : 201 : sigma%nomega_i =Sigp%nomegasi
933 : 201 : sigma%nomega4sd=Sigp%nomegasrd
934 : 201 : sigma%usepawu =usepawu
935 : :
936 : : !================================================
937 : : ! === Allocate arrays in the sigma_t datatype ===
938 : : !================================================
939 : 201 : b1gw=sigma%b1gw
940 : 201 : b2gw=sigma%b2gw
941 : :
942 : : ! hhartree(b1,b2,k,s)= <b1,k,s|T+v_{loc}+v_{nl}+v_{H}|b2,k,s>
943 [ - + + - : 106123 : ABI_CALLOC(sigma%hhartree, (b1gw:b2gw,b1gw:b2gw,sigma%nkibz,sigma%nsppol*sigma%nsig_ab))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + + + +
+ + + +
+ ]
944 : :
945 : : ! QP amplitudes and energies.
946 [ - + + - : 31612 : ABI_CALLOC(sigma%en_qp_diago, (sigma%nbnds,sigma%nkibz,sigma%nsppol))
- + + - -
+ + - + -
+ - - + -
+ - + + +
+ + + + ]
947 : :
948 : 201 : sigma%needs_eigvec_qp = sigp%gwcalctyp >= 10
949 [ + + ]: 201 : if (sigma%needs_eigvec_qp) then
950 [ - + + - : 73727 : ABI_CALLOC(sigma%eigvec_qp, (sigma%nbnds,sigma%nbnds,sigma%nkibz,sigma%nsppol))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + + + +
+ + + +
+ ]
951 : : end if
952 : :
953 [ - + + - : 12546 : ABI_CALLOC(sigma%vxcme, (b1gw:b2gw, sigma%nkibz, sigma%nsppol*sigma%nsig_ab))
- + + - -
+ + - + -
+ - - + -
+ - + + +
+ + + + ]
954 [ - + + - : 12345 : ABI_CALLOC(sigma%vUme, (b1gw:b2gw, sigma%nkibz, sigma%nsppol*sigma%nsig_ab))
- + + - -
+ + - - +
- + - + +
+ + + +
+ ]
955 [ - + + - : 12345 : ABI_CALLOC(sigma%sigxme, (b1gw:b2gw, sigma%nkibz, sigma%nsppol*sigma%nsig_ab))
- + + - -
+ + - - +
- + - + +
+ + + +
+ ]
956 [ - + + - : 12345 : ABI_CALLOC(sigma%sigxcnofme, (b1gw:b2gw, sigma%nkibz, sigma%nsppol*sigma%nsig_ab))
- + + - -
+ + - - +
- + - + +
+ + + +
+ ]
957 [ - + + - : 105922 : ABI_CALLOC(sigma%x_mat, (b1gw:b2gw, b1gw:b2gw, sigma%nkibz, sigma%nsppol*sigma%nsig_ab))
- + + - -
+ + - - +
+ - - + -
+ - + + +
+ + + + +
+ ]
958 [ - + + - : 61708 : ABI_CALLOC(sigma%sigcme, (b1gw:b2gw, sigma%nkibz, sigma%nomega_r, sigma%nsppol*sigma%nsig_ab))
- + + - +
+ + - - +
+ - + - +
+ - + - +
- + + + +
+ + + +
+ ]
959 [ - + + - : 61507 : ABI_CALLOC(sigma%sigxcme, (b1gw:b2gw, sigma%nkibz, sigma%nomega_r, sigma%nsppol*sigma%nsig_ab))
- + + - +
+ + - - +
+ - - + -
+ - + + +
+ + + + +
+ ]
960 [ - + + - : 12546 : ABI_CALLOC(sigma%sigcmee0, (b1gw:b2gw, sigma%nkibz, sigma%nsppol*sigma%nsig_ab))
- + + - -
+ + - + -
+ - - + -
+ - + + +
+ + + + ]
961 [ - + + - : 11226 : ABI_CALLOC(sigma%ze0, (b1gw:b2gw, sigma%nkibz, sigma%nsppol))
- + + - -
+ + - + -
+ - - + -
+ - + + +
+ + + + ]
962 [ - + + - : 12345 : ABI_CALLOC(sigma%dsigmee0, (b1gw:b2gw, sigma%nkibz, sigma%nsppol*sigma%nsig_ab))
- + + - -
+ + - - +
- + - + +
+ + + +
+ ]
963 [ - + + - : 12345 : ABI_CALLOC(sigma%sigmee, (b1gw:b2gw, sigma%nkibz, sigma%nsppol*sigma%nsig_ab))
- + + - -
+ + - - +
- + - + +
+ + + +
+ ]
964 [ - + + - : 11025 : ABI_CALLOC(sigma%degw, (b1gw:b2gw, sigma%nkibz, sigma%nsppol))
- + + - -
+ + - - +
- + - + +
+ + + +
+ ]
965 [ - + + - : 31411 : ABI_CALLOC(sigma%e0, (sigma%nbnds, sigma%nkibz, sigma%nsppol))
- + + - -
+ + - - +
- + - + +
+ + + +
+ ]
966 [ - + + - : 31612 : ABI_CALLOC(sigma%egw, (sigma%nbnds, sigma%nkibz, sigma%nsppol))
- + + - -
+ + - + -
+ - - + -
+ - + + +
+ + + + ]
967 [ - + + - : 2241 : ABI_CALLOC(sigma%e0gap, (sigma%nkibz, sigma%nsppol))
- + + - +
- + - - +
- + - + +
+ + + ]
968 [ - + + - : 2040 : ABI_CALLOC(sigma%degwgap, (sigma%nkibz, sigma%nsppol))
- + + - -
+ - + - +
+ + + + ]
969 [ - + + - : 2040 : ABI_CALLOC(sigma%egwgap, (sigma%nkibz, sigma%nsppol))
- + + - -
+ - + - +
+ + + + ]
970 : :
971 : : ! These quantities are used to evaluate $\Sigma(E)$ around the KS\QP eigenvalue
972 [ - + + - : 54312 : ABI_CALLOC(sigma%omega4sd, (b1gw:b2gw, sigma%nkibz, sigma%nomega4sd, sigma%nsppol))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + + + +
+ + + +
+ ]
973 [ - + + - : 60069 : ABI_CALLOC(sigma%sigcme4sd, (b1gw:b2gw, sigma%nkibz, sigma%nomega4sd, sigma%nsppol*sigma%nsig_ab))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + + + +
+ + + +
+ ]
974 [ - + + - : 59868 : ABI_CALLOC(sigma%sigxcme4sd, (b1gw:b2gw, sigma%nkibz, sigma%nomega4sd, sigma%nsppol*sigma%nsig_ab))
- + + - -
+ + - - +
+ - - + -
+ - + + +
+ + + + +
+ ]
975 : :
976 : : ! Mesh along the real axis.
977 : : !TODO Find better treatment
978 [ + + ]: 201 : if (sigma%nomega_r > 0) then
979 [ - + + - : 15 : ABI_MALLOC(sigma%omega_r, (sigma%nomega_r))
+ - - + -
+ - + ]
980 [ + + ]: 813 : sigma%omega_r(:)=Sigp%omega_r(:)
981 : : end if
982 : :
983 : : ! Analytic Continuation
984 : : ! FIXME omegasi should not be in Sigp% here we should construct the mesh
985 [ + + ]: 201 : if (mod10 == SIG_GW_AC) then
986 [ - + + - : 27 : ABI_MALLOC(sigma%omega_i, (sigma%nomega_i))
+ - + - -
+ - + -
+ ]
987 [ + - - + : 228 : sigma%omega_i = Sigp%omegasi
- - - - +
+ ]
988 [ - + + - : 11493 : ABI_CALLOC(sigma%sigcmesi, (b1gw:b2gw, sigma%nkibz, sigma%nomega_i, sigma%nsppol*sigma%nsig_ab))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + + + +
+ + + +
+ ]
989 [ - + + - : 11484 : ABI_CALLOC(sigma%sigxcmesi, (b1gw:b2gw, sigma%nkibz, sigma%nomega_i, sigma%nsppol*sigma%nsig_ab))
- + + - -
+ + - - +
+ - - + -
+ - + + +
+ + + + +
+ ]
990 : : end if
991 : :
992 : 201 : end subroutine sigma_init
993 : : !!***
994 : :
995 : : !----------------------------------------------------------------------
996 : :
997 : : !!****f* m_sigma/sigma_free
998 : : !! NAME
999 : : !! sigma_free
1000 : : !!
1001 : : !! FUNCTION
1002 : : !! Deallocate all associated pointers defined in the sigma_t data type.
1003 : : !!
1004 : : !! SOURCE
1005 : :
1006 : 201 : subroutine sigma_free(sigma)
1007 : :
1008 : : !Arguments ------------------------------------
1009 : : class(sigma_t),intent(inout) :: sigma
1010 : : ! *************************************************************************
1011 : :
1012 : : ! integer
1013 [ + - ]: 201 : ABI_SFREE(sigma%maxbnd)
1014 [ + - ]: 201 : ABI_SFREE(sigma%minbnd)
1015 : :
1016 : : ! real
1017 [ + - ]: 201 : ABI_SFREE(sigma%degwgap)
1018 [ + - ]: 201 : ABI_SFREE(sigma%egwgap)
1019 [ + - ]: 201 : ABI_SFREE(sigma%en_qp_diago)
1020 [ + - ]: 201 : ABI_SFREE(sigma%e0)
1021 [ + - ]: 201 : ABI_SFREE(sigma%e0gap)
1022 [ + + ]: 201 : ABI_SFREE(sigma%omega_r)
1023 [ + - ]: 201 : ABI_SFREE(sigma%kptgw)
1024 [ + - ]: 201 : ABI_SFREE(sigma%sigxme)
1025 [ + - ]: 201 : ABI_SFREE(sigma%sigxcnofme)
1026 [ + - ]: 201 : ABI_SFREE(sigma%x_mat)
1027 [ + - ]: 201 : ABI_SFREE(sigma%vxcme)
1028 [ + - ]: 201 : ABI_SFREE(sigma%vUme)
1029 : :
1030 : : ! complex
1031 [ + - ]: 201 : ABI_SFREE(sigma%degw)
1032 [ + - ]: 201 : ABI_SFREE(sigma%dsigmee0)
1033 [ + - ]: 201 : ABI_SFREE(sigma%egw)
1034 [ + + ]: 201 : ABI_SFREE(sigma%eigvec_qp)
1035 [ + + ]: 201 : ABI_SFREE(sigma%m_ks_to_qp)
1036 [ + - ]: 201 : ABI_SFREE(sigma%hhartree)
1037 [ + - ]: 201 : ABI_SFREE(sigma%sigcme)
1038 [ + - ]: 201 : ABI_SFREE(sigma%sigmee)
1039 [ + - ]: 201 : ABI_SFREE(sigma%sigcmee0)
1040 [ + + ]: 201 : ABI_SFREE(sigma%sigcmesi)
1041 [ + - ]: 201 : ABI_SFREE(sigma%sigcme4sd)
1042 [ + - ]: 201 : ABI_SFREE(sigma%sigxcme)
1043 [ + + ]: 201 : ABI_SFREE(sigma%sigxcmesi)
1044 [ + - ]: 201 : ABI_SFREE(sigma%sigxcme4sd)
1045 [ + - ]: 201 : ABI_SFREE(sigma%ze0)
1046 [ + + ]: 201 : ABI_SFREE(sigma%omega_i)
1047 [ + - ]: 201 : ABI_SFREE(sigma%omega4sd)
1048 : :
1049 : 201 : end subroutine sigma_free
1050 : : !!***
1051 : :
1052 : : !----------------------------------------------------------------------
1053 : :
1054 : : !!****f* m_sigma/sigma_get_exene
1055 : : !! NAME
1056 : : !! sigma_get_exene
1057 : : !!
1058 : : !! FUNCTION
1059 : : !! Compute exchange energy.
1060 : : !!
1061 : : !! INPUTS
1062 : : !! sigma<sigma_t>=Sigma results
1063 : : !! kmesh<kmesh_t>=BZ sampling.
1064 : : !! bands<band_t>=Bands with occupation factors
1065 : : !!
1066 : : !! SOURCE
1067 : :
1068 : 84 : real(dp) pure function sigma_get_exene(sigma, kmesh, bands) result(ex_energy)
1069 : :
1070 : : !Arguments ------------------------------------
1071 : : class(sigma_t),intent(in) :: sigma
1072 : : type(kmesh_t),intent(in) :: kmesh
1073 : : type(ebands_t),intent(in) :: bands
1074 : :
1075 : : !Local variables-------------------------------
1076 : : !scalars
1077 : : integer :: ik,ib,spin
1078 : : real(dp) :: wtk,occ_bks
1079 : : ! *************************************************************************
1080 : :
1081 : 84 : ex_energy = zero
1082 : :
1083 [ + + ]: 170 : do spin=1,sigma%nsppol
1084 [ + + ]: 638 : do ik=1,sigma%nkibz
1085 : 468 : wtk = kmesh%wt(ik)
1086 [ + + ]: 5224 : do ib=sigma%b1gw,sigma%b2gw
1087 : 4670 : occ_bks = bands%occ(ib,ik,spin)
1088 [ + + ]: 5138 : if (sigma%nsig_ab == 1) then
1089 : 4628 : ex_energy = ex_energy + half * occ_bks * wtk * sigma%sigxme(ib,ik,spin)
1090 : : else
1091 [ + + ]: 210 : ex_energy = ex_energy + half * occ_bks * wtk * SUM(sigma%sigxme(ib,ik,:))
1092 : : end if
1093 : : end do
1094 : : end do
1095 : : end do
1096 : :
1097 : 84 : end function sigma_get_exene
1098 : : !!***
1099 : :
1100 : : !----------------------------------------------------------------------
1101 : :
1102 : : !!****f* m_sigma/sigma_get_excene
1103 : : !! NAME
1104 : : !! sigma_get_excene
1105 : : !!
1106 : : !! FUNCTION
1107 : : !! Compute exchange correlation energy using MBB (nat. orb. functional approx.).
1108 : : !!
1109 : : !! INPUTS
1110 : : !! sigma<sigma_t>=Sigma results
1111 : : !! kmesh<kmesh_t>=BZ sampling.
1112 : : !! bands<band_t>=Bands with occupation factors
1113 : : !!
1114 : : !! SOURCE
1115 : :
1116 : 4 : real(dp) pure function sigma_get_excene(sigma, kmesh, bands) result(exc_energy)
1117 : :
1118 : : !Arguments ------------------------------------
1119 : : class(sigma_t),intent(in) :: sigma
1120 : : type(kmesh_t),intent(in) :: kmesh
1121 : : type(ebands_t),intent(in) :: bands
1122 : :
1123 : : !Local variables-------------------------------
1124 : : !scalars
1125 : : integer :: ik, ib, spin
1126 : : real(dp) :: wtk, occ_bks
1127 : : ! *************************************************************************
1128 : :
1129 : 4 : exc_energy = zero
1130 : :
1131 [ + + ]: 8 : do spin=1,sigma%nsppol
1132 [ + + ]: 32 : do ik=1,sigma%nkibz
1133 : 24 : wtk = kmesh%wt(ik)
1134 [ + + ]: 220 : do ib=sigma%b1gw,sigma%b2gw
1135 : 192 : occ_bks = bands%occ(ib,ik,spin)
1136 [ + - ]: 216 : if (sigma%nsig_ab==1) then
1137 [ + - ]: 192 : if (sigma%nsppol==1) then
1138 : 192 : exc_energy = exc_energy + sqrt( abs( half * occ_bks ) ) * wtk * sigma%sigxcnofme(ib,ik,spin) ! 2*sqrt(occ_i), occ in [0,2] -> [0,1].
1139 : : else
1140 : 0 : exc_energy = exc_energy + half * sqrt( abs( occ_bks ) ) * wtk * sigma%sigxcnofme(ib,ik,spin) ! 2*sqrt(occ_i), occ in [0,1] -> [0,1].
1141 : : end if
1142 : : else
1143 [ # # ]: 0 : exc_energy = exc_energy + half * sqrt( abs( occ_bks ) ) * wtk * SUM(sigma%sigxcnofme(ib,ik,:)) ! 2*sqrt(occ_i), occ in [0,1].
1144 : : end if
1145 : : end do
1146 : : end do
1147 : : end do
1148 : :
1149 : 4 : end function sigma_get_excene
1150 : : !!***
1151 : :
1152 : : !----------------------------------------------------------------------
1153 : :
1154 : : !!****f* ABINIT/sigma_get_haene
1155 : : !! NAME
1156 : : !! sigma_get_haene
1157 : : !!
1158 : : !! FUNCTION
1159 : : !! Compute the Hartree energy
1160 : : !!
1161 : : !! INPUTS
1162 : : !! Kmesh <kmesh_t>=Structure describing the k-point sampling.
1163 : : !! bands=<ebands_t>=Datatype gathering info on the QP energies (KS if one shot)
1164 : : !! eig(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=KS or QP energies for k-points, bands and spin
1165 : : !! occ(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=occupation numbers, for each k point in IBZ, each band and spin
1166 : : !! Mels
1167 : : !! %vhartr=matrix elements of $v_H$.
1168 : : !!
1169 : : !! OUTPUT
1170 : : !! Compute the Hartree energy on eh_energy
1171 : : !!
1172 : : !! SOURCE
1173 : :
1174 : 4 : real(dp) pure function sigma_get_haene(sigma, Mels, kmesh, bands) result(eh_energy)
1175 : :
1176 : : !Arguments ------------------------------------
1177 : : !scalars
1178 : : class(sigma_t),intent(in) :: sigma
1179 : : type(kmesh_t),intent(in) :: kmesh
1180 : : type(ebands_t),intent(in) :: bands
1181 : : type(melements_t),intent(in) :: Mels
1182 : :
1183 : : !Local variables-------------------------------
1184 : : integer :: ik,ib,spin
1185 : : real(dp) :: wtk,occ_bks
1186 : : ! *************************************************************************
1187 : :
1188 : 4 : eh_energy = zero
1189 : :
1190 [ + + ]: 8 : do spin=1,sigma%nsppol
1191 [ + + ]: 32 : do ik=1,sigma%nkibz
1192 : 24 : wtk = kmesh%wt(ik)
1193 [ + + ]: 220 : do ib=sigma%b1gw,sigma%b2gw
1194 : 192 : occ_bks = bands%occ(ib,ik,spin)
1195 [ + - ]: 216 : if (sigma%nsig_ab == 1) then ! Only closed-shell restricted is programed
1196 : 192 : eh_energy = eh_energy + occ_bks * wtk * Mels%vhartree(ib,ib,ik,spin)
1197 : : end if
1198 : : end do
1199 : : end do
1200 : : end do
1201 : :
1202 : 4 : eh_energy = half * eh_energy
1203 : :
1204 : 4 : end function sigma_get_haene
1205 : : !!***
1206 : :
1207 : : !----------------------------------------------------------------------
1208 : :
1209 : : !!****f* ABINIT/sigma_get_kiene
1210 : : !! NAME
1211 : : !! sigma_get_kiene
1212 : : !!
1213 : : !! FUNCTION
1214 : : !! Compute the kinetic energy
1215 : : !!
1216 : : !! INPUTS
1217 : : !! Kmesh <kmesh_t>=Structure describing the k-point sampling.
1218 : : !! bands=<ebands_t>=Datatype gathering info on the QP energies (KS if one shot)
1219 : : !! eig(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=KS or QP energies for k-points, bands and spin
1220 : : !! occ(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=occupation numbers, for each k point in IBZ, each band and spin
1221 : : !! Mels
1222 : : !! %kinetic=matrix elements of $T$.
1223 : : !!
1224 : : !! OUTPUT
1225 : : !! Compute the kinetic energy on ek_energy
1226 : : !!
1227 : : !! SOURCE
1228 : :
1229 : 4 : real(dp) pure function sigma_get_kiene(sigma, Mels, kmesh, bands) result(ek_energy)
1230 : :
1231 : : !Arguments ------------------------------------
1232 : : class(sigma_t),intent(in) :: sigma
1233 : : type(kmesh_t),intent(in) :: kmesh
1234 : : type(ebands_t),intent(in) :: bands
1235 : : type(melements_t),intent(in) :: Mels
1236 : :
1237 : : !Local variables-------------------------------
1238 : : !scalars
1239 : : integer :: ik, ib, spin
1240 : : real(dp) :: wtk, occ_bks
1241 : : ! *************************************************************************
1242 : :
1243 : 4 : ek_energy = zero
1244 : :
1245 [ + + ]: 8 : do spin=1,sigma%nsppol
1246 [ + + ]: 32 : do ik=1,sigma%nkibz
1247 : 24 : wtk = kmesh%wt(ik)
1248 [ + + ]: 220 : do ib=sigma%b1gw,sigma%b2gw
1249 : 192 : occ_bks = bands%occ(ib,ik,spin)
1250 [ + - ]: 216 : if (sigma%nsig_ab==1) then ! Only closed-shell restricted is programed
1251 : 192 : ek_energy = ek_energy + occ_bks * wtk * Mels%kinetic(ib,ib,ik,spin)
1252 : : end if
1253 : : end do
1254 : : end do
1255 : : end do
1256 : :
1257 : 4 : end function sigma_get_kiene
1258 : : !!***
1259 : :
1260 : : !----------------------------------------------------------------------
1261 : :
1262 : : !!****f* m_sigma/find_wpoles_for_cd
1263 : : !! NAME
1264 : : !! find_wpoles_for_cd
1265 : : !!
1266 : : !! FUNCTION
1267 : : !! Find the max frequency needed to account for all the poles of the GW self-energy
1268 : : !! in the case of contour deformation technique.
1269 : : !!
1270 : : !! INPUTS
1271 : : !! Sigp=sigparams_t
1272 : : !!
1273 : : !! OUTPUT
1274 : : !! omega_max
1275 : : !!
1276 : : !! SOURCE
1277 : :
1278 : : subroutine find_wpoles_for_cd(Sigp, sigma, Kmesh, ebands, omega_max)
1279 : :
1280 : : !Arguments ------------------------------------
1281 : : !scalars
1282 : : class(sigparams_t),intent(in) :: Sigp
1283 : : type(sigma_t),intent(in) :: sigma
1284 : : type(ebands_t),intent(in) :: ebands
1285 : : type(kmesh_t),intent(in) :: Kmesh
1286 : : real(dp),intent(out) :: omega_max
1287 : :
1288 : : !Local variables-------------------------------
1289 : : !scalars
1290 : : integer :: spin,ik_ibz,band_gr,bgw_start,bgw_stop,io,ioe0j
1291 : : integer :: ikgw,ikgw_ibz,ikgw_bz,band_gw,nomega_tot
1292 : : real(dp) :: e_green,e_screen,theta_mu_minus_e0i,e_qp,fact_sp
1293 : : !character(len=500) :: msg
1294 : : !arrays
1295 : : real(dp),allocatable :: omegame0i(:)
1296 : : ! *************************************************************************
1297 : :
1298 : : omega_max = smallest_real
1299 : : !
1300 : : ! === Normalization of theta_mu_minus_e0i ===
1301 : : ! * If nsppol==2, qp_occ $\in [0,1]$
1302 : : fact_sp=one
1303 : : if (ebands%nsppol==1) then
1304 : : fact_sp=half; if (ebands%nspinor==2) fact_sp=one
1305 : : end if
1306 : : !
1307 : : ! Total number of frequencies for sigma (Spectral function + mesh for the derivative).
1308 : : nomega_tot=sigma%nomega_r+sigma%nomega4sd
1309 : : ABI_MALLOC(omegame0i,(nomega_tot))
1310 : :
1311 : : ioe0j=sigma%nomega4sd/2+1
1312 : : !
1313 : : ! Loop over bands used to construct the Green function.
1314 : : do spin=1,ebands%nsppol
1315 : : do ik_ibz=1,ebands%nkpt
1316 : : do band_gr=1,ebands%nband(ik_ibz+(spin-1)*ebands%nkpt)
1317 : : e_green = ebands%eig(band_gr,ik_ibz,spin)
1318 : : theta_mu_minus_e0i= ebands%occ(band_gr,ik_ibz,spin)*fact_sp
1319 : : !
1320 : : ! Loop over GW states.
1321 : : do ikgw=1,Sigp%nkptgw
1322 : : bgw_start=Sigp%minbnd(ikgw,spin)
1323 : : bgw_stop =Sigp%minbnd(ikgw,spin)
1324 : : ikgw_bz =Sigp%kptgw2bz(ikgw_bz)
1325 : : ikgw_ibz =Kmesh%tab(ikgw_bz)
1326 : :
1327 : : do band_gw=bgw_start,bgw_stop
1328 : : e_qp = ebands%eig(band_gw,ikgw_ibz,spin)
1329 : : !
1330 : : ! Get frequencies $\omega$-\epsilon_in$ to evaluate $d\Sigma/dE$, note the spin
1331 : : ! subtract e_KS since we have stored e_KS+ Delta \omega in sigma%omega4sd, not required for AC
1332 : : if (sigma%nomega_r>0) omegame0i(1:sigma%nomega_r)=DBLE(Sigp%omega_r(1:sigma%nomega_r))-e_green
1333 : : do io=sigma%nomega_r+1,nomega_tot
1334 : : !omegame0i(io)=DBLE(sigma%omega4sd(band_gw,ikgw_ibz,io-sigma%nomega_r,spin)) - e_green
1335 : : !sigma%omega4sd(jb,ik_ibz,io,spin)=sigma%egw(jb,ik_ibz,spin)+Sigp%deltae*(io-ioe0j)
1336 : : omegame0i(io) = e_qp + Sigp%deltae*(io-ioe0j) - e_green
1337 : : end do
1338 : :
1339 : : do io=1,nomega_tot
1340 : : e_screen = ABS(omegame0i(io))
1341 : : if (omegame0i(io)>tol12) then
1342 : : !ket(spadc+ig,ios)=ket(spadc+ig,ios)+ct*(one-theta_mu_minus_e0i)
1343 : : if ( (one-theta_mu_minus_e0i) > tol12 ) omega_max = MAX(omega_max, e_screen)
1344 : : end if
1345 : : if (omegame0i(io)<-tol12) then
1346 : : !ket(spadc+ig,ios)=ket(spadc+ig,ios)-ct*theta_mu_minus_e0i
1347 : : if ( theta_mu_minus_e0i > tol12) omega_max = MAX(omega_max, e_screen)
1348 : : end if
1349 : : end do
1350 : :
1351 : : end do
1352 : : end do
1353 : : !
1354 : : end do
1355 : : end do
1356 : : end do
1357 : :
1358 : : ABI_FREE(omegame0i)
1359 : :
1360 : : end subroutine find_wpoles_for_cd
1361 : : !!***
1362 : :
1363 : : !----------------------------------------------------------------------
1364 : :
1365 : : !!****f* m_sigma/sigma_ncwrite
1366 : : !! NAME
1367 : : !! sigma_ncwrite
1368 : : !!
1369 : : !! FUNCTION
1370 : : !! Save the data stored in the sigma_t data type on a NETCDF file.
1371 : : !!
1372 : : !! INPUTS
1373 : : !! filename
1374 : : !!
1375 : : !! OUTPUT
1376 : : !!
1377 : : !! SOURCE
1378 : :
1379 : 168 : integer function sigma_ncwrite(sigma, Sigp, epsm1, dtset, ncid) result (ncerr)
1380 : :
1381 : : !Arguments ------------------------------------
1382 : : !scalars
1383 : : class(sigma_t),target,intent(in) :: sigma
1384 : : class(sigparams_t),target,intent(in) :: Sigp
1385 : : integer,intent(in) :: ncid
1386 : : type(epsm1_t),target,intent(in) :: epsm1
1387 : : type(dataset_type),intent(in) :: dtset
1388 : :
1389 : : !Local variables ---------------------------------------
1390 : : !scalars
1391 : : integer :: nbgw,ndim_sig,b1gw,b2gw,cplex
1392 : : !character(len=500) :: msg
1393 : : !arrays
1394 : 168 : real(dp),allocatable :: rdata2(:,:),rdata4(:,:,:,:),rdata5(:,:,:,:,:)
1395 : : ! *************************************************************************
1396 : :
1397 : : !@sigma_t
1398 : 168 : cplex=2; b1gw=sigma%b1gw; b2gw=sigma%b2gw; nbgw=b2gw-b1gw+1
1399 : 168 : ndim_sig=sigma%nsppol*sigma%nsig_ab
1400 : :
1401 : : ncerr = nctk_def_dims(ncid, [&
1402 : : nctkdim_t("cplex", cplex), nctkdim_t("b1gw", sigma%b1gw), nctkdim_t("b2gw", sigma%b2gw),&
1403 : : nctkdim_t("nbgw", nbgw), nctkdim_t("nkptgw", sigma%nkptgw), nctkdim_t("ndim_sig", ndim_sig), &
1404 : : nctkdim_t("nomega4sd", sigma%nomega4sd), nctkdim_t("nsig_ab", sigma%nsig_ab), &
1405 : : nctkdim_t("scr_nomega", epsm1%nomega) &
1406 [ + + ]: 1680 : ], defmode=.True.)
1407 [ - + ]: 168 : NCF_CHECK(ncerr)
1408 : :
1409 : : ! No. of real frequencies, might be zero.
1410 [ + + ]: 168 : if (sigma%nomega_r > 0) then
1411 [ - + ]: 5 : NCF_CHECK(nctk_def_dims(ncid, nctkdim_t("nomega_r", sigma%nomega_r)))
1412 : : end if
1413 : :
1414 : : ! No. of imaginary frequencies, might be zero.
1415 [ + + ]: 168 : if (sigma%nomega_i > 0) then
1416 [ - + ]: 4 : NCF_CHECK(nctk_def_dims(ncid, nctkdim_t("nomega_i", sigma%nomega_i)))
1417 : : end if
1418 : :
1419 : : ! =======================
1420 : : ! == Define variables ===
1421 : : ! =======================
1422 : : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: &
1423 : : 'sigma_nband', 'scr_nband', 'gwcalctyp', 'usepawu', "nfreqre", "nfreqim", "nfreqim_conv", &
1424 [ + + ]: 1680 : "gwcomp", "ppmodel"])
1425 [ - + ]: 168 : NCF_CHECK(ncerr)
1426 : :
1427 : : ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: &
1428 [ + + ]: 1512 : 'ecutwfn', 'ecuteps', 'ecutsigx', 'omegasrdmax', 'deltae', 'omegasrmax', 'scissor_ene', "gwencomp"])
1429 [ - + ]: 168 : NCF_CHECK(ncerr)
1430 : :
1431 : : ! TODO: Decrease size of file: Remove arrays whose size scale as mband ** 2
1432 : : ! especially those that are not commonly used e.g. hhartree.
1433 : : ncerr = nctk_def_arrays(ncid, [ &
1434 : : nctkarr_t("kptgw", "dp", "number_of_reduced_dimensions, nkptgw"),&
1435 : : nctkarr_t("minbnd", "i", "nkptgw, number_of_spins"),&
1436 : : nctkarr_t("maxbnd", "i", "nkptgw, number_of_spins"), &
1437 : : nctkarr_t('degwgap', "dp", 'number_of_kpoints, number_of_spins'),&
1438 : : nctkarr_t('egwgap', "dp", 'number_of_kpoints, number_of_spins'),&
1439 : : nctkarr_t('en_qp_diago', "dp",'max_number_of_states, number_of_kpoints, number_of_spins'),&
1440 : : nctkarr_t('e0', "dp", 'max_number_of_states, number_of_kpoints, number_of_spins'),&
1441 : : nctkarr_t('e0gap', "dp", 'number_of_kpoints, number_of_spins'),&
1442 : : nctkarr_t('sigxme', "dp", 'nbgw, number_of_kpoints, ndim_sig'),&
1443 : : nctkarr_t('vxcme', "dp", 'nbgw, number_of_kpoints, ndim_sig'),&
1444 : : nctkarr_t('degw', "dp", 'cplex, nbgw, number_of_kpoints, number_of_spins'),&
1445 : : nctkarr_t('dsigmee0', "dp", 'cplex, nbgw, number_of_kpoints, ndim_sig'),&
1446 : : nctkarr_t('egw', "dp",'cplex, max_number_of_states, number_of_kpoints, number_of_spins'),&
1447 : : nctkarr_t('hhartree', "dp",'cplex, nbgw, nbgw, number_of_kpoints, ndim_sig'),&
1448 : : nctkarr_t('sigmee', "dp", 'cplex, nbgw, number_of_kpoints, ndim_sig'),&
1449 : : nctkarr_t('sigcmee0', "dp",'cplex, nbgw, number_of_kpoints, ndim_sig'),&
1450 : : nctkarr_t('sigcme4sd', "dp",'cplex, nbgw, number_of_kpoints, nomega4sd, ndim_sig'),&
1451 : : nctkarr_t('sigxcme4sd', "dp", 'cplex, nbgw, number_of_kpoints, nomega4sd, ndim_sig'),&
1452 : : nctkarr_t('ze0',"dp", 'cplex, nbgw, number_of_kpoints, number_of_spins'),&
1453 : : nctkarr_t('omega4sd', "dp", 'cplex, nbgw, number_of_kpoints, nomega4sd, number_of_spins') &
1454 [ + + ]: 3528 : ])
1455 [ - + ]: 168 : NCF_CHECK(ncerr)
1456 : :
1457 [ + + ]: 168 : if (sigma%needs_eigvec_qp) then
1458 : : ncerr = nctk_def_arrays(ncid, [ &
1459 [ + + ]: 110 : nctkarr_t('eigvec_qp', "dp",'cplex, max_number_of_states, max_number_of_states, number_of_kpoints, number_of_spins')])
1460 [ - + ]: 55 : NCF_CHECK(ncerr)
1461 : : end if
1462 : :
1463 [ + - ]: 168 : if (sigma%usepawu == 0) then
1464 : 168 : ncerr = nctk_def_arrays(ncid, nctkarr_t("vUme", "dp", 'nbgw, number_of_kpoints, ndim_sig'))
1465 [ - + ]: 168 : NCF_CHECK(ncerr)
1466 : : end if
1467 : :
1468 [ + + ]: 168 : if (epsm1%nomega > 0) then
1469 : 133 : ncerr = nctk_def_arrays(ncid, nctkarr_t('scr_omega', "dp", 'cplex, scr_nomega'))
1470 [ - + ]: 133 : NCF_CHECK(ncerr)
1471 [ - + + - : 399 : ABI_MALLOC(rdata2, (2, epsm1%nomega))
+ - + - -
+ - + ]
1472 [ + - - + : 133 : rdata2 = c2r(epsm1%omega)
- - ]
1473 [ + + + + : 1984 : NCF_CHECK(nf90_put_var(ncid, vid('scr_omega'), rdata2 * Ha_eV))
- + - - -
- ]
1474 [ - + ]: 133 : ABI_FREE(rdata2)
1475 : : end if
1476 : :
1477 [ + + ]: 168 : if (sigma%nomega_r > 0) then
1478 : : ncerr = nctk_def_arrays(ncid, [&
1479 : : nctkarr_t('omega_r', "dp", "nomega_r"),&
1480 : : nctkarr_t('sigcme', "dp", 'cplex, nbgw, number_of_kpoints, nomega_r, ndim_sig'),&
1481 [ + + ]: 20 : nctkarr_t('sigxcme', "dp", 'cplex, nbgw, number_of_kpoints, nomega_r, ndim_sig')])
1482 [ - + ]: 5 : NCF_CHECK(ncerr)
1483 : : end if
1484 : :
1485 [ + + ]: 168 : if (sigma%nomega_i > 0) then
1486 : : ncerr = nctk_def_arrays(ncid, [&
1487 : : nctkarr_t('sigxcmesi', "dp", 'cplex, nbgw, number_of_kpoints, nomega_i, ndim_sig'),&
1488 : : nctkarr_t('sigcmesi', "dp",'cplex, nbgw, number_of_kpoints, nomega_i, ndim_sig'),&
1489 [ + + ]: 16 : nctkarr_t('omega_i', "dp", 'cplex, nomega_i')])
1490 [ - + ]: 4 : NCF_CHECK(ncerr)
1491 : : end if
1492 : :
1493 [ + + ]: 168 : if (allocated(sigma%m_ks_to_qp)) then
1494 : : ncerr = nctk_def_arrays(ncid, [nctkarr_t('m_ks_to_qp', "dp", &
1495 [ + + ]: 110 : "cplex, max_number_of_states, max_number_of_states, number_of_kpoints, number_of_spins")])
1496 [ - + ]: 55 : NCF_CHECK(ncerr)
1497 : : end if
1498 : :
1499 : : ! =====================
1500 : : ! === Start writing ===
1501 : : ! =====================
1502 [ - + ]: 168 : NCF_CHECK(nctk_set_datamode(ncid))
1503 [ - + ]: 168 : NCF_CHECK(nf90_put_var(ncid, vid('ecutwfn'), Sigp%ecutwfn))
1504 [ - + ]: 168 : NCF_CHECK(nf90_put_var(ncid, vid('ecuteps'), Sigp%ecuteps))
1505 [ - + ]: 168 : NCF_CHECK(nf90_put_var(ncid, vid('ecutsigx'), Sigp%ecutsigx))
1506 [ - + ]: 168 : NCF_CHECK(nf90_put_var(ncid, vid('sigma_nband'), Sigp%nbnds))
1507 [ - + ]: 168 : NCF_CHECK(nf90_put_var(ncid, vid('scr_nband'), epsm1%Hscr%nbnds_used))
1508 [ - + ]: 168 : NCF_CHECK(nf90_put_var(ncid, vid('gwcalctyp'), sigma%gwcalctyp))
1509 [ - + ]: 168 : NCF_CHECK(nf90_put_var(ncid, vid('usepawu'), sigma%usepawu))
1510 [ - + ]: 168 : NCF_CHECK(nf90_put_var(ncid, vid('nfreqre'), epsm1%nomega_r))
1511 [ - + ]: 168 : NCF_CHECK(nf90_put_var(ncid, vid('nfreqim'), epsm1%nomega_i))
1512 [ - + ]: 168 : NCF_CHECK(nf90_put_var(ncid, vid('nfreqim_conv'), epsm1%nomega_i_conv))
1513 : : ! Input variables (values given in the input file, not the quantities derived from them in the code).
1514 [ - + ]: 168 : NCF_CHECK(nf90_put_var(ncid, vid('gwcomp'), dtset%gwcomp))
1515 [ - + ]: 168 : NCF_CHECK(nf90_put_var(ncid, vid('ppmodel'), dtset%ppmodel))
1516 [ - + ]: 168 : NCF_CHECK(nf90_put_var(ncid, vid('gwencomp'), dtset%gwencomp))
1517 [ - + ]: 168 : NCF_CHECK(nf90_put_var(ncid, vid('kptgw'), sigma%kptgw))
1518 [ - + ]: 168 : NCF_CHECK(nf90_put_var(ncid, vid('minbnd'), sigma%minbnd))
1519 [ - + ]: 168 : NCF_CHECK(nf90_put_var(ncid, vid('maxbnd'),sigma%maxbnd))
1520 [ - + ]: 168 : NCF_CHECK(nf90_put_var(ncid, vid('omegasrdmax'), sigma%maxomega4sd * Ha_eV))
1521 [ - + ]: 168 : NCF_CHECK(nf90_put_var(ncid, vid('deltae'), sigma%deltae * Ha_eV))
1522 [ - + ]: 168 : NCF_CHECK(nf90_put_var(ncid, vid('omegasrmax'), sigma%maxomega_r * Ha_eV))
1523 [ - + ]: 168 : NCF_CHECK(nf90_put_var(ncid, vid('scissor_ene'), sigma%scissor_ene * Ha_eV))
1524 [ + + + + : 1338 : NCF_CHECK(nf90_put_var(ncid, vid('degwgap'), sigma%degwgap * Ha_eV))
- + - - -
- ]
1525 [ + + + + : 1338 : NCF_CHECK(nf90_put_var(ncid, vid('egwgap'), sigma%egwgap * Ha_eV))
- + - - -
- ]
1526 [ + + + + : 25900 : NCF_CHECK(nf90_put_var(ncid, vid('en_qp_diago'), sigma%en_qp_diago * Ha_eV))
+ + - + -
- - - -
- ]
1527 [ + + + + : 25900 : NCF_CHECK(nf90_put_var(ncid, vid('e0'), sigma%e0 * Ha_eV))
+ + - + -
- - - -
- ]
1528 [ + + + + : 1338 : NCF_CHECK(nf90_put_var(ncid, vid('e0gap'), sigma%e0gap * Ha_eV))
- + - - -
- ]
1529 : :
1530 [ + + ]: 168 : if (sigma%nomega_r > 0) then
1531 [ + + - + : 813 : NCF_CHECK(nf90_put_var(ncid, vid('omega_r'), sigma%omega_r * Ha_eV))
- - ]
1532 : : end if
1533 : :
1534 [ + + + + : 9760 : NCF_CHECK(nf90_put_var(ncid, vid('sigxme'), sigma%sigxme * Ha_eV))
+ + - + -
- - - -
- ]
1535 [ + + + + : 9760 : NCF_CHECK(nf90_put_var(ncid, vid('vxcme'), sigma%vxcme * Ha_eV))
+ + - + -
- - - -
- ]
1536 [ + + + + : 9760 : NCF_CHECK(nf90_put_var(ncid, vid('vUme'), sigma%vUme * Ha_eV))
+ + - + -
- - - -
- ]
1537 : :
1538 : : ! Have to transfer complex arrays
1539 [ - + + - : 840 : ABI_MALLOC(rdata4,(cplex, b1gw:b2gw, sigma%nkibz, sigma%nsppol))
- + + - -
+ + - + -
+ - - + -
+ - + ]
1540 [ + - + - : 168 : rdata4=c2r(sigma%degw)
+ - - + -
- ]
1541 [ + + + + : 22644 : NCF_CHECK(nf90_put_var(ncid, vid('degw'), rdata4 * Ha_eV))
+ + + + -
+ - - - -
- - - - ]
1542 [ - + ]: 168 : ABI_FREE(rdata4)
1543 : :
1544 [ - + + - : 840 : ABI_MALLOC(rdata4,(cplex, b1gw:b2gw, sigma%nkibz, sigma%nsppol*sigma%nsig_ab))
- + + - -
+ + - + -
+ - - + -
+ ]
1545 [ + - + - : 168 : rdata4 = c2r(sigma%dsigmee0)
+ - - + -
- ]
1546 [ - + ]: 168 : NCF_CHECK(nf90_put_var(ncid, vid('dsigmee0'), rdata4))
1547 [ - + ]: 168 : ABI_FREE(rdata4)
1548 : :
1549 [ - + + - : 840 : ABI_MALLOC(rdata4, (cplex, sigma%nbnds, sigma%nkibz, sigma%nsppol))
- + + - -
+ + - + -
+ - - + -
+ ]
1550 [ + - + - : 168 : rdata4 = c2r(sigma%egw)
+ - - + -
- ]
1551 [ + + + + : 75024 : NCF_CHECK(nf90_put_var(ncid, vid('egw'), rdata4 *Ha_eV))
+ + + + -
+ - - - -
- - - - ]
1552 [ - + ]: 168 : ABI_FREE(rdata4)
1553 : :
1554 [ + + ]: 168 : if (sigma%needs_eigvec_qp) then
1555 [ - + + - : 330 : ABI_MALLOC(rdata5, (cplex, sigma%nbnds, sigma%nbnds, sigma%nkibz, sigma%nsppol))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
1556 [ + - + - : 55 : rdata5 = c2r(sigma%eigvec_qp)
+ - + - -
+ - - ]
1557 [ - + ]: 55 : NCF_CHECK(nf90_put_var(ncid, vid('eigvec_qp'), rdata5))
1558 [ - + ]: 55 : ABI_FREE(rdata5)
1559 : : end if
1560 : :
1561 [ - + + - : 1008 : ABI_MALLOC(rdata5,(cplex, nbgw, nbgw, sigma%nkibz, sigma%nsppol * sigma%nsig_ab))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
1562 [ + - + - : 168 : rdata5 = c2r(sigma%hhartree)
+ - + - -
+ - - ]
1563 [ + + + + : 252766 : NCF_CHECK(nf90_put_var(ncid, vid('hhartree'), rdata5 * Ha_eV))
+ + + + +
+ - + - -
- - - - -
- - - ]
1564 [ - + ]: 168 : ABI_FREE(rdata5)
1565 : :
1566 [ + + ]: 168 : if (sigma%nomega_r > 0) then
1567 [ - + + - : 30 : ABI_MALLOC(rdata5,(cplex, nbgw, sigma%nkibz, sigma%nomega_r, sigma%nsppol*sigma%nsig_ab))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
1568 [ + - + - : 5 : rdata5 = c2r(sigma%sigcme)
+ - + - -
+ - - ]
1569 [ + + + + : 157566 : NCF_CHECK(nf90_put_var(ncid, vid('sigcme'), rdata5 * Ha_eV))
+ + + + +
+ - + - -
- - - - -
- - - ]
1570 [ - + ]: 5 : ABI_FREE(rdata5)
1571 : : end if
1572 : :
1573 [ - + + - : 840 : ABI_MALLOC(rdata4, (cplex, nbgw, sigma%nkibz, sigma%nsppol*sigma%nsig_ab))
- + + - -
+ + - + -
+ - - + -
+ ]
1574 [ + - + - : 168 : rdata4 = c2r(sigma%sigmee)
+ - - + -
- ]
1575 [ + + + + : 26376 : NCF_CHECK(nf90_put_var(ncid, vid('sigmee'), rdata4 * Ha_eV))
+ + + + -
+ - - - -
- - - - ]
1576 [ - + ]: 168 : ABI_FREE(rdata4)
1577 : :
1578 [ - + + - : 840 : ABI_MALLOC(rdata4, (cplex, nbgw, sigma%nkibz, sigma%nsppol*sigma%nsig_ab))
- + + - -
+ + - + -
+ - - + -
+ ]
1579 [ + - + - : 168 : rdata4 = c2r(sigma%sigcmee0)
+ - - + -
- ]
1580 [ + + + + : 26376 : NCF_CHECK(nf90_put_var(ncid, vid('sigcmee0'), rdata4 * Ha_eV))
+ + + + -
+ - - - -
- - - - ]
1581 [ - + ]: 168 : ABI_FREE(rdata4)
1582 : :
1583 [ + + ]: 168 : if (sigma%nomega_i > 0) then
1584 [ - + + - : 24 : ABI_MALLOC(rdata5, (cplex, nbgw, sigma%nkibz, sigma%nomega_i, sigma%nsppol*sigma%nsig_ab))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
1585 [ + - + - : 4 : rdata5 = c2r(sigma%sigcmesi)
+ - + - -
+ - - ]
1586 [ + + + + : 7558 : NCF_CHECK(nf90_put_var(ncid, vid('sigcmesi'), rdata5*Ha_eV))
+ + + + +
+ - + - -
- - - - -
- - - ]
1587 [ - + ]: 4 : ABI_FREE(rdata5)
1588 : : end if
1589 : :
1590 [ - + + - : 1008 : ABI_MALLOC(rdata5, (cplex, nbgw, sigma%nkibz, sigma%nomega4sd, sigma%nsppol*sigma%nsig_ab))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
1591 [ + - + - : 168 : rdata5 = c2r(sigma%sigcme4sd)
+ - + - -
+ - - ]
1592 [ + + + + : 125177 : NCF_CHECK(nf90_put_var(ncid, vid('sigcme4sd'), rdata5 * Ha_eV))
+ + + + +
+ - + - -
- - - - -
- - - ]
1593 [ - + ]: 168 : ABI_FREE(rdata5)
1594 : :
1595 [ + + ]: 168 : if (sigma%nomega_r > 0) then
1596 [ - + + - : 30 : ABI_MALLOC(rdata5,(cplex, nbgw, sigma%nkibz, sigma%nomega_r, sigma%nsppol*sigma%nsig_ab))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
1597 [ + - + - : 5 : rdata5 = c2r(sigma%sigxcme)
+ - + - -
+ - - ]
1598 [ + + + + : 157566 : NCF_CHECK(nf90_put_var(ncid, vid('sigxcme'), rdata5 * Ha_eV))
+ + + + +
+ - + - -
- - - - -
- - - ]
1599 [ - + ]: 5 : ABI_FREE(rdata5)
1600 : : end if
1601 : :
1602 [ + + ]: 168 : if (sigma%nomega_i > 0) then
1603 [ - + + - : 24 : ABI_MALLOC(rdata5,(cplex, nbgw, sigma%nkibz, sigma%nomega_i, sigma%nsppol*sigma%nsig_ab))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
1604 [ + - + - : 4 : rdata5 = c2r(sigma%sigxcmesi)
+ - + - -
+ - - ]
1605 [ + + + + : 7558 : NCF_CHECK(nf90_put_var(ncid, vid('sigxcmesi'), rdata5 * Ha_eV))
+ + + + +
+ - + - -
- - - - -
- - - ]
1606 [ - + ]: 4 : ABI_FREE(rdata5)
1607 : : end if
1608 : :
1609 [ + + ]: 168 : if (allocated(sigma%m_ks_to_qp)) then
1610 [ - + + - : 330 : ABI_MALLOC(rdata5,(cplex, sigma%nbnds, sigma%nbnds, sigma%nkibz, sigma%nsppol))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
1611 [ + - + - : 55 : rdata5 = c2r(sigma%m_ks_to_qp)
+ - + - -
+ - - ]
1612 [ - + ]: 55 : NCF_CHECK(nf90_put_var(ncid, vid('m_ks_to_qp'), rdata5))
1613 [ - + ]: 55 : ABI_FREE(rdata5)
1614 : : end if
1615 : :
1616 [ - + + - : 1008 : ABI_MALLOC(rdata5, (cplex, nbgw, sigma%nkibz, sigma%nomega4sd, sigma%nsppol*sigma%nsig_ab))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
1617 [ + - + - : 168 : rdata5 = c2r(sigma%sigxcme4sd)
+ - + - -
+ - - ]
1618 [ + + + + : 125177 : NCF_CHECK(nf90_put_var(ncid, vid('sigxcme4sd'), rdata5 * Ha_eV))
+ + + + +
+ - + - -
- - - - -
- - - ]
1619 [ - + ]: 168 : ABI_FREE(rdata5)
1620 : :
1621 [ - + + - : 840 : ABI_MALLOC(rdata4, (cplex, nbgw, sigma%nkibz, sigma%nsppol))
- + + - -
+ + - + -
+ - - + -
+ - + ]
1622 [ + - + - : 168 : rdata4 = c2r(sigma%ze0)
+ - - + -
- ]
1623 [ - + ]: 168 : NCF_CHECK(nf90_put_var(ncid, vid('ze0'), rdata4))
1624 [ - + ]: 168 : ABI_FREE(rdata4)
1625 : :
1626 [ + + ]: 168 : if (sigma%nomega_i > 0) then
1627 [ - + + - : 12 : ABI_MALLOC(rdata2, (cplex, sigma%nomega_i))
+ - - + -
+ - + ]
1628 [ + - - + : 4 : rdata2 = c2r(sigma%omega_i)
- - ]
1629 [ + + + + : 154 : NCF_CHECK(nf90_put_var(ncid, vid('omega_i'), rdata2 * Ha_eV))
- + - - -
- ]
1630 [ - + ]: 4 : ABI_FREE(rdata2)
1631 : : end if
1632 : :
1633 [ - + + - : 1008 : ABI_MALLOC(rdata5, (cplex, nbgw, sigma%nkibz, sigma%nomega4sd, sigma%nsppol))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
1634 [ + - + - : 168 : rdata5 = c2r(sigma%omega4sd)
+ - + - -
+ - - ]
1635 [ + + + + : 108944 : NCF_CHECK(nf90_put_var(ncid, vid('omega4sd'), rdata5 * Ha_eV))
+ + + + +
+ - + - -
- - - - -
- - - ]
1636 [ - + - + ]: 168 : ABI_FREE(rdata5)
1637 : :
1638 : :
1639 : : contains
1640 : 6654 : integer function vid(vname)
1641 : : character(len=*),intent(in) :: vname
1642 [ + - + - : 6654 : vid = nctk_idname(ncid, vname)
- - + - -
- + - - -
+ - - - +
- - - + -
- - + - -
- + - - -
+ - - - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - +
- ]
1643 : : end function vid
1644 : :
1645 : : end function sigma_ncwrite
1646 : : !!***
1647 : :
1648 : : !----------------------------------------------------------------------
1649 : :
1650 : : !!****f* m_sigma/sigma_distribute_bks
1651 : : !! NAME
1652 : : !! sigma_distribute_bks
1653 : : !!
1654 : : !! FUNCTION
1655 : : !! Distribute the loop over (b,k,s) used to calculate the self-energy matrix elements
1656 : : !! taking into account the MPI distribution of the wavefunctions and the use of
1657 : : !! symmetries to reduce the BZ sum to an appropriate irreducible wedge.
1658 : : !!
1659 : : !! INPUTS
1660 : : !! nsppol
1661 : : !! can_symmetrize(nsppol)=.TRUE if symmetries can be used to reduce the number of k-points to be summed.
1662 : : !! Kmesh<kmesh_t>
1663 : : !! Qmesh<kmesh_t>
1664 : : !! Ltg_kgw<littlegroup_t>
1665 : : !! Wfd(wfdgw_t)
1666 : : !! mg0(3)
1667 : : !! kptgw(3)
1668 : : !! [bks_mask(Wfd%mband,Kmesh%nbz,nsppol)]
1669 : : !! [got(Wfd%nproc)]=The number of tasks already assigned to the nodes.
1670 : : !! [global]=If true, an MPI global communication is performed such that each node will have the same table. Useful
1671 : : !! if for implementing algorithms in which each node needs to know the global distribution of the tasks, not only
1672 : : !! the task it has to complete. Defaults to .FALSE.
1673 : : !!
1674 : : !! OUTPUT
1675 : : !! my_nbks
1676 : : !! proc_distrb(Wfd%mband,Kmesh%nbz,nsppol)
1677 : : !!
1678 : : !! SIDE EFFECTS
1679 : : !! Wfd%bks_tab
1680 : : !!
1681 : : !! SOURCE
1682 : :
1683 : 1125 : subroutine sigma_distribute_bks(Wfd,Kmesh,Ltg_kgw,Qmesh,nsppol,can_symmetrize,kptgw,mg0,my_nbks,proc_distrb,got,bks_mask,global)
1684 : :
1685 : : !Arguments ------------------------------------
1686 : : !scalars
1687 : : class(wfdgw_t),intent(inout) :: Wfd
1688 : : integer,intent(in) :: nsppol
1689 : : integer,intent(out) :: my_nbks
1690 : : logical,optional,intent(in) :: global
1691 : : type(kmesh_t),intent(in) :: Kmesh,Qmesh
1692 : : type(littlegroup_t),intent(in) :: Ltg_kgw
1693 : : !arrays
1694 : : integer,intent(in) :: mg0(3)
1695 : : integer,optional,intent(inout) :: got(Wfd%nproc)
1696 : : integer,intent(out) :: proc_distrb(Wfd%mband,Kmesh%nbz,nsppol)
1697 : : real(dp),intent(in) :: kptgw(3)
1698 : : logical,intent(in) :: can_symmetrize(Wfd%nsppol)
1699 : : logical,optional,intent(in) :: bks_mask(Wfd%mband,Kmesh%nbz,nsppol)
1700 : :
1701 : : !Local variables-------------------------------
1702 : : !scalars
1703 : : integer :: ierr,ik_bz,ik_ibz,spin,iq_bz,my_nband
1704 : : !character(len=500) :: msg
1705 : : !arrays
1706 : : integer :: g0(3)
1707 : : real(dp) :: kgwmk(3)
1708 : 2240 : integer :: get_more(Wfd%nproc),my_band_list(Wfd%mband)
1709 : 2240 : logical :: bmask(Wfd%mband)
1710 : : !************************************************************************
1711 : :
1712 : 1120 : call wfd%update_bkstab()
1713 : :
1714 [ + + + + : 4336 : get_more=0; if (PRESENT(got)) get_more=got
+ + ]
1715 : :
1716 : : ! Different distribution of tasks depending whether symmetries can be used or not.
1717 [ + + + + : 1010966 : proc_distrb= xmpi_undefined_rank
+ + ]
1718 : :
1719 [ + + ]: 2254 : do spin=1,Wfd%nsppol
1720 : :
1721 [ + + ]: 2254 : if (can_symmetrize(spin)) then
1722 [ + + ]: 14934 : do ik_bz=1,Kmesh%nbz
1723 : 14704 : ik_ibz = Kmesh%tab(ik_bz)
1724 [ + + ]: 58816 : kgwmk = kptgw-Kmesh%bz(:,ik_bz) ! kptgw must be inside the BZ
1725 : 14704 : call findqg0(iq_bz,g0,kgwmk,Qmesh%nbz,Qmesh%bz,mG0) ! <- (mg0=mG0) Identify q_bz and G0 where q_bz+G0=k_gw-k_bz
1726 [ + + ]: 14934 : if (Ltg_kgw%ibzq(iq_bz)==1) then
1727 [ + + + + ]: 173696 : bmask=.FALSE.; bmask(1:Wfd%nband(ik_ibz,spin))=.TRUE.
1728 [ + + + + ]: 45896 : if (PRESENT(bks_mask)) bmask = bks_mask(:,ik_bz,spin)
1729 : 2728 : call wfd%distribute_bands(ik_ibz,spin,my_nband,my_band_list,got=get_more,bmask=bmask)
1730 [ + + + + ]: 47789 : if (my_nband>0) proc_distrb(my_band_list(1:my_nband),ik_bz,spin)=Wfd%my_rank
1731 : : end if
1732 : : end do
1733 : :
1734 : : else
1735 : : ! No symmetries for this spin. Divide the full BZ among procs.
1736 [ + + ]: 45964 : do ik_bz=1,Kmesh%nbz
1737 : 45060 : ik_ibz = Kmesh%tab(ik_bz)
1738 [ + + + + ]: 1167740 : bmask=.FALSE.; bmask(1:Wfd%nband(ik_ibz,spin))=.TRUE.
1739 [ + + + + ]: 371468 : if (PRESENT(bks_mask)) bmask = bks_mask(:,ik_bz,spin)
1740 : 45060 : call wfd%distribute_bands(ik_ibz,spin,my_nband,my_band_list,got=get_more,bmask=bmask)
1741 [ + + + + ]: 300774 : if (my_nband>0) proc_distrb(my_band_list(1:my_nband),ik_bz,spin)=Wfd%my_rank
1742 : : end do
1743 : : end if
1744 : : end do ! spin
1745 : :
1746 [ + + ]: 1120 : if (PRESENT(global)) then
1747 [ + - ]: 460 : if (global) then ! Each node will have the same table so that it will know how the tasks are distributed.
1748 [ + + + + : 453378 : proc_distrb = proc_distrb + 1
+ + ]
1749 [ + + + + : 453378 : where (proc_distrb == xmpi_undefined_rank + 1)
+ + + + ]
1750 : : proc_distrb = 0
1751 : : end where
1752 : 460 : call xmpi_sum(proc_distrb,Wfd%comm,ierr)
1753 [ + + + + : 1360134 : where (proc_distrb == 0)
+ + + + +
+ + + + +
+ + + + +
+ + + ]
1754 : : proc_distrb = xmpi_undefined_rank
1755 : : elsewhere
1756 : : proc_distrb = proc_distrb - 1
1757 : : end where
1758 : : !where (proc_distrb /= xmpi_undefined_rank)
1759 : : ! ltest = (ANY(proc_distrb == (/(ii,ii=0,Wfd%nproc-1)/)))
1760 : : !end where
1761 : : !if (.not.ltest) then
1762 : : ! write(std_out,*)proc_distrb
1763 : : ! ABI_BUG("Bug in the generation of proc_distrb table")
1764 : : !end if
1765 : : end if
1766 : : end if
1767 : :
1768 [ + + + + : 1010966 : my_nbks = COUNT(proc_distrb==Wfd%my_rank)
+ + + + ]
1769 [ + + + + ]: 2432 : if (PRESENT(got)) got=get_more
1770 : :
1771 [ + + + + ]: 1580 : end subroutine sigma_distribute_bks
1772 : : !!***
1773 : :
1774 : : !----------------------------------------------------------------------
1775 : :
1776 [ # # # # : 0 : end module m_sigma
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # #
# ]
1777 : : !!***
|
