Line data Source code
1 : !!****m* ABINIT/m_gstore_sigmaph
2 : !! NAME
3 : !! m_gstore_sigmaph
4 : !!
5 : !! FUNCTION
6 : !! Compute (diagonal) matrix elements of the e-ph self-energy (Fan Migdal + Debye Waller).
7 : !! in the KS basis using precomputed e-ph matrix elements.
8 : !! See also m_sigmaph, for a version in which the g-matrix elements are computed on-the-fly.
9 : !!
10 : !! COPYRIGHT
11 : !! Copyright (C) 2008-2025 ABINIT group (MG)
12 : !! This file is distributed under the terms of the
13 : !! GNU General Public License, see ~abinit/COPYING
14 : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : !!
16 : !! SOURCE
17 :
18 : #if defined HAVE_CONFIG_H
19 : #include "config.h"
20 : #endif
21 :
22 : #include "abi_common.h"
23 :
24 : module m_gstore_sigmaph
25 :
26 : use, intrinsic :: iso_c_binding
27 : use defs_basis
28 : use m_abicore
29 : use m_xmpi
30 : use m_errors
31 : use netcdf
32 : use m_nctk
33 : use m_ephtk
34 : use m_sigtk
35 :
36 : use m_io_tools, only : iomode_from_fname
37 : use m_numeric_tools, only : arth, c2r
38 : use m_time, only : cwtime, cwtime_report
39 : use m_fstrings, only : tolower, itoa, ftoa, sjoin, ktoa, ltoa, strcat, replace_ch0, yesno, string_in
40 : use m_special_funcs, only : gaussian
41 : use m_cgtools, only : cg_zgemm, cg_zdotc
42 : use m_kg, only : getph
43 : use defs_datatypes, only : pseudopotential_type
44 : use defs_abitypes, only : mpi_type
45 : use m_hdr, only : hdr_type, fform_from_ext
46 : use m_geometry, only : phdispl_cart2red_nmodes
47 : use m_ebands, only : ebands_t, gaps_t
48 : use m_kpts, only : kpts_timrev_from_kptopt, kpts_map
49 : use m_ioarr, only : read_rhor
50 : use m_fftcore, only : ngfft_seq
51 : use m_getgh1c, only : getgh1c, rf_transgrid_and_pack
52 : use m_ifc, only : ifc_type
53 : use m_dfpt_cgwf, only : stern_t
54 : use m_pawang, only : pawang_type
55 : use m_pawrad, only : pawrad_type
56 : use m_pawtab, only : pawtab_type
57 : use m_pawfgr, only : pawfgr_type
58 : use m_pawrhoij, only : pawrhoij_type
59 : use m_pawcprj, only : pawcprj_type, pawcprj_free
60 : use m_pstat, only : pstat_proc
61 : use m_occ, only : occ_be, occ_fd
62 : use m_lgroup, only : lgroup_t
63 : use m_dvdb, only : dvdb_t
64 : use m_crystal, only : crystal_t
65 : use m_hamiltonian, only : gs_hamiltonian_type, rf_hamiltonian_type
66 : use m_dtset, only : dataset_type
67 : use m_dtfil, only : datafiles_type
68 : use m_wfd, only : wfd_t
69 : use m_gstore, only : gstore_t, gqk_t, GSTORE_GTYPE_KS, GSTORE_GTYPE_GWPT, gstore_read_gtype
70 :
71 : implicit none
72 :
73 : private
74 : public :: gstore_sigmaph
75 :
76 : real(dp),private,parameter :: TOL_EDIFF = 0.001_dp * eV_Ha
77 : !!***
78 :
79 : !----------------------------------------------------------------------
80 :
81 : !!****t* m_gstore_sigmaph/sep_t
82 : !! NAME
83 : !! sep_t
84 : !!
85 : !! FUNCTION
86 : !! Container for the (diagonal) matrix elements of the electron-phonon self-energy
87 : !! in the KS representation i.e. Sigma_eph(omega, T, band, k, spin).
88 : !! Provides methods to compute QP corrections, spectral functions, QP linewidths and
89 : !! save the results to netcdf file.
90 : !!
91 : !! TODO
92 : !! Fix problem with spin parallelism and output of results.
93 : !!
94 : !! SOURCE
95 :
96 : type,public :: sep_t
97 :
98 : integer :: nwr = 0
99 : ! Number of frequency points along the real axis for Sigma(w) and spectral function A(w)
100 : ! Odd number so that the mesh is centered on the KS energy.
101 : ! The spectral function is computed only if nwr > 0 (taken from dtset%nfreqsp)
102 :
103 : integer :: ntemp = 0
104 : ! Number of temperatures.
105 :
106 : logical :: imag_only
107 :
108 : real(dp) :: wr_step
109 : ! Step of the linear mesh along the real axis (Ha units).
110 :
111 : complex(dp) :: ieta = zero
112 : ! Used to shift the poles in the complex plane (Ha units)
113 : ! Corresponds to `i eta` term in equations.
114 :
115 : real(dp),allocatable :: kTmesh(:)
116 : ! kTmesh(ntemp)
117 : ! List of temperatures (kT units).
118 :
119 : real(dp),allocatable :: mu_e(:)
120 : ! mu_e(ntemp)
121 : ! chemical potential of electrons for the different temperatures.
122 :
123 : complex(dp),allocatable :: vals_e0ks(:,:,:)
124 : ! Sigma_eph(omega=eKS, kT, band) for given (ikcalc, spin).
125 : ! Fan-Migdal + Debye-Waller
126 :
127 : complex(dp),allocatable :: fan_vals(:,:,:)
128 : ! (ntemp, nb_k, glob_nk)
129 : ! Fan-Migdal
130 :
131 : complex(dp),allocatable :: fan_stern_vals(:,:,:)
132 : ! (ntemp, nb_k, glob_nk)
133 : ! Fan-Migdal adiabatic Sternheimer part
134 :
135 : complex(dp),allocatable :: dvals_de0ks(:,:,:)
136 : ! (ntemp, nb_k, glob_nk)
137 : ! d Re Sigma_eph(omega, kT, band, kcalc) / d omega (omega=eKS)
138 :
139 : real(dp),allocatable :: dw_vals(:,:,:)
140 : ! dw_vals(ntemp, nb_k, glob_nk) for given (ikcalc, spin)
141 : ! Debye-Waller term (static).
142 :
143 : real(dp),allocatable :: dw_stern_vals(:,:,:)
144 : ! dw_stern_vals(ntemp, nb_k, glob_nk)
145 : ! Debye-Waller Sternheimer term (static) .
146 :
147 : complex(dp),allocatable :: vals_wr(:,:,:,:)
148 : ! vals_wr(nwr, ntemp, nb_k, glob_nk)
149 : ! Sigma_eph(omega, kT, band)
150 : ! enk_KS corresponds to nwr/2 + 1.
151 :
152 : integer :: phmesh_size
153 : ! Number of phonon frequencies in phonon mesh used for Eliashberg functions and
154 : ! and other omega-resolved quantities.
155 :
156 : real(dp),allocatable :: phmesh(:)
157 : ! phmesh(phmesh_size)
158 : ! phonon mesh in Ha.
159 :
160 : real(dp),allocatable :: gfw_vals(:,:,:,:)
161 : ! gfw_vals(phmesh_size, 3, nb_k, nkcalc)
162 : ! Generalized Eliashberg function a2F_{n,k,spin}(w)
163 : ! 1: |g(k,q)|^2 with delta(e_\nk - e_{m\kq})
164 : ! 2: Fan-Migdal in the adiabatic approximation
165 : ! 3: DW contribution in the adiabatic approximation.
166 : ! This array depends on (ikcalc, spin)
167 :
168 : integer :: a2f_ne = 0
169 : ! Number of points in a2f_emesh
170 :
171 : real(dp),allocatable :: a2f_emesh(:)
172 : ! a2f_emesh(a2f_ne)
173 : ! Energy mesh for electrons
174 :
175 : real(dp),allocatable :: a2few(:,:,:,:)
176 : ! a2few(a2f_ne, phmesh_size, nb_k, glob_nk)
177 : ! FM Eliashberg function a2f_\nk(e, w) = \sum_{mq} |g(k,q)|^2 delta(e - e_{m\kq}) delta(w - w_\qnu}
178 : ! This array depends on (ikcalc, spin) and is computed only if prteliash == 3
179 :
180 : real(dp),allocatable :: wrmesh_b(:,:,:)
181 : ! wrmesh_b(nwr, nb_k, glob_nk)
182 : ! Frequency mesh along the real axis (Ha units) used for the different bands
183 : ! Each mesh is **centered** on the corresponding KS energy.
184 :
185 : contains
186 :
187 : procedure :: gather_and_write_results => sep_gather_and_write_results
188 : ! Write main dimensions and header of sigmaph on a netcdf file.
189 :
190 : procedure :: free => sep_free
191 : ! Free dynamic memory
192 : end type sep_t
193 : !!***
194 :
195 : contains
196 : !!***
197 :
198 : !----------------------------------------------------------------------
199 :
200 : !!****f* m_gstore_sigmaph/gstore_sigmaph
201 : !! NAME
202 : !! gstore_sigmaph
203 : !!
204 : !! FUNCTION
205 : !! Compute diagonal matrix elements of the e-ph self-energy (Fan Migdal + Debye Waller).
206 : !! using precomputed e-ph matrix elements.
207 : !! See also m_sigmaph, for a version in which the g-matrix elements are computed on-the-fly.
208 : !!
209 : !! INPUTS
210 : !! wfk0_path=String with the path to the GS unperturbed WFK file.
211 : !! ngfft(18),ngfftf(18)=Coarse and Fine FFT meshes.
212 : !! dtset<dataset_type>=All input variables for this dataset.
213 : !! dtfil<datafiles_type>=Variables related to files.
214 : !! cryst: Crystalline structure
215 : !! ks_ebands<ebands_t>=The KS band structure (energies, occupancies, k-weights...)
216 : !! qp_ebands<ebands_t>=The QP band structure (energies, occupancies, k-weights...)
217 : !! dvdb<dbdb_type>=Database with the DFPT SCF potentials.
218 : !! ifc<ifc_type>=interatomic force constants and corresponding real space grid info.
219 : !! pawfgr <type(pawfgr_type)>=fine grid parameters and related data
220 : !! pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data.
221 : !! psps<pseudopotential_type>=Variables related to pseudopotentials.
222 : !!
223 : !! TODO:
224 : !! For off-diagonal terms, see https://journals.aps.org/prb/abstract/10.1103/PhysRevB.101.121102
225 : !!
226 : !! OUTPUT
227 : !! Write results to ab_out and netcdf files.
228 : !!
229 : !! SOURCE
230 :
231 7 : subroutine gstore_sigmaph(wfk0_path, ngfft, ngfftf, dtset, dtfil, cryst, ks_ebands, qp_ebands, wfk0_hdr, dvdb, ifc, &
232 7 : pawfgr, pawtab, psps, mpi_enreg, comm)
233 :
234 : !Arguments ------------------------------------
235 : !scalars
236 : character(len=*),intent(in) :: wfk0_path
237 : type(dataset_type),intent(in) :: dtset
238 : type(datafiles_type),intent(in) :: dtfil
239 : type(crystal_t),intent(in) :: cryst
240 : type(ebands_t),target,intent(in) :: ks_ebands, qp_ebands
241 : type(hdr_type),intent(in) :: wfk0_hdr
242 : type(dvdb_t),intent(inout) :: dvdb
243 : type(ifc_type),target,intent(in) :: ifc
244 : type(pseudopotential_type),intent(in) :: psps
245 : type(pawfgr_type),intent(in) :: pawfgr
246 : type(mpi_type),intent(inout) :: mpi_enreg
247 : integer,intent(in) :: comm
248 : !arrays
249 : integer,intent(in) :: ngfft(18),ngfftf(18)
250 : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
251 :
252 : !Local variables-------------------------------
253 : integer,parameter :: master = 0, cplex1 = 1, pawread0 = 0, ndat1 = 1, istwfk_1 = 1
254 : integer,parameter :: LOG_MODQ = 100, LOG_MODK = 1, g2mode_AA = 1, g2mode_KS_GWPT = 2
255 : integer :: n1, n2, n3, n4, n5, n6, nb_k, nb_kq, glob_nk, ntemp, cplex, my_npert, use_lgk, iw
256 : integer :: spin, my_is, my_ik, my_iq, my_ip, in_k, im_kq, ierr, gap_err, my_rank, ip1, ip2, nu, ipc, idir, ipert
257 : integer :: it, ik_ibz, ikq_ibz, band_k, band_kq, timrev_k, ii, ikcalc, natom, natom3, nsppol, nspden, nspinor, nkpt !,ik_bz
258 : integer :: isym_k,isym_kq,trev_k,trev_kq, with_cplex
259 : integer :: istwf_k, istwf_kq, npw_k, npw_kq, nkpg_kq, nfft, nfftf, mgfft, mgfftf, nkpg
260 : integer :: usecprj, mpw, ibsum_kq, band_me, u1_band, ncid, ncerr
261 : real(dp) :: wqnu, gkq2, weight_q, eig0nk, eig0mk, eig0mkq, ediff, gmod2, hmod2, gdw2, rfact, gdw2_stern !, rtmp !,nqnu,gkq2,gkq2_pf,
262 : real(dp) :: cpu_kk, wall_kk, gflops_kk, cpu_qq, wall_qq, gflops_qq, cpu_all, wall_all, gflops_all
263 : real(dp) :: estep
264 : logical :: q_is_gamma, intra_band, same_band, isirr_k, isirr_kq, stern_use_cache, print_time_kk, print_time_qq, gstore_from_file
265 : complex(dp) :: cfact !, sig_cplx
266 : character(len=5000) :: msg, qq_bz_string !, kk_string
267 : character(len=fnlen) :: path
268 : character(len=abi_slen) :: gtype
269 7 : type(gaps_t) :: gaps
270 7 : type(lgroup_t) :: lg_myk
271 28 : type(gstore_t) :: gstore
272 7 : type(sep_t) :: sigma
273 7 : type(hdr_type) :: pot_hdr
274 364 : type(crystal_t) :: pot_cryst
275 7 : type(wfd_t) :: wfd
276 : !type(u1_cache_t) :: u1c
277 : type(ebands_t),pointer :: ebands
278 7 : type(stern_t) :: stern
279 7 : type(gs_hamiltonian_type) :: gs_ham_kq
280 7 : type(rf_hamiltonian_type) :: rf_ham_kq
281 : !arrays
282 : integer :: gmax(3), g0_k(3), g0_kq(3), work_ngfft(18), units(2), my_kqmap(6)
283 21 : integer,allocatable :: phmodes_skip(:), gbound_kq(:,:), nband(:,:), wfd_istwfk(:), kg_kq(:,:) !, kg_k(:,:), gaussw_qnu(:)
284 14 : real(dp) :: kk(3), kk_ibz(3), kq_ibz(3), qpt(3), kq(3), fermie1_idir_ipert(3,cryst%natom), dotri(2)
285 14 : real(dp),allocatable :: vtrial(:,:), work(:,:,:,:), kinpw_k(:), kinpw_kq(:),kpg_kq(:,:),kpg_k(:,:)
286 21 : real(dp),allocatable :: ffnl_k(:,:,:,:),ffnl_kq(:,:,:,:),ph3d_k(:,:,:),ph3d_kq(:,:,:),v1scf(:,:,:,:)
287 7 : real(dp),allocatable :: dtw_weights(:,:),dwargs(:) !,alpha_mrta(:)
288 7 : real(dp),allocatable :: delta_e_minus_emkq(:) ! gkq2_lr(:,:,:)
289 14 : real(dp) :: displ_red_nu(2, 3, cryst%natom)
290 14 : real(dp),allocatable :: cg1s_kq(:,:,:,:), h1kets_kq_allperts(:,:,:,:)
291 14 : real(dp) :: vec_natom3(2, 3*cryst%natom) ! zpr_frohl_sphcorr(3*cryst%natom),
292 14 : real(dp),allocatable :: bra_kq(:,:), kets_k(:,:,:)
293 7 : real(dp),allocatable :: stern_ppb(:,:,:,:), stern_fan_t(:), stern_dw(:,:,:,:)
294 21 : real(dp),allocatable :: ph1d(:,:),vlocal(:,:,:,:),vlocal1(:,:,:,:,:), rfact_t(:), nqnu_t(:), f_mkq(:)
295 14 : complex(dp),allocatable :: cfact_t(:), cfact2_t(:), cfact_wr(:), tpp_red(:,:) !,fmw_frohl_sphcorr(:,:,:,:),
296 7 : logical,allocatable :: bks_mask(:,:,:),keep_ur(:,:,:)
297 7 : type(pawrhoij_type),allocatable :: pot_pawrhoij(:)
298 7 : type(pawcprj_type),allocatable :: cwaveprj0(:,:), cwaveprj(:,:)
299 : !----------------------------------------------------------------------
300 :
301 : !ebands => ks_ebands
302 7 : ebands => qp_ebands
303 :
304 21 : my_rank = xmpi_comm_rank(comm); units = [std_out, ab_out]
305 7 : natom = cryst%natom; natom3 = 3 * cryst%natom; nkpt = ebands%nkpt
306 7 : nsppol = dtset%nsppol; nspden = dtset%nspden; nspinor = dtset%nspinor
307 :
308 7 : call cwtime(cpu_all, wall_all, gflops_all, "start")
309 7 : call wrtout(units, " Computing Fan-Migdal + DW self-energy from GSTORE.nc", pre_newlines=1)
310 :
311 : ! gtype is only meaningful when reading a pre-existent GSTORE.nc file (gstore_read_gtype
312 : ! opens dtfil%filgstorein). When building gstore on the fly from ABIWAN.nc + GWAN.nc
313 : ! (see gstore%init_or_from_ncpath below), there is no file to read gtype from, and the
314 : ! interpolated construction always produces a single "atom"-representation g, never "gwpt".
315 7 : gtype = GSTORE_GTYPE_KS
316 7 : if (dtfil%filgstorein /= ABI_NOFILE) call gstore_read_gtype(dtfil%filgstorein, gtype, comm)
317 :
318 7 : with_cplex = 1
319 7 : if (gtype == GSTORE_GTYPE_GWPT) then
320 : ! Decide if self-energies should be computed with |g|^2 or g^KS g^GWPT.
321 11 : select case (dtset%gwpt_g2mode)
322 : case (g2mode_AA)
323 : with_cplex = 1
324 5 : call wrtout(units, " Using e-ph self-energy expression with |g|^2")
325 : case (g2mode_KS_GWPT)
326 1 : with_cplex = 2
327 1 : call wrtout(units, " Using e-ph self-energy expression with g^*_KS g_GWPT")
328 : case default
329 6 : ABI_ERROR(sjoin("Invalid dtset%gwpt_g2mode:", itoa(dtset%gwpt_g2mode)))
330 : end select
331 : else
332 : with_cplex = 1
333 1 : call wrtout(units, " Using e-ph self-energy expression with |g|^2")
334 : end if
335 :
336 : ! Init gstore and MPI grid from file and dtset (or, alternatively, via Wannier
337 : ! interpolation from ABIWAN.nc + GWAN.nc, see gstore_init_or_from_ncpath).
338 : ! The Fan-Migdal SE requires |g(k,q)|^2 as well as g2DW in the phonon representation.
339 : call gstore%init_or_from_ncpath(with_cplex, dtset, dtfil, wfk0_hdr, cryst, ebands, ifc, &
340 7 : "phonon", dtset%gstore_gname, .True., comm, gstore_from_file)
341 7 : if (gstore_from_file) then
342 7 : call wrtout(units, " Gstore built by reading a pre-existent GSTORE.nc file")
343 : else
344 0 : call wrtout(units, " Gstore built on the fly via Wannier interpolation (ABIWAN.nc + GWAN.nc)")
345 : end if
346 :
347 : ! Consistency check.
348 7 : ierr = 0
349 7 : if (gstore%qzone /= "bz") then
350 0 : ABI_ERROR_NOSTOP("gstore_sigmaph assumes qzone == `bz`", ierr)
351 : end if
352 7 : if (gstore%has_used_lgq /= 0) then
353 0 : ABI_ERROR_NOSTOP("gstore_sigmaph does not support use_lgq /=0.", ierr)
354 : end if
355 7 : if (ierr /= 0) then
356 : write(msg,'(a,i0,5a)')&
357 0 : 'Checking consistency of input data against itself gave ',ierr,' inconsistencies.',ch10,&
358 0 : 'The details of the problems can be FOUND ABOVE (or in output or log file), in an earlier WARNING.',ch10,&
359 0 : 'In parallel, the details might not even be printed there. Then, try running in sequential to see the details.'
360 0 : ABI_ERROR(msg)
361 : end if
362 :
363 : ! Check consistency of little group options
364 7 : ABI_CHECK(gstore%check_little_group(dtset, msg) == 0, msg)
365 :
366 7 : use_lgk = dtset%gstore_use_lgk
367 7 : if (gstore%has_used_lgk /= 0) use_lgk = gstore%has_used_lgk
368 7 : if (use_lgk == 0) then
369 0 : call wrtout(units, " Little group operations of the k-point won't be used to symmetry reduce the integral in q-space.")
370 : end if
371 : if (use_lgk /= 0) then
372 7 : call wrtout(units, " Little group operations of the k-point will be used to symmetry reduce the integral in q-space.")
373 : end if
374 :
375 : ! FFT meshes from input file, not necessarily equal to the ones found in the external files.
376 28 : nfftf = product(ngfftf(1:3)); mgfftf = maxval(ngfftf(1:3))
377 49 : nfft = product(ngfft(1:3)) ; mgfft = maxval(ngfft(1:3))
378 7 : n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3)
379 7 : n4 = ngfft(4); n5 = ngfft(5); n6 = ngfft(6)
380 :
381 : ! Initialize parameters in sigma object.
382 7 : sigma%imag_only = .False.
383 7 : sigma%ieta = + j_dpc * dtset%zcut
384 :
385 : ! Build (linear) mesh of K * temperatures. tsmesh(1:3) = [start, step, num]
386 7 : call dtset%get_ktmesh(sigma%ntemp, sigma%kTmesh)
387 :
388 : ! Compute the chemical potential at the different physical temperatures with Fermi-Dirac.
389 21 : ABI_MALLOC(sigma%mu_e, (sigma%ntemp))
390 20 : sigma%mu_e(:) = ebands%fermie
391 7 : if (dtset%eph_fermie == zero) then
392 7 : call ebands%get_muT_with_fd(sigma%ntemp, sigma%ktmesh, dtset%spinmagntarget, dtset%prtvol, sigma%mu_e, gstore%comm)
393 : end if
394 :
395 : ! Compute gaps.
396 7 : gaps = ebands%get_gaps(gap_err)
397 7 : if (gap_err /= 0) then
398 0 : ABI_ERROR("Cannot compute fundamental and direct gap (likely metal)")
399 : end if
400 :
401 7 : if (my_rank == master) then
402 : call gaps%print(units, kTmesh=sigma%ktmesh, mu_e=sigma%mu_e, &
403 7 : header="Gaps, band edges and relative position wrt Fermi level")
404 : end if
405 7 : call gaps%free()
406 :
407 : ! Frequency mesh for sigma(w) and spectral functions.
408 7 : call dtset%get_wrmesh_for_sigeph(sigma%nwr, sigma%wr_step)
409 :
410 : ! Setup a mask to skip accumulating the contribution of certain phonon modes.
411 7 : call ephtk_set_phmodes_skip(dtset%natom, dtset%eph_phrange, phmodes_skip)
412 :
413 : !==== Initialize most of the Hamiltonian (and derivative) ====
414 : ! 1) Allocate all arrays and initialize quantities that do not depend on k and spin.
415 : ! 2) Perform the setup needed for the non-local factors:
416 : !
417 : ! Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk.
418 : ! PAW: Initialize the overlap coefficients and allocate the Dij coefficients.
419 :
420 : ! Get one-dimensional structure factor information on the coarse grid.
421 21 : ABI_MALLOC(ph1d, (2,3*(2*mgfft+1)*natom))
422 7 : call getph(cryst%atindx, natom, n1, n2, n3, ph1d, cryst%xred)
423 :
424 7 : usecprj = 0
425 : call gs_ham_kq%init(psps, pawtab, nspinor, nsppol, nspden, natom,&
426 : dtset%typat, cryst%xred, nfft, mgfft, ngfft, cryst%rprimd, dtset%nloalg,&
427 : comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab, mpi_spintab=mpi_enreg%my_isppoltab,&
428 7 : usecprj=usecprj, ph1d=ph1d, nucdipmom=dtset%nucdipmom, gpu_option=dtset%gpu_option)
429 :
430 7 : if (dtset%eph_stern /= 0 .and. .not. sigma%imag_only) then
431 : ! Prepare call to Sternheimer solver.
432 : ! The static correction to FM_nk is:
433 : ! \sum_{qnu} (2n_qnu + 1) <H^1_{qnu} psi_nk| psi^1_{nk; qnu}>
434 :
435 : ! Allocate work space arrays.
436 : ! vtrial and vlocal are required for Sternheimer (H0). DFPT routines do not need it.
437 : ! Note nvloc in vlocal (we will select one/four spin components afterwards)
438 173635 : ABI_CALLOC(vtrial, (nfftf, nspden))
439 191793 : ABI_CALLOC(vlocal, (n4, n5, n6, gs_ham_kq%nvloc))
440 :
441 : ! Read the GS potential (vtrial) from input POT file.
442 : ! In principle one may store vtrial in the DVDB but getpot_filepath is simpler to implement.
443 7 : call wrtout(units, sjoin(" Reading GS KS potential for Sternheimer from: ", dtfil%filpotin))
444 : call read_rhor(dtfil%filpotin, cplex1, dtset%nspden, nfftf, ngfftf, pawread0, mpi_enreg, vtrial, pot_hdr, pot_pawrhoij, comm, &
445 7 : allow_interp=.True., want_varname="vtrial")
446 7 : pot_cryst = pot_hdr%get_crystal()
447 7 : if (gstore%cryst%compare(pot_cryst, header=" Comparing input crystal with POT crystal") /= 0) then
448 0 : ABI_ERROR("Crystal structure from WFK and POT do not agree! Check messages above!")
449 : end if
450 7 : call pot_cryst%free(); call pot_hdr%free()
451 :
452 : ! Initialize the wave function descriptor.
453 : ! Only wavefunctions for the symmetrical image of the k/k+q wavevectors treated by this MPI rank are stored.
454 28 : ABI_MALLOC(nband, (nkpt, nsppol))
455 35 : ABI_MALLOC(bks_mask, (dtset%mband, nkpt, nsppol))
456 28 : ABI_MALLOC(keep_ur, (dtset%mband, nkpt ,nsppol))
457 :
458 1385 : nband = dtset%mband; bks_mask = .False.; keep_ur = .False.
459 :
460 : ! Initialize bks_mask
461 7 : call gstore%fill_bks_mask(dtset%mband, nkpt, nsppol, bks_mask)
462 :
463 : ! mpw is the maximum number of plane-waves over k and k+q where k and k+q are in the BZ.
464 : ! we also need the max components of the G-spheres (k, k+q) in order to allocate the workspace array work
465 : ! that will be used to symmetrize the wavefunctions in G-space.
466 7 : call gstore%get_mpw_gmax(dtset%ecut, mpw, gmax)
467 :
468 : ! Init work_ngfft
469 28 : gmax = gmax + 4 ! FIXME: this is to account for umklapp, should also consider Gamma-only and istwfk
470 28 : gmax = 2*gmax + 1
471 :
472 7 : call ngfft_seq(work_ngfft, gmax)
473 : !write(std_out,*)"work_ngfft(1:3): ",work_ngfft(1:3)
474 35 : ABI_MALLOC(work, (2, work_ngfft(4), work_ngfft(5), work_ngfft(6)))
475 :
476 : ! Allocate PW-arrays. Note mpw in kg_kq
477 21 : ABI_MALLOC(kg_kq, (3, mpw))
478 :
479 : ! Impose istwfk=1 for all k points. This is also done in respfn (see inkpts)
480 : ! wfd_read_wfk will handle a possible conversion if WFK contains istwfk /= 1.
481 21 : ABI_MALLOC(wfd_istwfk, (nkpt))
482 48 : wfd_istwfk = 1
483 :
484 : ! TODO: Possible problem if nband < gstore%nband
485 : call wfd%init(cryst, pawtab, psps, keep_ur, dtset%mband, nband, nkpt, nsppol, bks_mask,&
486 : dtset%nspden, nspinor, dtset%ecut, dtset%ecutsm, dtset%dilatmx, wfd_istwfk, ebands%kptns, ngfft,&
487 7 : dtset%nloalg, dtset%prtvol, dtset%pawprtvol, comm)
488 :
489 14 : call wfd%print([std_out], header="Wavefunctions for Sternheimer.")
490 7 : call pstat_proc%print(_PSTAT_ARGS_)
491 :
492 7 : ABI_FREE(nband)
493 7 : ABI_FREE(bks_mask)
494 7 : ABI_FREE(keep_ur)
495 7 : ABI_FREE(wfd_istwfk)
496 :
497 : ! Read wavefunctions.
498 7 : call wfd%read_wfk(wfk0_path, iomode_from_fname(wfk0_path))
499 :
500 : ! Open the DVDB file
501 7 : call dvdb%open_read(ngfftf, xmpi_comm_self)
502 7 : ABI_CHECK(dvdb%has_fields("pot1", msg), msg)
503 :
504 : ! Activate parallelism over perturbations
505 7 : call gstore%set_perts_distrib(cryst, dvdb, my_npert)
506 :
507 : ! Prepare FT interpolation of DFPT potentials.
508 7 : call dvdb%ftinterp_setup(dtset%ddb_ngqpt, gstore%qptopt, 1, dtset%ddb_shiftq, nfftf, ngfftf, xmpi_comm_self)
509 :
510 28 : ABI_MALLOC(tpp_red, (natom3, natom3))
511 28 : ABI_MALLOC(gbound_kq, (2*wfd%mgfft+8, 2))
512 : end if ! eph_stern /= 0
513 :
514 : ! Compute phonon frequency mesh.
515 7 : call ifc%get_phmesh(dtset%ph_wstep, sigma%phmesh_size, sigma%phmesh)
516 :
517 7 : sigma%a2f_ne = 0
518 7 : if (dtset%prteliash == 3) then
519 : ! TODO: dosdeltae should have a default value.
520 : ! TODO: Use logmesh/double mesh for electrons?
521 0 : estep = dtset%dosdeltae; if (estep <= zero) estep = 0.05 * eV_Ha
522 0 : sigma%a2f_ne = nint((maxval(ebands%eig) - minval(ebands%eig)) / estep) + 1
523 0 : if (my_rank == master) then
524 0 : write(std_out, *)" Computing a2f with ", sigma%a2f_ne, " points for electrons and ", sigma%phmesh_size, " points for phonons."
525 0 : write(std_out, *)" doseltae:", estep, ", tsmear:", dtset%tsmear
526 : end if
527 0 : ABI_MALLOC(sigma%a2f_emesh, (sigma%a2f_ne))
528 0 : sigma%a2f_emesh = arth(minval(ebands%eig), estep, sigma%a2f_ne)
529 : end if
530 :
531 : ! Allocate workspace arrays for Eliashberg functions
532 7 : if (dtset%prteliash /= 0) then
533 0 : ABI_MALLOC(dtw_weights, (sigma%phmesh_size, 2))
534 0 : ABI_MALLOC(dwargs, (sigma%phmesh_size))
535 0 : if (sigma%a2f_ne > 0) then
536 0 : ABI_MALLOC(delta_e_minus_emkq, (sigma%a2f_ne))
537 : end if
538 : end if
539 :
540 : ! Allocate work space arrays used inside the loops. Then we are ready to go!
541 7 : ntemp = sigma%ntemp
542 21 : ABI_MALLOC(nqnu_t, (ntemp))
543 14 : ABI_MALLOC(f_mkq, (ntemp))
544 21 : ABI_MALLOC(cfact_t, (ntemp))
545 14 : ABI_MALLOC(cfact2_t, (ntemp))
546 14 : ABI_MALLOC(rfact_t, (ntemp))
547 14 : ABI_MALLOC(stern_fan_t, (ntemp))
548 :
549 7 : call pstat_proc%print(_PSTAT_ARGS_)
550 :
551 : ! Create netcdf file (only master works, HDF5 + MPI-IO is handled afterwards by reopening the file inside ncwrite_comm)
552 7 : path = strcat(dtfil%filnam_ds(4), "_GSEPH.nc")
553 7 : if (my_rank == master) then
554 : ! Master creates the netcdf file used to store the results of the calculation.
555 7 : NCF_CHECK(nctk_open_create(ncid, path, xmpi_comm_self))
556 : !NCF_CHECK(wfk_hdr%ncwrite(ncid, fform_from_ext("GSEPH.nc"), nc_define=.True.))
557 7 : NCF_CHECK(cryst%ncwrite(ncid))
558 7 : NCF_CHECK(ebands%ncwrite(ncid))
559 :
560 : ! Add dimensions.
561 : ncerr = nctk_def_dims(ncid, [ &
562 : nctkdim_t("nsppol", nsppol), nctkdim_t("ntemp", ntemp), nctkdim_t("natom3", 3 * natom3), &
563 : !nctkdim_t("glob_nk", sigma%glob_nk), nctkdim_t("max_nbcalc", sigma%max_nbcalc), &
564 : nctkdim_t("phmesh_size", sigma%phmesh_size) &
565 : !nctkdim_t("nqibz", sigma%nqibz), nctkdim_t("nqbz", sigma%nqbz)
566 : ], &
567 35 : defmode=.True.)
568 7 : NCF_CHECK(ncerr)
569 :
570 7 : if (sigma%nwr > 0) then
571 2 : NCF_CHECK(nctk_def_dims(ncid, [nctkdim_t("nwr", sigma%nwr)]))
572 : end if
573 7 : if (dtset%prteliash == 3) then
574 0 : NCF_CHECK(nctk_def_dims(ncid, [nctkdim_t("a2f_ne", sigma%a2f_ne)]))
575 : end if
576 :
577 : !ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: &
578 : ! "symsigma", "nbsum", "bsum_start", "bsum_stop", "symdynmat", &
579 : ! "ph_intmeth", "eph_intmeth", "qint_method", "eph_transport", &
580 : ! "imag_only", "symv1scf", "dvdb_add_lr", "mrta", "ibte_prep", "eph_prtscratew", "eph_ahc_type"])
581 : !NCF_CHECK(ncerr)
582 : !ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: &
583 : ! "eta", "wr_step", "eph_fsewin", "eph_fsmear", "eph_extrael", "eph_fermie", &
584 : ! "ph_wstep", "ph_smear", "eph_phwinfact"])
585 : !NCF_CHECK(ncerr)
586 :
587 : ! Define arrays with results.
588 : ncerr = nctk_def_arrays(ncid, [ &
589 : !nctkarr_t("ngqpt", "int", "three"), &
590 : !nctkarr_t("eph_ngqpt_fine", "int", "three"), &
591 : !nctkarr_t("eph_phrange", "int", "two"), &
592 : !nctkarr_t("eph_phrange_w", "dp", "two"), &
593 : !nctkarr_t("ddb_ngqpt", "int", "three"), &
594 : !nctkarr_t("ph_ngqpt", "int", "three"), &
595 : !!nctkarr_t("frohl_params", "dp", "four"), &
596 : !nctkarr_t("bstart_ks", "int", "glob_nk, nsppol"), &
597 : !nctkarr_t("nbcalc_ks", "int", "glob_nk, nsppol"), &
598 : !nctkarr_t("kcalc", "dp", "three, glob_nk"), &
599 : !nctkarr_t("kcalc2ibz", "int", "glob_nk, six"), &
600 : nctkarr_t("kTmesh", "dp", "ntemp"), &
601 : nctkarr_t("mu_e", "dp", "ntemp") &
602 21 : ])
603 7 : NCF_CHECK(ncerr)
604 :
605 : ! ======================================================
606 : ! Write data that do not depend on the (kpt, spin) loop.
607 : ! ======================================================
608 7 : NCF_CHECK(nctk_set_datamode(ncid))
609 :
610 7 : ii = 0; if (sigma%imag_only) ii = 1
611 : !ncerr = nctk_write_iscalars(ncid, [character(len=nctk_slen) :: &
612 : ! "symsigma", "nbsum", "bsum_start", "bsum_stop", &
613 : ! "symdynmat", "ph_intmeth", "eph_intmeth", "qint_method", &
614 : ! "eph_transport", "imag_only", "symv1scf", "dvdb_add_lr", "mrta", "ibte_prep", "eph_prtscratew", "eph_ahc_type"], &
615 : ! [self%symsigma, self%nbsum, self%bsum_start, self%bsum_stop, &
616 : ! dtset%symdynmat, dtset%ph_intmeth, dtset%eph_intmeth, self%qint_method, dtset%eph_transport, ii, &
617 : ! dtset%symv1scf, dtset%dvdb_add_lr, self%mrta, dtset%ibte_prep, dtset%eph_prtscratew, dtset%eph_ahc_type])
618 : !NCF_CHECK(ncerr)
619 : !ncerr = nctk_write_dpscalars(ncid, [character(len=nctk_slen) :: &
620 : ! "eta", "wr_step", "eph_fsewin", "eph_fsmear", "eph_extrael", "eph_fermie", "ph_wstep", "ph_smear", "eph_phwinfact"], &
621 : ! [aimag(self%ieta), self%wr_step, dtset%eph_fsewin, dtset%eph_fsmear, dtset%eph_extrael, dtset%eph_fermie, &
622 : ! dtset%ph_wstep, dtset%ph_smear, dtset%eph_phwinfact])
623 : !NCF_CHECK(ncerr)
624 :
625 : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ngqpt"), sigma%ngqpt))
626 : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "eph_ngqpt_fine"), dtset%eph_ngqpt_fine))
627 : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ddb_ngqpt"), dtset%ddb_ngqpt))
628 : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ph_ngqpt"), dtset%ph_ngqpt))
629 : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "sigma_ngkpt"), dtset%sigma_ngkpt))
630 : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "sigma_erange"), dtset%sigma_erange))
631 : !!NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "frohl_params"), dtset%frohl_params))
632 : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "eph_phrange"), dtset%eph_phrange))
633 : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "eph_phrange_w"), dtset%eph_phrange_w))
634 : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "bstart_ks"), sigma%bstart_ks))
635 : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "nbcalc_ks"), sigma%nbcalc_ks))
636 : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "kcalc"), sigma%kcalc))
637 : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "kcalc2ibz"), sigma%kcalc2ibz))
638 7 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "kTmesh"), sigma%kTmesh))
639 7 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "mu_e"), sigma%mu_e))
640 : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "phmesh"), sigma%phmesh))
641 : !if (dtset%prteliash == 3) then
642 : ! NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "a2f_emesh"), sigma%a2f_emesh))
643 : !end if
644 : !NCF_CHECK(nf90_close(ncid))
645 : end if ! master
646 :
647 7 : call xmpi_barrier(comm)
648 :
649 : ! Now reopen the file inside ncwrite_comm to perform parallel-IO (required for k-point parallelism).
650 : !if (self%ncwrite_comm%value /= xmpi_comm_null) then
651 : ! NCF_CHECK(nctk_open_modify(ncid, path, self%ncwrite_comm%value))
652 : ! NCF_CHECK(nctk_set_datamode(ncid))
653 : !end if
654 :
655 7 : call wrtout(std_out, " Begin computation of the self-energy matrix elements.")
656 :
657 7 : if (dtset%userib /= 0) then
658 0 : call wrtout(units, sjoin(" userib /= 0 => Include only one q-point in the integration. qpt", ktoa(dtset%qptn)))
659 : ! Notes:
660 : ! 1) Each q-point contribution is weighted by (multiplicity / nqbz), not by 1 / nqbz.
661 : ! 2) When use_lgk is enabled, there is no guarantee that the q-point
662 : ! specified in dtset%qptn belongs to the zone (IBZ_k). One shouls check the log file for messages
663 : end if
664 :
665 : ! Loop over collinear spins.
666 14 : do my_is=1,gstore%my_nspins
667 7 : associate (gqk => gstore%gqk(my_is), cryst => gstore%cryst)
668 7 : spin = gstore%my_spins(my_is); nb_k = gqk%nb_k; nb_kq = gqk%nb_kq; glob_nk = gqk%glob_nk
669 :
670 : ! Allocate arrays for self-energy matrix elements.
671 207 : ABI_CALLOC(sigma%vals_e0ks, (ntemp, nb_k, glob_nk))
672 200 : ABI_CALLOC(sigma%dvals_de0ks, (ntemp, nb_k, glob_nk))
673 200 : ABI_CALLOC(sigma%fan_vals, (ntemp, nb_k, glob_nk))
674 200 : ABI_CALLOC(sigma%fan_stern_vals, (ntemp, nb_k, glob_nk))
675 207 : ABI_CALLOC(sigma%dw_vals, (ntemp, nb_k, glob_nk))
676 200 : ABI_CALLOC(sigma%dw_stern_vals, (ntemp, nb_k, glob_nk))
677 :
678 : ! Prepare computation of Sigma_{nk}(w) and spectral function.
679 7 : if (sigma%nwr > 0) then
680 3243 : ABI_CALLOC(sigma%vals_wr, (sigma%nwr, ntemp, nb_k, glob_nk))
681 814 : ABI_CALLOC(sigma%wrmesh_b, (sigma%nwr, nb_k, glob_nk))
682 3 : ABI_MALLOC(cfact_wr, (sigma%nwr))
683 : end if
684 :
685 : ! Prepare calculation of generalized Eliashberg functions
686 : ! prteliash == 0 deactivates computation (default).
687 7 : if (dtset%prteliash /= 0) then
688 0 : ABI_MALLOC(sigma%gfw_vals, (sigma%phmesh_size, 3, nb_k, glob_nk))
689 : end if
690 :
691 7 : if (dtset%prteliash == 3) then
692 0 : ABI_CALLOC(sigma%a2few, (sigma%a2f_ne, sigma%phmesh_size, nb_k, glob_nk))
693 : end if
694 :
695 4405 : ABI_CALLOC(stern_dw, (2, natom3, natom3, nb_k))
696 4398 : ABI_CALLOC(stern_ppb, (2, natom3, natom3, nb_k))
697 :
698 : ! Loop over my k-points in |n,k>.
699 7 : ABI_CHECK(gqk%my_nq /= 0, "gqm%my_nk cannot be zero here!")
700 :
701 19 : do my_ik=1,gqk%my_nk
702 48 : kk = gqk%my_kpts(:, my_ik)
703 12 : print_time_kk = my_rank == 0 .and. (my_ik <= LOG_MODK .or. mod(my_ik, LOG_MODK) == 0)
704 12 : if (print_time_kk) call cwtime(cpu_kk, wall_kk, gflops_kk, "start")
705 :
706 12 : ik_ibz = gqk%my_k2ibz(1, my_ik); isym_k = gqk%my_k2ibz(2, my_ik)
707 48 : trev_k = gqk%my_k2ibz(6, my_ik); g0_k = gqk%my_k2ibz(3:5, my_ik)
708 48 : isirr_k = (isym_k == 1 .and. trev_k == 0 .and. all(g0_k == 0))
709 0 : ABI_CHECK(isirr_k, "For the time being the k-point in Sigma_{nk} must be in the IBZ")
710 : kk_ibz = ebands%kptns(:,ik_ibz)
711 :
712 : ! Will store results in sigma% using glob_ik index.
713 12 : ikcalc = gqk%my_k2glob(my_ik)
714 :
715 : ! Compute the little group of the k-point so that we can sum g(k,q) only for q in the IBZ_k.
716 12 : if (use_lgk /= 0) then
717 12 : timrev_k = kpts_timrev_from_kptopt(ebands%kptopt)
718 12 : call lg_myk%init(cryst, kk, timrev_k, gstore%nqbz, gstore%qbz, gstore%nqibz, gstore%qibz, xmpi_comm_self)
719 : end if
720 :
721 12 : if (sigma%nwr > 0) then
722 : ! Prepare computation of Sigma_{nk}(w) and spectral function.
723 : ! Build linear mesh **centered** around the KS energy.
724 5 : do in_k=1,nb_k
725 4 : band_k = in_k + gqk%bstart_k - 1
726 4 : eig0nk = ebands%eig(band_k, ik_ibz, spin) - sigma%wr_step * (sigma%nwr / 2)
727 809 : sigma%wrmesh_b(:,in_k,ikcalc) = arth(eig0nk, sigma%wr_step, sigma%nwr)
728 : end do
729 : end if
730 :
731 12 : if (dtset%eph_stern /= 0) then
732 12 : npw_k = wfd%npwarr(ik_ibz); istwf_k = wfd%istwfk(ik_ibz)
733 : call gs_ham_kq%eph_setup_k("k", kk, istwfk_1, npw_k, wfd%kdata(ik_ibz)%kg_k, dtset, cryst, psps, & ! in
734 12 : nkpg, kpg_k, ffnl_k, kinpw_k, ph3d_k, xmpi_comm_self) ! out
735 : end if
736 :
737 : ! Sum over my q-points.
738 12 : ABI_CHECK(gqk%my_nq /= 0, "gqm%my_nq cannot be zero here!")
739 :
740 220 : do my_iq=1,gqk%my_nq
741 1040 : call gqk%myqpt(my_iq, gstore, weight_q, qpt); q_is_gamma = sum(qpt**2) < tol14
742 :
743 : ! weight_q is computed here. It depends whether we are summing over the full BZ or IBZ_k.
744 : ! IMPORTANT: We cannot cycle is my_iq == 1 as this is the iteration in which we broadcast stern_dw if eph_stern /= 0.
745 : ! Also weight_q should be set to zero if q is not in the IBZ_k when my_iq == 1.
746 208 : weight_q = one / gstore%nqbz
747 208 : if (use_lgk /= 0) then
748 208 : ii = lg_myk%findq_ibzk(qpt); if (ii == -1 .and. my_iq /= 1) cycle
749 0 : weight_q = zero
750 69 : if (ii /= -1) weight_q = lg_myk%weights(ii)
751 : end if
752 :
753 : ! Select contribution from a certaing q-point
754 : ! Don't cycle as this will interfere with parallelism over q-points and DW
755 69 : if (dtset%userib /= 0) then
756 0 : if (any(abs(qpt - dtset%qptn) > tol14)) then
757 0 : weight_q = zero
758 : else
759 0 : write(msg, *) "weight for qpt is", weight_q, " with multiplicity:", weight_q * gstore%nqbz
760 0 : call wrtout(std_out, msg)
761 : end if
762 : end if
763 :
764 : !iq_bz = gqk%my_q2bz(my_iq); qq_is_gamma = sum(qq_bz**2) < tol14
765 69 : qq_bz_string = ktoa(qpt)
766 :
767 69 : print_time_qq = my_rank == 0 .and. (my_iq <= LOG_MODQ .or. mod(my_iq, LOG_MODQ) == 0)
768 69 : if (print_time_qq) then
769 69 : call cwtime(cpu_qq, wall_qq, gflops_qq, "start")
770 69 : call inds2str(0, sjoin(" Computing Sigma_eph for qq_bz:", qq_bz_string), my_iq, gqk%my_nq, gqk%glob_nq, msg)
771 69 : call wrtout(std_out, sjoin(msg, ", and spin:", itoa(spin)), pre_newlines=1)
772 : end if
773 :
774 : ! Find the image of k+q in the IBZ.
775 276 : kq = kk + qpt
776 69 : if (kpts_map("symrel", ebands%kptopt, cryst, gstore%krank_ibz, 1, kq, my_kqmap) /= 0) then
777 0 : ABI_ERROR(sjoin("Cannot map k+q to IBZ with k+q:", ktoa(kq)))
778 : end if
779 69 : ikq_ibz = my_kqmap(1)
780 :
781 69 : ikq_ibz = my_kqmap(1); isym_kq = my_kqmap(2)
782 276 : trev_kq = my_kqmap(6); g0_kq = my_kqmap(3:5)
783 217 : isirr_kq = (isym_kq == 1 .and. trev_kq == 0 .and. all(g0_kq == 0))
784 276 : kq_ibz = ebands%kptns(:, ikq_ibz)
785 :
786 69 : if (dtset%eph_stern /= 0 .and. .not. sigma%imag_only) then
787 : ! Activate Sternheimer.
788 : ! NB: Assume adiabatic AHC expression to compute the contribution of states above gqk%nb.
789 :
790 : ! Get istwf_kq, npw_kq, kg_kq for k+q.
791 : call wfd%get_gvec_gbound(cryst%gmet, dtset%ecut, kq, ikq_ibz, isirr_kq, dtset%nloalg, & ! in
792 69 : istwf_kq, npw_kq, kg_kq, nkpg_kq, kpg_kq, gbound_kq) ! out
793 :
794 : call gs_ham_kq%eph_setup_k("kq", kq, istwfk_1, npw_kq, kg_kq, dtset, cryst, psps, & ! in
795 69 : nkpg, kpg_kq, ffnl_kq, kinpw_kq, ph3d_kq, xmpi_comm_self) ! out
796 :
797 : ! Fourier interpolation of the DFPT potentials.
798 69 : call dvdb%get_ftqbz(qpt, cplex, nfftf, ngfftf, v1scf, gqk%pert_comm%value)
799 :
800 : ! Build array with GS wavefunctions cg_kq at k+q to prepare call to dfpt_cgwf.
801 : !stern_use_cache = merge(.True., .False., dtset%eph_stern == 1)
802 69 : stern_use_cache = .False.
803 621 : fermie1_idir_ipert = zero ! FIXME: This is needed for metals.
804 : call stern%init(dtset, npw_k, npw_kq, nspinor, nb_kq, nb_kq, fermie1_idir_ipert, &
805 69 : stern_use_cache, work_ngfft, mpi_enreg, xmpi_comm_self)
806 :
807 207 : ABI_MALLOC(bra_kq, (2, npw_kq*nspinor))
808 938 : do ibsum_kq=1, nb_kq
809 : ! Reconstruct u_kq(G) from the IBZ image.
810 : call wfd%rotate_cg(ibsum_kq, ndat1, spin, kq_ibz, npw_kq, kg_kq, istwf_kq, &
811 869 : cryst, my_kqmap, gbound_kq, work_ngfft, work, bra_kq)
812 2993642 : stern%cgq(:,:,ibsum_kq) = bra_kq
813 : end do
814 69 : ABI_FREE(bra_kq)
815 :
816 : ! Loop over all 3*natom perturbations (Each core prepares its own potentials)
817 : ! In the inner loop, we calculate H1 * psi_k, stored in h1kets_kq on the k+q sphere.
818 : ! Allocate vlocal1 with correct cplex. Note nvloc
819 483 : ABI_MALLOC_OR_DIE(vlocal1, (cplex*n4, n5, n6, gs_ham_kq%nvloc, gqk%my_npert), ierr)
820 7760587 : ABI_CALLOC(cg1s_kq, (2, npw_kq*nspinor, natom3, nb_k))
821 :
822 : ! h1kets_kq are MPI distributed inside pert_comm but we need off-diagonal pp' terms --> collect results.
823 7760518 : ABI_CALLOC(h1kets_kq_allperts, (2, npw_kq*nspinor, natom3, nb_k))
824 :
825 483 : do my_ip=1, gqk%my_npert
826 : !print *, "my_ip:", my_ip
827 414 : ipc = gqk%my_pertcases(my_ip); idir = mod(ipc-1, 3) + 1; ipert = (ipc - idir) / 3 + 1
828 :
829 : ! Set up local potential vlocal1 with proper dimensioning, from vtrial1 taking into account the spin.
830 : ! Each CPU prepares its own potentials.
831 : call rf_transgrid_and_pack(spin, nspden, psps%usepaw, cplex, nfftf, nfft, ngfft, gs_ham_kq%nvloc, &
832 414 : pawfgr, mpi_enreg, vtrial, v1scf(:,:,:,my_ip), vlocal, vlocal1(:,:,:,:,my_ip))
833 :
834 : ! Continue to initialize the Hamiltonian (call it here to support dfpt_cgwf Sternheimer).
835 414 : call gs_ham_kq%load_spin(spin, vlocal=vlocal, with_nonlocal=.true.)
836 :
837 : ! Prepare application of the NL part.
838 414 : call rf_ham_kq%init(cplex, gs_ham_kq, ipert, has_e1kbsc=.true.)
839 414 : call rf_ham_kq%load_spin(spin, vlocal1=vlocal1(:,:,:,:,my_ip), with_nonlocal=.true.)
840 :
841 1656 : ABI_MALLOC(kets_k, (2, npw_k*nspinor, nb_k))
842 :
843 2706 : do in_k=1,nb_k
844 2292 : band_k = in_k + gqk%bstart_k - 1
845 30696 : stern%bands_treated_now(:) = 0; stern%bands_treated_now(band_k) = 1
846 30696 : stern%rank_band = 0; u1_band = band_k; band_me = band_k
847 :
848 : ! Init entry in cg1s_kq, either from cache or with zeros.
849 7759860 : cg1s_kq(:,:,ipc,in_k) = zero
850 2292 : call wfd%copy_cg(band_k, ik_ibz, spin, kets_k(1, 1, in_k))
851 :
852 : !print *, "Stern for band_k", band_k, " with nb_kq:", gqk%nb_kq
853 : ! NOTE: Here we are using the KS energies instead of the QP ones
854 : call stern%solve(u1_band, band_me, idir, ipert, qpt, gs_ham_kq, rf_ham_kq, &
855 : ks_ebands%eig(:,ik_ibz,spin), ks_ebands%eig(:,ikq_ibz,spin), &
856 2292 : kets_k(:,:,in_k), cwaveprj0, cg1s_kq(:,:,ipc,in_k), cwaveprj, msg, ierr)
857 2292 : ABI_CHECK(ierr == 0, msg)
858 :
859 : ! Store H(1) applied to GS wavefunction Psi_nk(0)
860 7760274 : h1kets_kq_allperts(:,:,ipc,in_k) = stern%gh1c_n
861 : end do ! in_k
862 :
863 414 : ABI_FREE(kets_k)
864 483 : call rf_ham_kq%free()
865 : end do ! my_ip (loop over my perturbations)
866 :
867 69 : call xmpi_sum(h1kets_kq_allperts, gqk%pert_comm%value, ierr)
868 69 : call xmpi_sum(cg1s_kq, gqk%pert_comm%value, ierr)
869 :
870 : ! Compute S_pp' = <D_{qp} vscf u_nk|u'_{nk+q p'}>
871 : ! Compute <D^q_p H psi_nk | D^q_p' psi_nk> and store it in stern_ppb
872 451 : do in_k=1,nb_k
873 : !call xmpi_allgather(h1kets_kq(:,:,:,in_k), 2*npw_kq*nspinor*gqk%my_npert, &
874 : ! h1kets_kq_allperts(:,:,:,in_k), gqk%pert_comm%value, ierr)
875 :
876 : call cg_zgemm("C", "N", npw_kq*nspinor, natom3, natom3, &
877 382 : h1kets_kq_allperts(:,:,:,in_k), cg1s_kq(:,:,:,in_k), stern_ppb(:,:,:,in_k))
878 :
879 : ! Save data for Debye-Waller
880 8203 : if (q_is_gamma) stern_dw(:,:,:,in_k) = stern_ppb(:,:,:,in_k)
881 : end do
882 :
883 69 : ABI_FREE(cg1s_kq)
884 69 : ABI_FREE(v1scf)
885 69 : ABI_FREE(vlocal1)
886 69 : ABI_FREE(h1kets_kq_allperts)
887 69 : call stern%free()
888 :
889 : ! Broadcast stern_dw here
890 138 : if (my_iq == 1) then
891 12 : if (gqk%qpt_comm%me == master) then
892 12 : ABI_CHECK(q_is_gamma, "Master in qpt_comm% should have Gamma when my_iq == 1")
893 : end if
894 12 : call xmpi_bcast(stern_dw, master, gqk%qpt_comm%value, ierr)
895 : end if
896 : end if ! eph_stern
897 :
898 : ! Sum over my phonon modes.
899 483 : do my_ip=1,gqk%my_npert
900 414 : nu = my_ip + gqk%my_pert_start - 1; wqnu = gqk%my_wnuq(my_ip, my_iq)
901 :
902 : ! Ignore unstable modes or modes that should be skipped.
903 414 : if (ephtk_skip_phmode(nu, wqnu, phmodes_skip, dtset%eph_phrange_w)) cycle
904 :
905 1026 : nqnu_t(:) = occ_be(wqnu, sigma%kTmesh, zero)
906 :
907 378 : if (dtset%prteliash /= 0) then
908 0 : dwargs = sigma%phmesh - wqnu
909 0 : dtw_weights(:, 1) = gaussian(dwargs, dtset%ph_smear)
910 : end if
911 :
912 378 : if (dtset%eph_stern /= 0) then
913 : ! Compute T_pp'(q,nu) matrix in reduced coordinates.
914 378 : call phdispl_cart2red_nmodes(natom, 1, cryst%gprimd, gqk%my_displ_cart(:,:,:,my_ip,my_iq), displ_red_nu)
915 378 : call sigtk_dw_tpp_red(natom, displ_red_nu, tpp_red)
916 : end if
917 :
918 : ! Sum over bands in |m,k+q>.
919 5217 : do im_kq=1,gqk%nb_kq
920 4770 : band_kq = im_kq + gqk%bstart_kq - 1
921 4770 : eig0mkq = ebands%eig(band_kq, ikq_ibz, spin)
922 4770 : eig0mk = ebands%eig(band_kq, ik_ibz, spin)
923 :
924 4770 : if (dtset%prteliash == 3) then
925 0 : delta_e_minus_emkq = gaussian(sigma%a2f_emesh - eig0mkq, dtset%tsmear)
926 : end if
927 :
928 : ! Compute electronic occupations for all Temps (note mu_e(it) Fermi level)
929 13590 : do it=1,ntemp
930 13590 : f_mkq(it) = occ_fd(eig0mkq, sigma%kTmesh(it), sigma%mu_e(it))
931 : end do
932 :
933 : ! Loop over the n band index in |n,k>.
934 31068 : do in_k=1,nb_k
935 25920 : band_k = in_k + gqk%bstart_k - 1
936 25920 : eig0nk = ebands%eig(band_k, ik_ibz, spin)
937 25920 : ediff = eig0nk - eig0mk
938 : intra_band = q_is_gamma .and. ediff <= TOL_EDIFF
939 25920 : same_band = band_k == band_kq
940 :
941 : ! The frequency dependent part evaluated at eig0nk for all T.
942 25920 : if (dtset%eph_ahc_type == 1) then
943 : cfact_t(:) = (nqnu_t + f_mkq ) / (eig0nk - eig0mkq + wqnu + sigma%ieta) + &
944 68040 : (nqnu_t - f_mkq + one) / (eig0nk - eig0mkq - wqnu + sigma%ieta)
945 : else
946 0 : cfact_t(:) = (two * nqnu_t + one) / (eig0nk - eig0mkq + sigma%ieta)
947 : end if
948 :
949 : ! Note the weight_q included in gkq2
950 25920 : if (with_cplex == 1) then
951 24120 : gkq2 = weight_q * gqk%my_g2(my_ip, im_kq, my_iq, in_k, my_ik)
952 :
953 : else
954 : gkq2 = weight_q * real(conjg(gqk%my_g_ks(my_ip, im_kq, my_iq, in_k, my_ik)) * &
955 1800 : gqk%my_g (my_ip, im_kq, my_iq, in_k, my_ik))
956 : end if
957 : !print *, "gkq2: ", gkq2
958 :
959 68040 : cfact_t = cfact_t * gkq2
960 :
961 : ! Compute contribution to Fan-Migdal for M > nb_kq
962 25920 : if (dtset%eph_stern /= 0 .and. im_kq == 1) then
963 : ! sum_{pp'} d_p* Stern_{pp'} d_p' with d = displ_red_nu and S = stern_ppb(:,:,:,in_k)
964 39672 : vec_natom3 = zero
965 2088 : call cg_zgemm("N", "N", natom3, natom3, 1, stern_ppb(:,:,:,in_k), displ_red_nu, vec_natom3)
966 2088 : dotri = cg_zdotc(natom3, displ_red_nu, vec_natom3)
967 : !write(std_out, *)"dotri:", dotri
968 2088 : rfact = dotri(1)
969 2088 : rfact = rfact * weight_q / (two * wqnu)
970 7344 : stern_fan_t = (two * nqnu_t(:) + one) * rfact
971 :
972 5256 : sigma%fan_stern_vals(:, in_k, ikcalc) = sigma%fan_stern_vals(:, in_k, ikcalc) + stern_fan_t
973 7344 : cfact_t = cfact_t + stern_fan_t
974 : end if
975 :
976 68040 : sigma%vals_e0ks(:, in_k, ikcalc) = sigma%vals_e0ks(:, in_k, ikcalc) + cfact_t
977 68040 : sigma%fan_vals(:, in_k, ikcalc) = sigma%fan_vals(:, in_k, ikcalc) + cfact_t
978 :
979 : ! Derivative of FM sigma at eig0nk for all T.
980 : ! Accumulate d(Re Sigma) / dw(w=eKS) for state in_k
981 : !cfact(x) = (nqnu_t + f_mkq ) / (x - eig0mkq + wqnu + sigma%ieta) + &
982 : ! (nqnu_t - f_mkq + one) / (x - eig0mkq - wqnu + sigma%ieta)
983 25920 : gmod2 = (eig0nk - eig0mkq + wqnu) ** 2
984 25920 : hmod2 = (eig0nk - eig0mkq - wqnu) ** 2
985 : rfact_t(:) = (nqnu_t + f_mkq ) * (-gmod2 + aimag(sigma%ieta)**2) / (gmod2 + aimag(sigma%ieta)**2) ** 2 + &
986 68040 : (nqnu_t - f_mkq + one) * (-hmod2 + aimag(sigma%ieta)**2) / (hmod2 + aimag(sigma%ieta)**2) ** 2
987 :
988 68040 : sigma%dvals_de0ks(:, in_k, ikcalc) = sigma%dvals_de0ks(:, in_k, ikcalc) + gkq2 * rfact_t
989 :
990 : ! Accumulate Sigma(w) for state |nk> if spectral function is wanted.
991 25920 : if (sigma%nwr > 0) then
992 : ! Zcut version
993 18000 : do it=1,ntemp
994 : cfact_wr(:) = (nqnu_t(it) + f_mkq(it) ) / (sigma%wrmesh_b(:,in_k, ikcalc) - eig0mkq + wqnu + sigma%ieta) + &
995 2908800 : (nqnu_t(it) - f_mkq(it) + one) / (sigma%wrmesh_b(:,in_k, ikcalc) - eig0mkq - wqnu + sigma%ieta)
996 2908800 : cfact_wr(:) = gkq2 * cfact_wr(:)
997 :
998 : !if (intra_band .and. sigma%frohl_model == 1) then
999 : ! ! Add Frohlich correction to Sigma_nk(w)
1000 : ! cfact_wr(:) = zero; if (same_band) cfact_wr(:) = fmw_frohl_sphcorr(:,nu,it,in_k)
1001 : !end if
1002 :
1003 2912400 : sigma%vals_wr(:,it,in_k,ikcalc) = sigma%vals_wr(:,it,in_k,ikcalc) + cfact_wr(:)
1004 :
1005 : ! Add static term from Sternheimer to Sigma(w) as well.
1006 : !if (dtset%eph_stern /= 0) then
1007 : ! !sigma%vals_wr(:,it,in_k,ikcalc) = sigma%vals_wr(:,it,in_k,ikcalc) + rtmp
1008 : !end if
1009 : end do
1010 : end if ! nwr > 0
1011 :
1012 25920 : gdw2 = gqk%my_gdw2(my_ip, im_kq, my_iq, in_k, my_ik)
1013 :
1014 : ! Accumulate DW for each T, add it to Sigma(e0) and Sigma(w) as well
1015 : ! - (2 n_{q\nu} + 1) * gdw2 / (e_nk - e_mk)
1016 25920 : if (abs(ediff) > EPHTK_WTOL) then
1017 55008 : cfact_t(:) = - weight_q * gdw2 * (two * nqnu_t + one) / (ediff + sigma%ieta)
1018 : else
1019 13032 : cfact_t(:) = zero
1020 : end if
1021 :
1022 25920 : if (dtset%eph_stern /= 0 .and. im_kq == 1) then
1023 : ! Compute DW term for M > nb_kq.
1024 2088 : cfact = zero
1025 14616 : do ip2=1,natom3
1026 89784 : do ip1=1,natom3
1027 87696 : cfact = cfact + tpp_red(ip1, ip2) * cmplx(stern_dw(1,ip1,ip2,in_k), stern_dw(2,ip1,ip2,in_k), kind=dp)
1028 : end do
1029 : end do
1030 : ! There's no 1/two here because I don't symmetrize the expression.
1031 : ! TODO: Test symmetrization, real quantity? add support for the different Eliashberg functions with Stern
1032 2088 : gdw2_stern = real(cfact) / (four * wqnu)
1033 :
1034 : ! Add contribution due to the Sternheimer. ediff is absorbed in Sternheimer.
1035 7344 : cfact2_t = - weight_q * gdw2_stern * (two * nqnu_t(:) + one)
1036 7344 : cfact_t = cfact_t + cfact2_t
1037 5256 : sigma%dw_stern_vals(:, in_k, ikcalc) = sigma%dw_stern_vals(:, in_k, ikcalc) + real(cfact2_t)
1038 : end if
1039 :
1040 68040 : sigma%dw_vals(:, in_k, ikcalc) = sigma%dw_vals(:, in_k, ikcalc) + real(cfact_t)
1041 68040 : sigma%vals_e0ks(:, in_k, ikcalc) = sigma%vals_e0ks(:, in_k, ikcalc) + real(cfact_t)
1042 :
1043 25920 : if (sigma%nwr > 0) then
1044 : ! Add static DW term to Sigma(w).
1045 18000 : do it=1,ntemp
1046 2912400 : sigma%vals_wr(:, it, in_k, ikcalc) = sigma%vals_wr(:, it, in_k, ikcalc) + real(cfact_t(it))
1047 : end do
1048 : end if
1049 :
1050 : ! Optionally, accumulate contribution to Eliashberg functions
1051 30690 : if (dtset%prteliash /= 0) then
1052 : ! EPH strength with delta(e_{nk} - e_{m\kq})
1053 : !rfact = gaussian(eig0nk - eig0mkq, dtset%tsmear)
1054 : !sigma%gf_nnuq(in_k, nu, iq_ibz_k, 1) = sigma%gf_nnuq(in_k, nu, iq_ibz_k, 1) + &
1055 : ! rfact * (gkq_nu(1, in_k, nu) ** 2 + gkq_nu(2, in_k, nu) ** 2)
1056 :
1057 : !! Treat contribution to Eliashberg function due to Fan term.
1058 : !if (ediff > wqnu) then
1059 : ! rfact = one / ediff
1060 : !else
1061 : ! ! Non adiabatic regime --> Add complex shift.
1062 : ! ! Note however that the expression for this flavor of Eliashberg function relies on adiabaticity.
1063 : ! rfact = real(one / (ediff + sigma%ieta))
1064 : !end if
1065 :
1066 : !gf_val = gkq_nu(1, in_k, nu) ** 2 + gkq_nu(2, in_k, nu) ** 2
1067 : !if (intra_band .and. sigma%frohl_model == 1) then
1068 : ! gf_val = zero; if (same_band) gf_val = zpr_frohl_sphcorr(nu) * (four_pi / three * q0rad ** 3)
1069 : !end if
1070 :
1071 : !sigma%gf_nnuq(in_k, nu, iq_ibz_k, 2) = sigma%gf_nnuq(in_k, nu, iq_ibz_k, 2) + gf_val * rfact
1072 : ! TODO: Add Sternheimer contribution
1073 :
1074 :
1075 : ! Optionally, accumulate DW contribution to Eliashberg functions.
1076 : !if (abs(ediff) > EPHTK_WTOL) then
1077 : ! sigma%gf_nnuq(ib_k, nu, iq_ibz_k, 3) = sigma%gf_nnuq(ib_k, nu, iq_ibz_k, 3) - gdw2 / ediff
1078 : !end if
1079 :
1080 0 : if (dtset%prteliash == 3) then
1081 : ! Accumulate: |g(k,q)|^2 delta(e - e_{m\kq}) delta(w - w_\qnu}
1082 0 : do iw=1,sigma%phmesh_size
1083 : sigma%a2few(:, iw, in_k, ikcalc) = sigma%a2few(:, iw, in_k, ikcalc) + &
1084 0 : delta_e_minus_emkq(:) * dtw_weights(iw, 1) * gkq2
1085 : end do
1086 : end if
1087 : end if ! prteliash /= 0
1088 :
1089 : end do ! in_k
1090 : end do ! im_kq
1091 : end do ! my_ip
1092 :
1093 69 : ABI_SFREE(kpg_kq)
1094 69 : ABI_SFREE(ffnl_kq)
1095 69 : ABI_SFREE(kinpw_kq)
1096 69 : ABI_SFREE(ph3d_kq)
1097 :
1098 81 : if (print_time_qq) then
1099 69 : call inds2str(2, "My q-point", my_iq, gqk%my_nq, gqk%glob_nq, msg)
1100 69 : call cwtime_report(msg, cpu_qq, wall_qq, gflops_qq); if (my_iq == LOG_MODQ) call wrtout(std_out, "...", do_flush=.True.)
1101 : end if
1102 : end do ! my_iq
1103 :
1104 12 : ABI_SFREE(kpg_k)
1105 12 : ABI_SFREE(ffnl_k)
1106 12 : ABI_SFREE(kinpw_k)
1107 12 : ABI_SFREE(ph3d_k)
1108 12 : call lg_myk%free()
1109 :
1110 19 : if (print_time_kk) then
1111 12 : call inds2str(3, "My k-point", my_ik, gqk%my_nk, gqk%glob_nk, msg)
1112 12 : call cwtime_report(msg, cpu_kk, wall_kk, gflops_kk); if (my_ik == LOG_MODK) call wrtout(std_out, "...", do_flush=.True.)
1113 : end if
1114 : end do ! my_ik
1115 :
1116 7 : ABI_SFREE(stern_ppb)
1117 7 : ABI_SFREE(stern_dw)
1118 7 : ABI_SFREE(cfact_wr)
1119 :
1120 14 : call sigma%gather_and_write_results(ncid, gstore, gqk, dtset, ebands)
1121 : end associate
1122 : end do ! my_is
1123 :
1124 7 : call cwtime_report(" gstore_sigmaph full calculation", cpu_all, wall_all, gflops_all, end_str=ch10)
1125 :
1126 7 : ABI_FREE(nqnu_t)
1127 7 : ABI_FREE(f_mkq)
1128 7 : ABI_FREE(cfact_t)
1129 7 : ABI_FREE(cfact2_t)
1130 7 : ABI_FREE(rfact_t)
1131 7 : ABI_FREE(stern_fan_t)
1132 7 : ABI_FREE(phmodes_skip)
1133 7 : ABI_FREE(ph1d)
1134 7 : ABI_SFREE(vtrial)
1135 : ABI_SFREE(cfact_wr)
1136 : ABI_SFREE(vtrial)
1137 7 : ABI_SFREE(vlocal)
1138 7 : ABI_SFREE(kg_kq)
1139 7 : ABI_SFREE(gbound_kq)
1140 7 : ABI_SFREE(tpp_red)
1141 7 : ABI_SFREE(work)
1142 7 : ABI_SFREE(dtw_weights)
1143 7 : ABI_SFREE(dwargs)
1144 7 : ABI_SFREE(delta_e_minus_emkq)
1145 :
1146 21 : call wfd%free(); call gstore%free(); call sigma%free(); call gs_ham_kq%free()
1147 :
1148 : contains
1149 :
1150 150 : subroutine inds2str(level, prefix, my_ik, my_nk, nk_tot, out_str)
1151 : character(len=*),intent(in) :: prefix
1152 : integer,intent(in) :: level, my_ik, my_nk, nk_tot
1153 : character(len=*),intent(out) :: out_str
1154 :
1155 150 : out_str = sjoin(prefix, itoa(my_ik), "/", itoa(my_nk), "[", itoa(nk_tot), "]")
1156 846 : out_str = repeat(' ', 4 * level) // trim(out_str)
1157 150 : end subroutine inds2str
1158 :
1159 : end subroutine gstore_sigmaph
1160 : !!***
1161 :
1162 : !----------------------------------------------------------------------
1163 :
1164 : !!****f* m_gstore_sigmaph/sep_gather_and_write_results
1165 : !! NAME
1166 : !! sep_gather_and_write_results
1167 : !!
1168 : !! FUNCTION
1169 : !! Collect results for a given spin, average results in the degenerate subspace.
1170 : !! Finally, write results to ab_out and netcdf file.
1171 : !!
1172 : !! INPUTS
1173 : !!
1174 : !! OUTPUT
1175 : !!
1176 : !! SOURCE
1177 :
1178 7 : subroutine sep_gather_and_write_results(sigma, root_ncid, gstore, gqk, dtset, ebands)
1179 :
1180 : !Arguments ------------------------------------
1181 : class(sep_t),intent(inout) :: sigma
1182 : integer,intent(in) :: root_ncid
1183 : type(gstore_t),intent(in) :: gstore
1184 : type(gqk_t),intent(in) :: gqk
1185 : type(ebands_t),intent(in) :: ebands
1186 : type(dataset_type),intent(in) :: dtset
1187 :
1188 : !Local variables-------------------------------
1189 : integer,parameter :: max_ntemp = 50, master = 0
1190 : integer :: it, in_k, ikcalc, ik_bz, spin, ierr, bstart_k, bstop_k, cnt, ndeg, spin_ncid, ncerr
1191 : integer :: band_k,ik_ibz,ib_val,ib_cond,jj,ideg,ii,iw, nstates !, nb_k
1192 : logical :: changed_k, iwrite
1193 : real(dp) :: ravg,kse,kse_prev,dw,fan0,ks_gap,kse_val,kse_cond,qpe_oms,qpe_oms_val,qpe_oms_cond
1194 : real(dp) :: ravg2 ! invsig2fmts, tau
1195 : complex(dp) :: sig0c,zc,qpe,qpe_prev,qpe_val,qpe_cond,cavg1,cavg2,cavg3,cavg4
1196 : character(len=500) :: this_gtype ! msg
1197 7 : type(degtab_t) :: degtab
1198 : !arrays
1199 7 : integer,allocatable :: degblock(:,:)
1200 : real(dp) :: kcalc(3)
1201 14 : real(dp) :: qp_gaps(sigma%ntemp),qpoms_gaps(sigma%ntemp), gfw_avg(sigma%phmesh_size, 3)
1202 7 : real(dp),allocatable :: aw(:,:,:), a2few_avg(:,:) !, gather_srate(:,:,:,:), grp_srate(:,:,:,:)
1203 14 : real(dp) :: ks_enes(gqk%nb_k), ze0_vals(sigma%ntemp, gqk%nb_k)
1204 14 : complex(dp) :: qpoms_enes(sigma%ntemp, gqk%nb_k),qp_enes(sigma%ntemp, gqk%nb_k)
1205 : !! *************************************************************************
1206 :
1207 7 : spin = gqk%spin
1208 :
1209 : !call wrtout(std_out, "Summing partial sigma terms inside qgk%comm.", do_flush=.True.)
1210 7 : call xmpi_sum(sigma%vals_e0ks, gqk%comm%value, ierr)
1211 7 : call xmpi_sum(sigma%dvals_de0ks, gqk%comm%value, ierr)
1212 7 : call xmpi_sum(sigma%fan_vals, gqk%comm%value, ierr)
1213 7 : call xmpi_sum(sigma%fan_stern_vals, gqk%comm%value, ierr)
1214 7 : call xmpi_sum(sigma%dw_vals, gqk%comm%value, ierr)
1215 7 : call xmpi_sum(sigma%dw_stern_vals, gqk%comm%value, ierr)
1216 7 : if (sigma%nwr > 0) call xmpi_sum(sigma%vals_wr, gqk%comm%value, ierr)
1217 : !call wrtout(std_out, "Sum completed.", do_flush=.True.)
1218 :
1219 : ! Only procs inside ncwrite_comm perform IO (ab_out and ncid)
1220 7 : iwrite = gqk%comm%me == 0; if (.not. iwrite) return
1221 :
1222 7 : this_gtype = "KS"
1223 7 : if (gstore%gtype == GSTORE_GTYPE_GWPT .and. dtset%gstore_gname == "gvals") this_gtype = "GWPT"
1224 7 : if (gstore%gtype == GSTORE_GTYPE_GWPT .and. dtset%gstore_gname == "gvals_ks") this_gtype = "KS"
1225 :
1226 : ! Write self-energy matrix elements for this spin
1227 : ! NB: Only master writes
1228 : ! Create hdf group for this spin.
1229 7 : NCF_CHECK(nf90_def_grp(root_ncid, strcat("gqk", "_spin", itoa(spin)), spin_ncid))
1230 : !NCF_CHECK(nf90_inq_ncid(root_ncid, strcat("data_spin", itoa(spin)), spin_ncid))
1231 :
1232 : ! Define dimensions and arrays inside group at runtime
1233 : ncerr = nctk_def_dims(spin_ncid, [ &
1234 : nctkdim_t("glob_nk", gqk%glob_nk), &
1235 : nctkdim_t("nb_kq", gqk%nb_kq), &
1236 : nctkdim_t("nb_k", gqk%nb_k) &
1237 28 : ], defmode=.True.)
1238 7 : NCF_CHECK(ncerr)
1239 :
1240 : ncerr = nctk_def_arrays(spin_ncid, [ &
1241 : nctkarr_t("vals_e0ks", "dp", "two, ntemp, nb_k, glob_nk"), &
1242 : nctkarr_t("fan_vals", "dp", "two, ntemp, nb_k, glob_nk"), &
1243 : nctkarr_t("fan_stern_vals", "dp", "two, ntemp, nb_k, glob_nk"), &
1244 : nctkarr_t("dvals_de0ks", "dp", "two, ntemp, nb_k, glob_nk"), &
1245 : nctkarr_t("dw_vals", "dp", "ntemp, nb_k, glob_nk"), &
1246 : nctkarr_t("dw_stern_vals", "dp", "ntemp, nb_k, glob_nk"), &
1247 : nctkarr_t("qpoms_enes", "dp", "two, ntemp, nb_k, glob_nk"), &
1248 : nctkarr_t("qp_enes", "dp", "two, ntemp, nb_k, glob_nk"), &
1249 : nctkarr_t("ze0_vals", "dp", "ntemp, nb_k, glob_nk"), &
1250 : nctkarr_t("ks_enes", "dp", "nb_k, glob_nk"), &
1251 : nctkarr_t("ks_gaps", "dp", "nb_k, glob_nk"), &
1252 : nctkarr_t("qpoms_gaps", "dp", "ntemp, nb_k, glob_nk"), &
1253 : nctkarr_t("qp_gaps", "dp", "ntemp, nb_k, glob_nk") &
1254 : !nctkarr_t("phmesh", "dp", "phmesh_size"), &
1255 : !nctkarr_t("vcar_calc", "dp", "three, max_nbcalc, glob_nk, nsppol") &
1256 98 : ])
1257 7 : NCF_CHECK(ncerr)
1258 :
1259 7 : if (sigma%nwr > 0) then
1260 : ! Make room for the spectral function.
1261 : ncerr = nctk_def_arrays(spin_ncid, [ &
1262 : nctkarr_t("wrmesh_b", "dp", "nwr, nb_k, glob_nk"), &
1263 : nctkarr_t("vals_wr", "dp", "two, nwr, ntemp, nb_k, glob_nk"), &
1264 : nctkarr_t("spfunc_wr", "dp", "nwr, ntemp, nb_k, glob_nk") &
1265 4 : ])
1266 1 : NCF_CHECK(ncerr)
1267 : end if
1268 :
1269 7 : if (dtset%prteliash /= 0) then
1270 : ncerr = nctk_def_arrays(spin_ncid, [ &
1271 : nctkarr_t("gfw_vals", "dp", "phmesh_size, three, nb_k, glob_nk") &
1272 0 : ])
1273 0 : NCF_CHECK(ncerr)
1274 0 : if (dtset%prteliash == 3) then
1275 : ncerr = nctk_def_arrays(spin_ncid, [ &
1276 : nctkarr_t("a2f_emesh", "dp", "a2f_ne"), &
1277 : nctkarr_t("a2few", "dp", "a2f_ne, phmesh_size, nb_k, glob_nk") &
1278 0 : ])
1279 0 : NCF_CHECK(ncerr)
1280 : end if
1281 : end if
1282 :
1283 : ! Write data.
1284 7 : NCF_CHECK(nctk_set_datamode(spin_ncid))
1285 7 : NCF_CHECK(nf90_put_var(spin_ncid, vid_spin("vals_e0ks"), c2r(sigma%vals_e0ks)))
1286 7 : NCF_CHECK(nf90_put_var(spin_ncid, vid_spin("fan_vals"), c2r(sigma%fan_vals)))
1287 7 : NCF_CHECK(nf90_put_var(spin_ncid, vid_spin("fan_stern_vals"), c2r(sigma%fan_stern_vals)))
1288 7 : NCF_CHECK(nf90_put_var(spin_ncid, vid_spin("dvals_de0ks"), c2r(sigma%dvals_de0ks)))
1289 7 : NCF_CHECK(nf90_put_var(spin_ncid, vid_spin("dw_vals"), sigma%dw_vals))
1290 7 : NCF_CHECK(nf90_put_var(spin_ncid, vid_spin("dw_stern_vals"), sigma%dw_stern_vals))
1291 7 : if (dtset%prteliash /=0) then
1292 0 : NCF_CHECK(nf90_put_var(spin_ncid, vid_spin("gfw_vals"), sigma%gfw_vals))
1293 : end if
1294 7 : if (dtset%prteliash == 3) then
1295 0 : NCF_CHECK(nf90_put_var(spin_ncid, vid_spin("a2f_emesh"), sigma%a2f_emesh))
1296 : end if
1297 :
1298 : ! Write legend.
1299 7 : if (spin == 1) then
1300 7 : write(ab_out,"(a)")repeat("=", 80)
1301 7 : write(ab_out,"(a)")" Final results in eV."
1302 7 : write(ab_out,"(a)")" Notations:"
1303 7 : write(ab_out,"(a)")" eKS: Kohn-Sham energy. eQP: quasi-particle energy."
1304 7 : write(ab_out,"(a)")" eQP - eKS: Difference between the QP and the KS energy."
1305 7 : write(ab_out,"(a)")" SE1(eKS): Real part of the self-energy computed at the KS energy, SE2 for imaginary part."
1306 7 : write(ab_out,"(a)")" Z(eKS): Renormalization factor."
1307 7 : write(ab_out,"(a)")" FAN: Real part of the Fan term at eKS. DW: Debye-Waller term."
1308 7 : write(ab_out,"(a)")" DeKS: KS energy difference between this band and band-1, DeQP same meaning but for eQP."
1309 7 : write(ab_out,"(a)")" OTMS: On-the-mass-shell approximation with eQP ~= eKS + Sigma(omega=eKS)"
1310 7 : write(ab_out,"(a)")" TAU(eKS): Lifetime in femtoseconds computed at the KS energy."
1311 7 : write(ab_out,"(a)")" mu_e: Fermi level for given (T, nelect)"
1312 7 : write(ab_out,"(a)")" "
1313 7 : write(ab_out,"(a)")" "
1314 7 : write(ab_out,"(2a)")" Using g(k,q) of type: ", trim(this_gtype)
1315 : !write(ab_out,"(2a)")" Treatment of gg: ", dtset%gwpt_g2mode
1316 7 : write(ab_out,"(a)")" "
1317 7 : write(ab_out,"(a)")" "
1318 : end if
1319 :
1320 : ! Compute QP energies and Gaps (Note that I'm assuming a non-magnetic semiconductor!)
1321 7 : ib_val = nint(ebands%nelect / (two / ebands%nspinor)); ib_cond = ib_val + 1
1322 :
1323 7 : if (sigma%a2f_ne > 0) then
1324 0 : ABI_MALLOC(a2few_avg, (sigma%a2f_ne, sigma%phmesh_size))
1325 : end if
1326 :
1327 19 : do ikcalc=1,gqk%glob_nk
1328 12 : ik_bz = gstore%kglob2bz(ikcalc, spin)
1329 12 : ik_ibz = gstore%kbz2ibz(1, ik_bz)
1330 48 : kcalc = gstore%kbz(:, ik_bz)
1331 :
1332 12 : if (dtset%symsigma == +1) then
1333 : ! Average self-energy matrix elements in the degenerate subspace.
1334 12 : bstart_k = gqk%bstart_k; bstop_k = gqk%bstop_k
1335 12 : call ebands%enclose_degbands(ik_ibz, spin, bstart_k, bstop_k, changed_k, dtset%symsigma_de, degblock=degblock)
1336 12 : bstart_k = gqk%bstart_k; bstop_k = gqk%bstop_k
1337 : !if (changed_k) then
1338 : ! ABI_WARNING("Changed")
1339 : !end if
1340 :
1341 : ! Store band indices used for averaging (shifted by bstart_k)
1342 12 : ndeg = size(degblock, dim=2)
1343 70 : ABI_MALLOC(degtab%bids, (ndeg))
1344 :
1345 46 : do ii=1,ndeg
1346 : ! Make sure boundaries are within the input nk states.
1347 : ! In principle the nk states should be initialized so that all degenerate states are included.
1348 34 : degblock(1, ii) = max(degblock(1, ii), bstart_k)
1349 34 : degblock(2, ii) = min(degblock(2, ii), bstop_k)
1350 34 : cnt = degblock(2, ii) - degblock(1, ii) + 1
1351 102 : ABI_MALLOC(degtab%bids(ii)%vals, (cnt))
1352 : degtab%bids(ii)%vals = [(jj, jj= &
1353 : degblock(1, ii) - bstart_k + 1, &
1354 284 : degblock(2, ii) - bstart_k + 1)]
1355 : end do
1356 :
1357 : ! Average self-energy matrix elements in the degenerate subspace.
1358 46 : do ideg=1,size(degtab%bids)
1359 12 : associate (bids => degtab%bids(ideg)%vals)
1360 34 : nstates = size(bids)
1361 :
1362 : ! Symmetrize Eliashberg functions
1363 34 : if (dtset%prteliash > 0) then
1364 0 : gfw_avg = sum(sigma%gfw_vals(:, :, bids(:), ikcalc), dim=3) / nstates
1365 0 : do ii=1,nstates
1366 0 : sigma%gfw_vals(:, :, bids(ii), ikcalc) = gfw_avg
1367 : end do
1368 0 : if (sigma%a2f_ne > 0) then
1369 0 : a2few_avg = sum(sigma%a2few(:, :, bids(:), ikcalc), dim=3) / nstates
1370 0 : do ii=1,nstates
1371 0 : sigma%a2few(:, :, bids(ii), ikcalc) = a2few_avg
1372 : end do
1373 : end if
1374 : end if
1375 :
1376 114 : do it=1,sigma%ntemp
1377 : ! Average QP(T) and Z(T).
1378 138 : cavg1 = sum(sigma%vals_e0ks(it, bids(:), ikcalc)) / nstates
1379 138 : cavg2 = sum(sigma%dvals_de0ks(it, bids(:), ikcalc)) / nstates
1380 138 : cavg3 = sum(sigma%fan_vals(it, bids(:), ikcalc)) / nstates
1381 138 : cavg4 = sum(sigma%fan_stern_vals(it, bids(:), ikcalc)) / nstates
1382 138 : ravg = sum(sigma%dw_vals(it, bids(:), ikcalc)) / nstates
1383 138 : ravg2 = sum(sigma%dw_stern_vals(it, bids(:), ikcalc)) / nstates
1384 :
1385 138 : do ii=1,nstates
1386 92 : sigma%vals_e0ks(it, bids(ii), ikcalc) = cavg1
1387 92 : sigma%dvals_de0ks(it, bids(ii), ikcalc) = cavg2
1388 92 : sigma%fan_vals(it, bids(ii), ikcalc) = cavg3
1389 92 : sigma%fan_stern_vals(it, bids(ii), ikcalc) = cavg4
1390 92 : sigma%dw_vals(it, bids(ii), ikcalc) = ravg
1391 138 : sigma%dw_stern_vals(it, bids(ii), ikcalc) = ravg2
1392 : end do ! ii
1393 :
1394 80 : if (sigma%nwr > 0) then
1395 : ! Average Sigma(omega, T)
1396 1616 : do iw=1,sigma%nwr
1397 4824 : cavg1 = sum(sigma%vals_wr(iw, it, bids(:), ikcalc)) / nstates
1398 4832 : do ii=1,nstates
1399 4824 : sigma%vals_wr(iw, it, bids(ii), ikcalc) = cavg1
1400 : end do
1401 : end do
1402 : end if
1403 : end do ! it
1404 : end associate
1405 : end do ! ideg
1406 :
1407 12 : call degtab%free()
1408 12 : ABI_FREE(degblock)
1409 : end if ! symsigma == +1
1410 :
1411 12 : kse_val = huge(one) * tol6; kse_cond = huge(one) * tol6
1412 332 : qp_enes = huge(one) * tol6; qpoms_enes = huge(one) * tol6
1413 240 : ks_enes = huge(one) * tol6; ze0_vals = huge(one) * tol6
1414 48 : ks_gap = -one; qpoms_gaps = -one; qp_gaps = -one
1415 :
1416 : ! Loop over temperatures.
1417 30 : do it=1,sigma%ntemp
1418 : ! Write header.
1419 18 : if (it <= max_ntemp) then
1420 18 : if (ebands%nsppol == 1) then
1421 : write(ab_out,"(3a,f6.1,a,f8.3)") &
1422 18 : "K-point: ", trim(ktoa(kcalc)), ", T: ", sigma%kTmesh(it) / kb_HaK, &
1423 36 : " [K], mu_e: ", sigma%mu_e(it) * Ha_eV
1424 : else
1425 : write(ab_out,"(3a,i1,a,f6.1,a,f8.3)") &
1426 0 : "K-point: ", trim(ktoa(kcalc)), ", spin: ", spin, ", T: ",sigma%kTmesh(it) / kb_HaK, &
1427 0 : " [K], mu_e: ", sigma%mu_e(it) * Ha_eV
1428 : end if
1429 18 : if (sigma%imag_only) then
1430 0 : write(ab_out,"(a)")" B eKS SE2(eKS) TAU(eKS) DeKS"
1431 : else
1432 18 : write(ab_out,"(a)")" B eKS eQP eQP-eKS SE1(eKS) SE2(eKS) Z(eKS) FAN(eKS) DW DeKS DeQP"
1433 : end if
1434 : end if
1435 :
1436 : ! Loop over band n_k for this k-point and spin.
1437 110 : do in_k=1,gqk%nb_k
1438 92 : band_k = in_k + bstart_k - 1
1439 92 : kse = ebands%eig(band_k, ik_ibz, spin)
1440 92 : ks_enes(in_k) = kse
1441 92 : sig0c = sigma%vals_e0ks(it, in_k, ikcalc)
1442 92 : dw = sigma%dw_vals(it, in_k, ikcalc)
1443 92 : fan0 = real(sig0c) - dw
1444 : ! Compute QP energies with On-the-Mass-Shell approximation and first renormalization i.e. Z(eKS)
1445 : ! TODO: Note that here I use the full Sigma including the imaginary part
1446 : !zc = one / (one - sigma%dvals_de0ks(it, in_k))
1447 92 : zc = one / (one - real(sigma%dvals_de0ks(it, in_k, ikcalc)))
1448 92 : ze0_vals(it, in_k) = real(zc)
1449 92 : qpe = kse + real(zc) * real(sig0c)
1450 92 : qpe_oms = kse + real(sig0c)
1451 92 : if (in_k == 1) then
1452 18 : kse_prev = kse; qpe_prev = qpe
1453 : end if
1454 92 : if (band_k == ib_val) then
1455 18 : kse_val = kse; qpe_val = qpe; qpe_oms_val = qpe_oms
1456 : end if
1457 92 : if (band_k == ib_cond) then
1458 18 : kse_cond = kse; qpe_cond = qpe; qpe_oms_cond = qpe_oms
1459 : end if
1460 :
1461 92 : if (it <= max_ntemp) then
1462 92 : if (sigma%imag_only) then
1463 : ! 1/tau = 2 Imag(Sigma)
1464 : !invsig2fmts = Time_Sec * 1e+15 / two
1465 : !tau = 999999.0_dp
1466 : !if (abs(aimag(sig0c)) > tol16) tau = invsig2fmts / abs(aimag(sig0c))
1467 : !tau = min(tau, 999999.0_dp)
1468 : !write(ab_out, "(i4,2(f8.3,1x),f8.1,1x,f8.3)") &
1469 : ! band_k, kse * Ha_eV, aimag(sig0c) * Ha_eV, tau, (kse - kse_prev) * Ha_eV
1470 : else
1471 : write(ab_out, "(i4, 10(f8.3,1x))") &
1472 92 : band_k, kse * Ha_eV, real(qpe) * Ha_eV, (real(qpe) - kse) * Ha_eV, &
1473 184 : real(sig0c) * Ha_eV, aimag(sig0c) * Ha_eV, real(zc), &
1474 184 : fan0 * Ha_eV, dw * Ha_eV, (kse - kse_prev) * Ha_eV, real(qpe - qpe_prev) * Ha_eV
1475 : end if
1476 : end if
1477 :
1478 92 : if (in_k > 1) then
1479 74 : kse_prev = kse; qpe_prev = qpe
1480 : end if
1481 92 : qpoms_enes(it, in_k) = qpe_oms
1482 92 : qp_enes(it, in_k) = qpe
1483 110 : if (kse_val /= huge(one) * tol6 .and. kse_cond /= huge(one) * tol6) then
1484 : ! We have enough states to compute the gap.
1485 56 : if (it == 1) ks_gap = kse_cond - kse_val
1486 56 : qpoms_gaps(it) = qpe_oms_cond - qpe_oms_val
1487 56 : qp_gaps(it) = real(qpe_cond - qpe_val)
1488 : end if
1489 : end do ! in_k
1490 :
1491 : ! Print KS and QP gaps.
1492 30 : if (it <= max_ntemp) then
1493 18 : if (.not. sigma%imag_only) then
1494 18 : if (kse_val /= huge(one) * tol6 .and. kse_cond /= huge(one) * tol6) then
1495 18 : write(ab_out, "(a)")" "
1496 18 : write(ab_out, "(a,f8.3,1x,2(a,i0),a)")" KS gap: ",ks_gap * Ha_eV, &
1497 36 : "(assuming bval:", ib_val, " ==> bcond:", ib_cond, ")"
1498 18 : write(ab_out, "(2(a,f8.3),a)")" QP gap: ",qp_gaps(it) * Ha_eV," (OTMS: ",qpoms_gaps(it) * Ha_eV, ")"
1499 18 : write(ab_out, "(2(a,f8.3),a)")" QP_gap - KS_gap: ",(qp_gaps(it) - ks_gap) * Ha_eV,&
1500 36 : " (OTMS: ",(qpoms_gaps(it) - ks_gap) * Ha_eV, ")"
1501 18 : write(ab_out, "(a)")" "
1502 : end if
1503 : else
1504 0 : if (kse_val /= huge(one) * tol6 .and. kse_cond /= huge(one) * tol6) then
1505 0 : write(ab_out, "(a)")" "
1506 0 : write(ab_out, "(a,f8.3,1x,2(a,i0),a)")" KS gap: ",ks_gap * Ha_eV, "(assuming bval:",ib_val," ==> bcond:",ib_cond,")"
1507 0 : write(ab_out, "(a)")" "
1508 : end if
1509 : end if
1510 18 : write(ab_out, "(a)")repeat("=", 92)
1511 : end if
1512 :
1513 : end do ! it
1514 :
1515 : ! Dump QP energies and gaps for this kpoint.
1516 60 : NCF_CHECK(nf90_put_var(spin_ncid, vid_spin("qpoms_enes"), c2r(qpoms_enes), start=[1,1,1,ikcalc]))
1517 60 : NCF_CHECK(nf90_put_var(spin_ncid, vid_spin("qp_enes"), c2r(qp_enes), start=[1,1,1,ikcalc]))
1518 48 : NCF_CHECK(nf90_put_var(spin_ncid, vid_spin("ze0_vals"), ze0_vals, start=[1,1,ikcalc]))
1519 36 : NCF_CHECK(nf90_put_var(spin_ncid, vid_spin("ks_enes"), ks_enes, start=[1,ikcalc]))
1520 36 : NCF_CHECK(nf90_put_var(spin_ncid, vid_spin("ks_gaps"), ks_gap, start=[1,ikcalc]))
1521 36 : NCF_CHECK(nf90_put_var(spin_ncid, vid_spin("qpoms_gaps"), qpoms_gaps, start=[1,ikcalc]))
1522 36 : NCF_CHECK(nf90_put_var(spin_ncid, vid_spin("qp_gaps"), qp_gaps, start=[1,ikcalc]))
1523 :
1524 12 : if (sigma%nwr > 0) then
1525 4 : NCF_CHECK(nf90_put_var(spin_ncid, vid_spin("wrmesh_b"), sigma%wrmesh_b(:,:,ikcalc), start=[1,1,ikcalc]))
1526 6 : NCF_CHECK(nf90_put_var(spin_ncid, vid_spin("vals_wr"), c2r(sigma%vals_wr(:,:,:,ikcalc)), start=[1,1,1,1,ikcalc]))
1527 :
1528 : ! Compute spectral function. A = -1/pi [Im Sigma(ww)] / ([ww - ee - Re Sigma(ww)] ** 2 + Im Sigma(ww) ** 2])
1529 5 : ABI_MALLOC(aw, (sigma%nwr, sigma%ntemp, gqk%nb_k))
1530 5 : do in_k=1,gqk%nb_k
1531 4 : band_k = in_k + bstart_k - 1
1532 4 : kse = ebands%eig(band_k, ik_ibz, spin)
1533 21 : do it=1,sigma%ntemp
1534 : aw(:, it, in_k) = -piinv * aimag(sigma%vals_wr(:, it, in_k, ikcalc)) / &
1535 3236 : ((sigma%wrmesh_b(:, in_k, ikcalc) - kse - real(sigma%vals_wr(:, it, in_k, ikcalc))) ** 2 + aimag(sigma%vals_wr(:, it, in_k, ikcalc)) ** 2)
1536 : end do
1537 : end do
1538 5 : NCF_CHECK(nf90_put_var(spin_ncid, vid_spin("spfunc_wr"), aw, start=[1, 1, 1, ikcalc]))
1539 1 : ABI_FREE(aw)
1540 : end if
1541 :
1542 19 : if (dtset%prtvol > 0 .and. (ikcalc == 1)) then
1543 1 : if (allocated(sigma%gfw_vals)) then
1544 0 : write(ab_out, "(2a)")" omega and Eliashberg function gf_{nk}(omega) for testing purposes:"
1545 0 : iw = (sigma%phmesh_size / 2)
1546 0 : do in_k=1,min(gqk%nb_k, 5)
1547 0 : band_k = in_k + bstart_k - 1
1548 0 : write(ab_out, "(a, i0)")"For band:", band_k
1549 0 : do jj=0,1
1550 0 : write(ab_out, "(4(f8.3,2x))")sigma%phmesh(iw+jj), (sigma%gfw_vals(iw+jj, ii, in_k, ikcalc), ii=1,3)
1551 : end do
1552 : end do
1553 0 : write(ab_out, "(a)")ch10
1554 : end if
1555 :
1556 1 : if (sigma%nwr >= 3) then
1557 1 : write(ab_out, "(2a)")ch10," omega and Sigma_nk(omega, T=1) in eV for testing purposes:"
1558 1 : it = 1; iw = (sigma%nwr / 2)
1559 5 : do in_k=1,min(gqk%nb_k, 5)
1560 4 : band_k = in_k + bstart_k - 1
1561 4 : write(ab_out, "(a, i0)")"For band:", band_k
1562 13 : do ii=0,1
1563 12 : write(ab_out, "(3(f8.3,2x))")sigma%wrmesh_b(iw+ii, in_k, ikcalc) * Ha_eV, sigma%vals_wr(iw+ii, it, in_k, ikcalc) * Ha_eV
1564 : end do
1565 : end do
1566 1 : write(ab_out, "(a)")ch10
1567 : end if
1568 : end if
1569 :
1570 : end do ! ikcalc
1571 :
1572 7 : if (sigma%ntemp > max_ntemp) then
1573 0 : write(ab_out, "(a,i0,a)")" No more than ", max_ntemp, " temperatures are written to the main output file."
1574 0 : write(ab_out, "(2a)")" Please use the GSEPH.nc file and AbiPy to analyze the results.",ch10
1575 : end if
1576 :
1577 14 : ABI_SFREE(a2few_avg)
1578 :
1579 : contains
1580 129 : integer function vid_spin(var_name)
1581 : character(len=*),intent(in) :: var_name
1582 129 : vid_spin = nctk_idname(spin_ncid, var_name)
1583 : end function vid_spin
1584 :
1585 : end subroutine sep_gather_and_write_results
1586 : !!***
1587 :
1588 : !----------------------------------------------------------------------
1589 :
1590 : !!****f* m_gstore_sigmaph/sep_free
1591 : !! NAME
1592 : !! sep_free
1593 : !!
1594 : !! FUNCTION
1595 : !! Free dynamic memory.
1596 : !!
1597 : !! SOURCE
1598 :
1599 7 : subroutine sep_free(sigma)
1600 :
1601 : !Arguments ------------------------------------
1602 : class(sep_t),intent(inout) :: sigma
1603 : ! *********************************************************************
1604 :
1605 7 : ABI_SFREE(sigma%kTmesh)
1606 7 : ABI_SFREE(sigma%mu_e)
1607 7 : ABI_SFREE(sigma%vals_e0ks)
1608 7 : ABI_SFREE(sigma%dvals_de0ks)
1609 7 : ABI_SFREE(sigma%fan_vals)
1610 7 : ABI_SFREE(sigma%fan_stern_vals)
1611 7 : ABI_SFREE(sigma%dw_vals)
1612 7 : ABI_SFREE(sigma%dw_stern_vals)
1613 7 : ABI_SFREE(sigma%vals_wr)
1614 7 : ABI_SFREE(sigma%wrmesh_b)
1615 7 : ABI_SFREE(sigma%phmesh)
1616 7 : ABI_SFREE(sigma%gfw_vals)
1617 7 : ABI_SFREE(sigma%a2f_emesh)
1618 7 : ABI_SFREE(sigma%a2few)
1619 :
1620 7 : end subroutine sep_free
1621 : !!***
1622 :
1623 0 : end module m_gstore_sigmaph
1624 : !!***
|
