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