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