Branch data Line data Source code
1 : : !!****m* ABINIT/m_migdal_eliashberg
2 : : !! NAME
3 : : !! m_migdal_eliashberg
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 2008-2026 ABINIT group (MG)
9 : : !! This file is distributed under the terms of the
10 : : !! GNU General Public License, see ~abinit/COPYING
11 : : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : : !!
13 : : !! SOURCE
14 : :
15 : : #if defined HAVE_CONFIG_H
16 : : #include "config.h"
17 : : #endif
18 : :
19 : : #include "abi_common.h"
20 : :
21 : : module m_migdal_eliashberg
22 : :
23 : : use defs_basis
24 : : use m_abicore
25 : : use m_xmpi
26 : : use m_errors
27 : : use m_krank
28 : : use m_htetra
29 : : use netcdf
30 : : use m_nctk
31 : : use m_crystal
32 : : use m_dtset
33 : : use m_dtfil
34 : :
35 : : use m_time, only : cwtime, cwtime_report, sec2str
36 : : use m_fstrings, only : strcat, sjoin, itoa, ftoa, ktoa, ltoa
37 : : use m_copy, only : alloc_copy
38 : : use m_numeric_tools, only : simpson_int
39 : : use m_special_funcs, only : gaussian
40 : : use m_ebands, only : ebands_t, edos_t
41 : : use m_bz_mesh, only : kpath_t
42 : : use m_geometry, only : phdispl_cart2red_nmodes
43 : : use m_ephtk, only : ephtk_gkknu_from_atm, EPHTK_WTOL
44 : : use m_gstore, only : gstore_t
45 : :
46 : : implicit none
47 : :
48 : : private
49 : :
50 : : public :: migdal_eliashberg_iso
51 : : !public :: migdal_eliashberg_aniso
52 : :
53 : : ! Internal policy for mode-resolved lambda on a q-path. When enabled, all
54 : : ! modes in a degenerate phonon subspace are assigned the subspace-averaged
55 : : ! lambda. The sum over the subspace is therefore preserved.
56 : : logical,parameter :: average_lambda_qpath_degenerate = .True.
57 : : real(dp),parameter :: lambda_qpath_degen_tol = tol6
58 : : !!***
59 : :
60 : : !----------------------------------------------------------------------
61 : :
62 : : !!****t* m_migdal_eliashberg/iso_solver_t
63 : : !! NAME
64 : : !! iso_solver_t
65 : : !!
66 : : !! FUNCTION
67 : : !!
68 : : !! NOTES
69 : : !!
70 : : !! SOURCE
71 : :
72 : : type, public :: iso_solver_t
73 : :
74 : : integer :: ntemp
75 : :
76 : : integer :: max_niter = -1
77 : : ! Maximum number of iterations.
78 : :
79 : : integer :: max_nmix = 4
80 : :
81 : : integer :: ncid
82 : :
83 : : integer :: comm = xmpi_undefined
84 : :
85 : : !integer :: niw = -1
86 : : ! Number of Matsubara frequencies.
87 : :
88 : : !real(dp) :: kt = -one
89 : : ! K * T in Ha.
90 : :
91 : : real(dp) :: tolerance = -one
92 : :
93 : : real(dp),allocatable :: zeta_iw(:), delta_iw(:)
94 : : real(dp),allocatable :: prev_zeta_iw(:), prev_delta_iw(:)
95 : : real(dp),allocatable :: delta_iw_mix(:,:)
96 : :
97 : : contains
98 : :
99 : : procedure :: free => iso_solver_free
100 : : ! Free dynamic memory
101 : :
102 : : procedure :: solve => iso_solver_solve
103 : :
104 : : end type iso_solver_t
105 : : !!***
106 : :
107 : : contains
108 : : !!***
109 : :
110 : :
111 : : !----------------------------------------------------------------------
112 : :
113 : : !!****f* m_migdal_eliashberg/migdal_eliashberg_iso
114 : : !! NAME
115 : : !! migdal_eliashber_iso
116 : : !!
117 : : !! FUNCTION
118 : : !!
119 : : !! INPUTS
120 : : !!
121 : : !! OUTPUT
122 : : !!
123 : : !! SOURCE
124 : :
125 : 0 : subroutine migdal_eliashberg_iso(gstore, dtset, dtfil)
126 : :
127 : : !Arguments ------------------------------------
128 : : !scalars
129 : : type(dataset_type),intent(in) :: dtset
130 : : type(datafiles_type),intent(in) :: dtfil
131 : : type(gstore_t),target,intent(inout) :: gstore
132 : :
133 : : !Local variables-------------------------------
134 : : !scalars
135 : : integer,parameter :: master = 0
136 : : integer :: nproc, my_rank, ierr, itemp, ntemp, niw, ncid, iw, spin
137 : : integer :: edos_intmeth
138 : : !integer :: spin, natom3, cnt !, band, ib, nb, my_ik, my_iq, my_is
139 : : !integer :: ik_ibz, ik_bz, ebands_timrev, iq_bz, iq_ibz !, ikq_ibz, ikq_bz
140 : : !integer :: ncid, spin_ncid, ncerr, gstore_fform
141 : : integer :: phmesh_size, units(2) !, iw
142 : : real(dp) :: kt, wmax, cpu, wall, gflops, edos_step, edos_broad, lambda_iso, omega_log, omega_2, alpha !, sigma, ecut, eshift, eig0nk
143 : : character(len=500) :: msg
144 : : class(crystal_t),pointer :: cryst
145 : : class(ebands_t),pointer :: ebands
146 : 0 : type(iso_solver_t) :: iso
147 : 0 : type(edos_t) :: edos
148 : : !arrays
149 : 0 : real(dp),allocatable :: ktmesh(:), lambda_ij(:), imag_w(:), imag_2w(:), phmesh(:), a2fw(:), a2fw_raw(:)
150 : 0 : real(dp),allocatable :: a2f_1mom(:), a2f_1mom_int(:)
151 : 0 : real(dp),allocatable :: phfreq_qibz(:,:), phlambda_qibz(:,:,:)
152 : 0 : real(dp),allocatable :: qpath(:,:), phfreq_qpath(:,:), phdispl_cart_qpath(:,:,:,:), phlambda_qpath(:,:,:)
153 : : !----------------------------------------------------------------------
154 : :
155 : 0 : nproc = xmpi_comm_size(gstore%comm); my_rank = xmpi_comm_rank(gstore%comm)
156 [ # # ]: 0 : units = [std_out, ab_out]
157 : :
158 : 0 : call cwtime(cpu, wall, gflops, "start")
159 : 0 : call wrtout(units, " Solving isotropic Migdal-Eliashberg equations on the imaginary axis", pre_newlines=2)
160 : :
161 : 0 : cryst => gstore%cryst; ebands => gstore%ebands
162 : : !natom3 = 3 * cryst%natom; nsppol = ebands%nsppol
163 : :
164 : : ! Consistency check
165 : 0 : ierr = 0
166 [ # # ]: 0 : ABI_CHECK_NOSTOP(gstore%kzone == "bz", "gstore_kzone == 'bz' is required", ierr)
167 [ # # # # : 0 : ABI_CHECK_NOSTOP(any(gstore%qzone == ["bz ", "ibz"]), "gstore_qzone must be 'bz' or 'ibz'", ierr)
# # ]
168 [ # # ]: 0 : ABI_CHECK(ierr == 0, "Wrong gstore object for migdal_eliashberg_iso. See messages above")
169 : :
170 : : ! Compute electron DOS.
171 : 0 : call dtset%get_edos_params(edos_intmeth, edos_step, edos_broad)
172 : 0 : edos = ebands%get_edos(cryst, edos_intmeth, edos_step, edos_broad, gstore%comm)
173 : :
174 : : ! A disentangled Wannier Hamiltonian generally spans only a subspace of the
175 : : ! original bands. In this case ebands%nelect still describes the complete
176 : : ! ab-initio manifold, so locating eF by integrating the DOS can fail (e.g. a
177 : : ! four-band model with nelect = 8 appears completely filled). The electronic
178 : : ! delta functions used below are centered at ebands%fermie, copied from the
179 : : ! ab-initio bands, hence evaluate N(eF) at the same chemical potential.
180 [ # # # # ]: 0 : if (edos%ief == 0 .and. gstore%has_wannier) then
181 : 0 : iw = int((ebands%fermie - edos%mesh(1)) / edos%step) + 1
182 [ # # # # ]: 0 : ABI_CHECK(iw >= 1 .and. iw < edos%nw, "The ab-initio Fermi level lies outside the energy range of the Wannier-interpolated bands")
183 : 0 : alpha = (ebands%fermie - edos%mesh(iw)) / edos%step
184 [ # # ]: 0 : do spin=0,edos%nsppol
185 : 0 : edos%gef(spin) = (one - alpha) * edos%dos(iw,spin) + alpha * edos%dos(iw+1,spin)
186 : 0 : edos%ghf(spin) = edos%gef(spin)
187 : : end do
188 : 0 : edos%ief = iw
189 : 0 : edos%ihf = iw
190 : 0 : call wrtout(units, " Using the ab-initio Fermi level to evaluate the DOS of the Wannier band subspace.")
191 : : end if
192 : :
193 : : !! Store DOS per spin channel
194 : : !n0(:) = edos%gef(1:edos%nsppol)
195 [ # # ]: 0 : if (my_rank == master) then
196 : 0 : call edos%print(units)
197 : : !path = strcat(dtfil%filnam_ds(4), "_EDOS")
198 : : !call wrtout(ab_out, sjoin("- Writing electron DOS to file:", path, ch10))
199 : : !call edos%write(path)
200 : : end if
201 : :
202 : : ! Compute phonon frequency mesh.
203 : 0 : call gstore%ifc%get_phmesh(dtset%ph_wstep, phmesh_size, phmesh)
204 : :
205 : : ! Compute Eliashberg function a2F(w)
206 [ # # # # : 0 : ABI_MALLOC(a2fw, (phmesh_size))
# # # # #
# # # ]
207 : 0 : call get_a2fw(gstore, edos%gef(0), phmesh_size, phmesh, a2fw, phfreq_qibz, phlambda_qibz)
208 : :
209 : : ! Compute mode-resolved lambda on the phonon q-path when the matrix
210 : : ! elements can be evaluated at arbitrary q with Wannier interpolation.
211 : 0 : call get_lambda_qpath_wan(gstore, dtset, edos%gef(0), qpath, phfreq_qpath, phdispl_cart_qpath, phlambda_qpath)
212 : :
213 : : ! Save results.
214 : 0 : ncid = nctk_noid
215 [ # # ]: 0 : if (my_rank == master) then
216 : 0 : call alloc_copy(a2fw, a2fw_raw)
217 [ # # ]: 0 : ABI_CHECK(edos%gef(0) > zero, "The electronic DOS at the Fermi level must be positive")
218 [ # # ]: 0 : a2fw = a2fw / (edos%gef(0) / two)
219 : :
220 [ # # # # : 0 : ABI_MALLOC(a2f_1mom, (phmesh_size))
# # # # ]
221 [ # # # # : 0 : ABI_MALLOC(a2f_1mom_int, (phmesh_size))
# # # # ]
222 [ # # ]: 0 : a2f_1mom = zero
223 [ # # # # ]: 0 : where (phmesh > tol12) a2f_1mom = a2fw / phmesh
224 : 0 : call simpson_int(phmesh_size, dtset%ph_wstep, a2f_1mom, a2f_1mom_int)
225 : 0 : lambda_iso = two * a2f_1mom_int(phmesh_size)
226 : 0 : write(msg, "(a,es16.8)")" Isotropic lambda from a2F(w): ", lambda_iso
227 : 0 : call wrtout(units, msg)
228 : :
229 [ # # ]: 0 : ABI_CHECK(lambda_iso > zero, "Cannot compute omega_log because the isotropic lambda is not positive")
230 [ # # ]: 0 : a2f_1mom = zero
231 [ # # ]: 0 : do iw=1,phmesh_size
232 [ # # ]: 0 : if (phmesh(iw) > tol12) a2f_1mom(iw) = a2fw(iw) * log(phmesh(iw)) / phmesh(iw)
233 : : end do
234 : 0 : call simpson_int(phmesh_size, dtset%ph_wstep, a2f_1mom, a2f_1mom_int)
235 : 0 : omega_log = exp(two * a2f_1mom_int(phmesh_size) / lambda_iso)
236 : 0 : write(msg, "(a,es16.8,a,es16.8,a)")" Isotropic omega_log from a2F(w): ", omega_log, &
237 : 0 : " (Ha), ", omega_log * Ha_K, " (K)"
238 : 0 : call wrtout(units, msg)
239 : :
240 : : ! Allen-Dynes square-root second moment:
241 : : ! omega_2^2 = (2 / lambda) integral dw w a2F(w).
242 [ # # # # : 0 : a2f_1mom = phmesh * a2fw
# # # # ]
243 : 0 : call simpson_int(phmesh_size, dtset%ph_wstep, a2f_1mom, a2f_1mom_int)
244 : 0 : omega_2 = sqrt(two * a2f_1mom_int(phmesh_size) / lambda_iso)
245 : 0 : write(msg, "(a,es16.8,a,es16.8,a)")" Isotropic omega_2 from a2F(w): ", omega_2, &
246 : 0 : " (Ha), ", omega_2 * Ha_K, " (K)"
247 : 0 : call wrtout(units, msg)
248 [ # # ]: 0 : ABI_FREE(a2f_1mom)
249 : 0 : ABI_FREE(a2f_1mom_int)
250 : :
251 [ # # ]: 0 : NCF_CHECK(nctk_open_create(ncid, strcat(dtfil%filnam_ds(4), "_ISOME.nc") , xmpi_comm_self))
252 [ # # ]: 0 : NCF_CHECK(cryst%ncwrite(ncid))
253 [ # # ]: 0 : NCF_CHECK(ebands%ncwrite(ncid))
254 [ # # ]: 0 : NCF_CHECK(edos%ncwrite(ncid))
255 : 0 : call isome_ncwrite_spectral(ncid, dtset, gstore, phmesh_size, phmesh, a2fw_raw, a2fw, edos%gef(0), omega_2)
256 : 0 : call isome_ncwrite_qibz(ncid, gstore, phfreq_qibz, phlambda_qibz)
257 [ # # ]: 0 : if (allocated(qpath)) then
258 : 0 : call isome_ncwrite_qpath(ncid, qpath, phfreq_qpath, phdispl_cart_qpath, phlambda_qpath)
259 : : end if
260 [ # # ]: 0 : ABI_FREE(a2fw_raw)
261 : : end if
262 : :
263 [ # # ]: 0 : ABI_SFREE(qpath)
264 [ # # ]: 0 : ABI_SFREE(phfreq_qpath)
265 [ # # ]: 0 : ABI_SFREE(phdispl_cart_qpath)
266 [ # # ]: 0 : ABI_SFREE(phlambda_qpath)
267 [ # # ]: 0 : ABI_FREE(phfreq_qibz)
268 [ # # ]: 0 : ABI_FREE(phlambda_qibz)
269 [ # # ]: 0 : ABI_FREE(a2fw)
270 [ # # ]: 0 : ABI_FREE(phmesh)
271 : 0 : call edos%free()
272 : :
273 : 0 : call dtset%get_ktmesh(ntemp, ktmesh)
274 : :
275 : : !NVHPC and LLVM don't like using this constructor because allocatable arrays aren't set.
276 : : #if defined FC_NVHPC || defined FC_LLVM
277 : : iso%ntemp=ntemp
278 : : iso%max_niter=10
279 : : iso%tolerance=tol10
280 : : iso%ncid=ncid
281 : : iso%comm=gstore%comm
282 : : #else
283 : 0 : iso = iso_solver_t(ntemp=ntemp, max_niter=10, tolerance=tol10, ncid=ncid, comm=gstore%comm)
284 : : #endif
285 : :
286 [ # # ]: 0 : do itemp=1,ntemp
287 : : ! Generate Matsubara mesh for this T with cutoff wmax.
288 : 0 : kt = ktmesh(itemp)
289 : 0 : wmax = one
290 [ # # ]: 0 : call matsubara_mesh("bosons", kt, wmax, niw, imag_w)
291 : :
292 : : ! Compute lambda(w_i - w_j)
293 [ # # # # : 0 : ABI_MALLOC(lambda_ij, (2 * niw))
# # # # #
# # # #
# ]
294 [ # # # # : 0 : ABI_MALLOC(imag_2w, (2 * niw))
# # # # ]
295 : :
296 : : !call wrtout(std_out, " Computing lambda_iso_iw...")
297 : : !call get_lambda_iso_iw(gstore, 2 * niw, imag_2w, lambda_ij)
298 [ # # ]: 0 : ABI_FREE(imag_2w)
299 : :
300 : : !call iso%solve(itemp, kt, niw, imag_w, lambda_ij)
301 [ # # ]: 0 : ABI_FREE(lambda_ij)
302 [ # # ]: 0 : ABI_FREE(imag_w)
303 : : end do ! itemp
304 : :
305 [ # # ]: 0 : if (my_rank == master) then
306 [ # # ]: 0 : NCF_CHECK(nf90_close(ncid))
307 : : end if
308 : :
309 [ # # ]: 0 : ABI_FREE(ktmesh)
310 : 0 : call iso%free()
311 : :
312 : 0 : call cwtime_report(" migdal_eliashberg_iso:", cpu, wall, gflops)
313 : :
314 [ # # # # : 0 : end subroutine migdal_eliashberg_iso
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # ]
315 : : !!***
316 : :
317 : : !----------------------------------------------------------------------
318 : :
319 : : !!****f* m_migdal_eliashberg/isome_ncwrite_spectral
320 : : !! NAME
321 : : !! isome_ncwrite_spectral
322 : : !!
323 : : !! FUNCTION
324 : : !! Write the first revision of the spectral section of the ISOME.nc file.
325 : : !! The raw matrix-element spectral sum is stored together with the
326 : : !! provisionally DOS-normalized Eliashberg function so that the normalization
327 : : !! can be audited without recomputing the electron-phonon matrix elements.
328 : : !!
329 : : !! INPUTS
330 : : !! ncid=NetCDF file identifier, opened on xmpi_comm_self by the master rank.
331 : : !! dtset<dataset_type>=Input variables.
332 : : !! gstore<gstore_t>=Electron-phonon matrix-element container.
333 : : !! nomega=Number of points in the phonon-frequency mesh.
334 : : !! omega(nomega)=Phonon-frequency mesh in Hartree.
335 : : !! a2f_raw(nomega)=Raw matrix-element spectral sum.
336 : : !! a2f(nomega)=DOS-normalized Eliashberg function.
337 : : !! edos_fermie=Total electronic DOS at the Fermi level.
338 : : !! omega_2=Allen-Dynes square-root second moment in Hartree.
339 : : !!
340 : : !! SOURCE
341 : :
342 : 0 : subroutine isome_ncwrite_spectral(ncid, dtset, gstore, nomega, omega, a2f_raw, a2f, edos_fermie, omega_2)
343 : :
344 : : !Arguments ------------------------------------
345 : : !scalars
346 : : integer,intent(in) :: ncid, nomega
347 : : real(dp),intent(in) :: edos_fermie, omega_2
348 : : type(dataset_type),intent(in) :: dtset
349 : : type(gstore_t),intent(in) :: gstore
350 : : !arrays
351 : : real(dp),intent(in) :: omega(nomega), a2f_raw(nomega), a2f(nomega)
352 : :
353 : : !Local variables-------------------------------
354 : : !scalars
355 : : integer :: ncerr
356 : : !----------------------------------------------------------------------
357 : :
358 : 0 : ncerr = nctk_def_dims(ncid, nctkdim_t("a2f_nomega", nomega), defmode=.True.)
359 [ # # ]: 0 : NCF_CHECK(ncerr)
360 : :
361 : : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: &
362 [ # # ]: 0 : "isome_schema_version", "eph_intmeth", "ph_intmeth"])
363 [ # # ]: 0 : NCF_CHECK(ncerr)
364 : : ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: &
365 [ # # ]: 0 : "eph_fsmear", "ph_smear", "ph_wstep", "a2f_edos_fermie", "a2f_dos_normalization", "omega_2"])
366 [ # # ]: 0 : NCF_CHECK(ncerr)
367 : :
368 : : ncerr = nctk_def_arrays(ncid, [ &
369 : : nctkarr_t("gstore_ngqpt", "int", "three"), &
370 : : nctkarr_t("eph_ngqpt_fine", "int", "three"), &
371 : : nctkarr_t("ddb_ngqpt", "int", "three"), &
372 : : nctkarr_t("a2f_mesh", "dp", "a2f_nomega"), &
373 : : nctkarr_t("a2f_values_raw", "dp", "a2f_nomega"), &
374 [ # # ]: 0 : nctkarr_t("a2f_values", "dp", "a2f_nomega")])
375 [ # # ]: 0 : NCF_CHECK(ncerr)
376 : :
377 [ # # ]: 0 : NCF_CHECK(nctk_set_datamode(ncid))
378 : : ncerr = nctk_write_iscalars(ncid, &
379 : : [character(len=nctk_slen) :: "isome_schema_version", "eph_intmeth", "ph_intmeth"], &
380 [ # # # # ]: 0 : [3, dtset%eph_intmeth, dtset%ph_intmeth])
381 [ # # ]: 0 : NCF_CHECK(ncerr)
382 : : ncerr = nctk_write_dpscalars(ncid, &
383 : : [character(len=nctk_slen) :: &
384 : : "eph_fsmear", "ph_smear", "ph_wstep", "a2f_edos_fermie", "a2f_dos_normalization", "omega_2"], &
385 [ # # # # ]: 0 : [dtset%eph_fsmear, dtset%ph_smear, dtset%ph_wstep, edos_fermie, edos_fermie / two, omega_2])
386 [ # # ]: 0 : NCF_CHECK(ncerr)
387 : :
388 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gstore_ngqpt"), gstore%ngqpt))
389 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "eph_ngqpt_fine"), dtset%eph_ngqpt_fine))
390 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ddb_ngqpt"), dtset%ddb_ngqpt))
391 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "a2f_mesh"), omega))
392 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "a2f_values_raw"), a2f_raw))
393 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "a2f_values"), a2f))
394 : :
395 : 0 : end subroutine isome_ncwrite_spectral
396 : : !!***
397 : :
398 : : !----------------------------------------------------------------------
399 : :
400 : : !!****f* m_migdal_eliashberg/isome_ncwrite_qibz
401 : : !! NAME
402 : : !! isome_ncwrite_qibz
403 : : !!
404 : : !! FUNCTION
405 : : !! Write phonon frequencies and the normalized mode-resolved lambda(q,nu) on
406 : : !! the phonon IBZ. No averaging is applied inside degenerate subspaces;
407 : : !! their gauge-invariant contribution enters a2F through the mode sum.
408 : : !!
409 : : !! SOURCE
410 : :
411 [ # # # # ]: 0 : subroutine isome_ncwrite_qibz(ncid, gstore, phfreq, phlambda)
412 : :
413 : : !Arguments ------------------------------------
414 : : integer,intent(in) :: ncid
415 : : type(gstore_t),intent(in) :: gstore
416 : : real(dp),intent(in) :: phfreq(:,:), phlambda(:,:,:)
417 : :
418 : : !Local variables-------------------------------
419 : : integer :: ncerr, natom3
420 : : !----------------------------------------------------------------------
421 : :
422 : 0 : natom3 = size(phfreq, dim=1)
423 [ # # ]: 0 : ABI_CHECK_IEQ(size(phfreq, dim=2), gstore%nqibz, "Invalid phfreq q-IBZ dimension")
424 [ # # ]: 0 : ABI_CHECK_IEQ(size(phlambda, dim=1), natom3, "Invalid phlambda mode dimension")
425 [ # # ]: 0 : ABI_CHECK_IEQ(size(phlambda, dim=2), gstore%nqibz, "Invalid phlambda q-IBZ dimension")
426 [ # # ]: 0 : ABI_CHECK_IEQ(size(phlambda, dim=3), gstore%nsppol, "Invalid phlambda spin dimension")
427 : :
428 : : ncerr = nctk_def_dims(ncid, [nctkdim_t("isome_nqibz", gstore%nqibz), &
429 [ # # ]: 0 : nctkdim_t("isome_natom3", natom3)], defmode=.True.)
430 [ # # ]: 0 : NCF_CHECK(ncerr)
431 : : ncerr = nctk_def_arrays(ncid, [ &
432 : : nctkarr_t("qibz", "dp", "number_of_reduced_dimensions, isome_nqibz"), &
433 : : nctkarr_t("wtq", "dp", "isome_nqibz"), &
434 : : nctkarr_t("phfreq_qibz", "dp", "isome_natom3, isome_nqibz"), &
435 [ # # ]: 0 : nctkarr_t("phlambda_qibz", "dp", "isome_natom3, isome_nqibz, number_of_spins")])
436 [ # # ]: 0 : NCF_CHECK(ncerr)
437 : :
438 [ # # ]: 0 : NCF_CHECK(nctk_set_datamode(ncid))
439 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "qibz"), gstore%qibz))
440 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "wtq"), gstore%wtq))
441 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "phfreq_qibz"), phfreq))
442 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "phlambda_qibz"), phlambda))
443 : :
444 : 0 : end subroutine isome_ncwrite_qibz
445 : : !!***
446 : :
447 : : !----------------------------------------------------------------------
448 : :
449 : : !!****f* m_migdal_eliashberg/isome_ncwrite_qpath
450 : : !! NAME
451 : : !! isome_ncwrite_qpath
452 : : !!
453 : : !! FUNCTION
454 : : !! Write phonon frequencies, displacements, and mode-resolved
455 : : !! electron-phonon coupling lambda(q,nu) along the input phonon path.
456 : : !!
457 : : !! SOURCE
458 : :
459 [ # # # # : 0 : subroutine isome_ncwrite_qpath(ncid, qpath, phfreq, phdispl_cart, phlambda)
# # # # ]
460 : :
461 : : !Arguments ------------------------------------
462 : : !scalars
463 : : integer,intent(in) :: ncid
464 : : !arrays
465 : : real(dp),intent(in) :: qpath(:,:), phfreq(:,:), phdispl_cart(:,:,:,:), phlambda(:,:,:)
466 : :
467 : : !Local variables-------------------------------
468 : : !scalars
469 : : integer :: ncerr, nqpath, natom3, nsppol
470 : : !----------------------------------------------------------------------
471 : :
472 : 0 : nqpath = size(qpath, dim=2)
473 : 0 : natom3 = size(phfreq, dim=1)
474 : 0 : nsppol = size(phlambda, dim=3)
475 : :
476 : : ncerr = nctk_def_dims(ncid, [ &
477 : : nctkdim_t("isome_nqpath", nqpath), &
478 : : nctkdim_t("number_of_phonon_modes", natom3), &
479 [ # # ]: 0 : nctkdim_t("isome_nsppol", nsppol)], defmode=.True.)
480 [ # # ]: 0 : NCF_CHECK(ncerr)
481 : :
482 [ # # ]: 0 : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "phlambda_qpath_average_degenerate"])
483 [ # # ]: 0 : NCF_CHECK(ncerr)
484 [ # # ]: 0 : ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "phlambda_qpath_degen_tol"])
485 [ # # ]: 0 : NCF_CHECK(ncerr)
486 : :
487 : : ncerr = nctk_def_arrays(ncid, [ &
488 : : nctkarr_t("qpath", "dp", "number_of_reduced_dimensions, isome_nqpath"), &
489 : : nctkarr_t("phfreq_qpath", "dp", "number_of_phonon_modes, isome_nqpath"), &
490 : : nctkarr_t("phdispl_cart_qpath", "dp", "two, number_of_phonon_modes, number_of_phonon_modes, isome_nqpath"), &
491 [ # # ]: 0 : nctkarr_t("phlambda_qpath", "dp", "number_of_phonon_modes, isome_nqpath, isome_nsppol")])
492 [ # # ]: 0 : NCF_CHECK(ncerr)
493 : :
494 [ # # ]: 0 : NCF_CHECK(nctk_set_datamode(ncid))
495 : : ncerr = nctk_write_iscalars(ncid, &
496 : : [character(len=nctk_slen) :: "phlambda_qpath_average_degenerate"], &
497 [ # # ]: 0 : [merge(1, 0, average_lambda_qpath_degenerate)])
498 [ # # ]: 0 : NCF_CHECK(ncerr)
499 [ # # ]: 0 : ncerr = nctk_write_dpscalars(ncid, [character(len=nctk_slen) :: "phlambda_qpath_degen_tol"], [lambda_qpath_degen_tol])
500 [ # # ]: 0 : NCF_CHECK(ncerr)
501 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "qpath"), qpath))
502 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "phfreq_qpath"), phfreq))
503 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "phdispl_cart_qpath"), phdispl_cart))
504 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "phlambda_qpath"), phlambda))
505 : :
506 : 0 : end subroutine isome_ncwrite_qpath
507 : : !!***
508 : :
509 : : !----------------------------------------------------------------------
510 : :
511 : : !!****f* m_migdal_eliashberg/get_lambda_qpath_wan
512 : : !! NAME
513 : : !! get_lambda_qpath_wan
514 : : !!
515 : : !! FUNCTION
516 : : !! Compute lambda(q,nu) on the path defined by ph_qpath and ph_ndivsm.
517 : : !! The path construction is the same as in ifc_mkphbs. Arbitrary-q
518 : : !! electronic energies and matrix elements are obtained from the Wannier
519 : : !! Hamiltonian and GWAN interpolation, respectively.
520 : : !!
521 : : !! SOURCE
522 : :
523 : 0 : subroutine get_lambda_qpath_wan(gstore, dtset, edos_fermie, qpoints, phfreq, phdispl_cart, phlambda)
524 : :
525 : : !Arguments ------------------------------------
526 : : !scalars
527 : : class(gstore_t),intent(inout) :: gstore
528 : : type(dataset_type),intent(in) :: dtset
529 : : real(dp),intent(in) :: edos_fermie
530 : : !arrays
531 : : real(dp),allocatable,intent(out) :: qpoints(:,:), phfreq(:,:), phdispl_cart(:,:,:,:), phlambda(:,:,:)
532 : :
533 : : !Local variables-------------------------------
534 : : !scalars
535 : : integer,parameter :: master = 0
536 : : integer :: my_is, my_ik, my_ip, iq, ipc, nu, in_k, im_kq, ierr, my_rank
537 : : integer :: natom, natom3, nwan, nqpath, spin, ik_start, ik_stop, ikb, nkb, nk_batch
538 : : real(dp) :: weight_k, g2, fs_weight, spin_factor, cpu, wall, gflops
539 : : logical :: has_gwan
540 : 0 : type(kpath_t) :: qpath
541 : : character(len=500) :: msg
542 : : !arrays
543 : : integer :: units(2)
544 : 0 : real(dp),allocatable :: displ_cart4(:,:,:,:), displ_red4(:,:,:,:), displ_red(:,:,:)
545 : 0 : real(dp),allocatable :: gatm_real(:,:,:,:,:), gnu_real(:,:,:,:,:)
546 : 0 : real(dp),allocatable :: eig_k(:,:), eig_kq(:,:)
547 : 0 : complex(dp),allocatable :: intp_gatm(:,:,:,:), gatm_full(:,:,:,:), g_req(:,:,:,:)
548 : : !----------------------------------------------------------------------
549 : :
550 [ # # ]: 0 : units = [std_out, ab_out]
551 : :
552 [ # # # # ]: 0 : if (dtset%ph_nqpath <= 0 .or. dtset%ph_ndivsm <= 0) return
553 : :
554 : 0 : has_gwan = gstore%my_nspins > 0
555 [ # # ]: 0 : if (has_gwan) has_gwan = allocated(gstore%gqk(1)%wan%grpe_wwp)
556 [ # # ]: 0 : if (.not. has_gwan) then
557 [ # # ]: 0 : call wrtout([std_out, ab_out], " Skipping lambda(q,nu): GWAN interpolation is not available for this GSTORE.")
558 : 0 : return
559 : : end if
560 : :
561 [ # # ]: 0 : ABI_CHECK(dtset%eph_fsmear > zero, "lambda(q,nu) along a path requires a positive eph_fsmear for Gaussian Fermi-surface integration")
562 [ # # ]: 0 : ABI_CHECK(edos_fermie > zero, "The electronic DOS at the Fermi level must be positive")
563 [ # # ]: 0 : ABI_CHECK(gstore%kzone == "bz", "get_lambda_qpath_wan requires gstore_kzone = 'bz'")
564 : 0 : call cwtime(cpu, wall, gflops, "start")
565 : :
566 : 0 : natom = gstore%cryst%natom; natom3 = 3 * natom
567 : 0 : call qpath%init(dtset%ph_qpath(:,1:dtset%ph_nqpath), gstore%cryst%gprimd, dtset%ph_ndivsm)
568 : 0 : nqpath = qpath%npts
569 : :
570 [ # # # # : 0 : ABI_MALLOC(qpoints, (3, nqpath))
# # # # #
# # # #
# ]
571 [ # # # # : 0 : ABI_MALLOC(phfreq, (natom3, nqpath))
# # # # #
# # # # #
# # # # ]
572 [ # # # # : 0 : ABI_MALLOC(phdispl_cart, (2, natom3, natom3, nqpath))
# # # # #
# # # # #
# # # # #
# # # ]
573 [ # # # # : 0 : ABI_CALLOC(phlambda, (natom3, nqpath, gstore%nsppol))
# # # # #
# # # # #
# # # # #
# # # # #
# # # # ]
574 [ # # # # : 0 : ABI_MALLOC(displ_cart4, (2, 3, natom, natom3))
# # # # #
# # # # #
# # ]
575 [ # # # # : 0 : ABI_MALLOC(displ_red4, (2, 3, natom, natom3))
# # # # #
# # # ]
576 [ # # # # : 0 : ABI_MALLOC(displ_red, (2, natom3, natom3))
# # # # #
# # # # #
# # ]
577 : :
578 [ # # # # : 0 : qpoints = qpath%points
# # # # #
# # # # #
# # ]
579 [ # # ]: 0 : do iq=1,nqpath
580 : 0 : call gstore%ifc%fourq(gstore%cryst, qpoints(:,iq), phfreq(:,iq), displ_cart4, out_displ_red=displ_red4)
581 [ # # ]: 0 : phdispl_cart(:,:,:,iq) = reshape(displ_cart4, [2, natom3, natom3])
582 : : end do
583 : :
584 : 0 : write(msg, "(a,i0,a,es12.4,a)")" Computing lambda(q,nu) at ", nqpath, &
585 : 0 : " path points with eph_fsmear: ", dtset%eph_fsmear * Ha_meV, " meV"
586 : 0 : call wrtout(units, msg, pre_newlines=1)
587 : :
588 [ # # ]: 0 : do my_is=1,gstore%my_nspins
589 : 0 : associate (gqk => gstore%gqk(my_is))
590 : 0 : spin = gqk%spin; nwan = gqk%wan%nwan
591 : 0 : ABI_CHECK_IEQ(nwan, gqk%nb_k, "Wannier and GSTORE band dimensions differ")
592 [ # # ]: 0 : ABI_CHECK_IEQ(nwan, gqk%nb_kq, "Wannier and GSTORE band dimensions differ")
593 : :
594 : 0 : nk_batch = gqk%wan%eph_kbatch_size(gqk%my_nk, natom3)
595 [ # # # # : 0 : ABI_MALLOC(intp_gatm, (nwan, nwan, gqk%my_npert, nk_batch))
# # # # #
# # # # #
# # # # #
# # # #
# ]
596 [ # # # # : 0 : ABI_MALLOC(gatm_full, (nwan, nwan, natom3, nk_batch))
# # # # #
# # # # #
# # # # #
# # # #
# ]
597 [ # # # # : 0 : ABI_MALLOC(g_req, (gqk%wan%nr_e, nwan, nwan, gqk%my_npert))
# # # # #
# # # # #
# # # # #
# # # #
# ]
598 [ # # # # : 0 : ABI_MALLOC(gatm_real, (2, nwan, nwan, 1, natom3))
# # # # #
# # # # #
# # # # #
# ]
599 [ # # # # : 0 : ABI_MALLOC(gnu_real, (2, nwan, nwan, 1, natom3))
# # # # #
# # # # #
# # ]
600 [ # # # # : 0 : ABI_MALLOC(eig_k, (nwan, nk_batch))
# # # # #
# # # # #
# # ]
601 [ # # # # : 0 : ABI_MALLOC(eig_kq, (nwan, nk_batch))
# # # # #
# # # ]
602 : :
603 [ # # ]: 0 : do iq=1,nqpath
604 [ # # ]: 0 : if (gqk%qpt_comm%skip(iq)) cycle
605 : 0 : call gstore%ifc%fourq(gstore%cryst, qpoints(:,iq), phfreq(:,iq), displ_cart4, out_displ_red=displ_red4)
606 [ # # ]: 0 : displ_red = reshape(displ_red4, [2, natom3, natom3])
607 : 0 : call gqk%wan%prepare_eph_q(qpoints(:,iq), g_req)
608 : :
609 [ # # # # : 0 : do ik_start=1,gqk%my_nk,nk_batch
# # # # ]
610 : 0 : ik_stop = min(gqk%my_nk, ik_start + nk_batch - 1); nkb = ik_stop - ik_start + 1
611 : :
612 : : ! Interpolate one bounded k block and reuse the returned electronic
613 : : ! energies in the Fermi-surface weights.
614 : : call gqk%wan%interp_eph_manyk_from_q(gstore%cryst, nkb, gqk%my_kpts(:,ik_start:ik_stop), &
615 : : qpoints(:,iq), g_req, intp_gatm(:,:,:,1:nkb), &
616 : 0 : out_eigens_k=eig_k(:,1:nkb), out_eigens_kq=eig_kq(:,1:nkb))
617 : :
618 [ # # # # : 0 : gatm_full(:,:,:,1:nkb) = czero
# # # # ]
619 [ # # ]: 0 : do ikb=1,nkb
620 [ # # ]: 0 : do ipc=1,gqk%my_npert
621 [ # # # # ]: 0 : gatm_full(:,:,gqk%my_pertcases(ipc),ikb) = intp_gatm(:,:,ipc,ikb)
622 : : end do
623 : : end do
624 [ # # ]: 0 : if (gqk%pert_comm%nproc > 1) call xmpi_sum(gatm_full(:,:,:,1:nkb), gqk%pert_comm%value, ierr)
625 : :
626 [ # # ]: 0 : do ikb=1,nkb
627 : 0 : my_ik = ik_start + ikb - 1
628 : 0 : weight_k = gqk%my_wtk(my_ik)
629 [ # # # # : 0 : gatm_real(1,:,:,1,:) = real(gatm_full(:,:,:,ikb), kind=dp)
# # ]
630 [ # # # # : 0 : gatm_real(2,:,:,1,:) = aimag(gatm_full(:,:,:,ikb))
# # ]
631 : 0 : call ephtk_gkknu_from_atm(nwan, nwan, 1, natom, gatm_real, phfreq(:,iq), displ_red, gnu_real)
632 : :
633 [ # # ]: 0 : do my_ip=1,gqk%my_npert
634 : 0 : nu = gqk%my_pertcases(my_ip)
635 [ # # ]: 0 : if (phfreq(nu,iq) < EPHTK_WTOL) cycle
636 [ # # ]: 0 : do in_k=1,nwan
637 [ # # ]: 0 : do im_kq=1,nwan
638 : 0 : g2 = gnu_real(1,im_kq,in_k,1,nu)**2 + gnu_real(2,im_kq,in_k,1,nu)**2
639 : : fs_weight = gaussian(eig_k(in_k,ikb) - gstore%ebands%fermie, dtset%eph_fsmear) * &
640 : 0 : gaussian(eig_kq(im_kq,ikb) - gstore%ebands%fermie, dtset%eph_fsmear)
641 : 0 : phlambda(nu,iq,spin) = phlambda(nu,iq,spin) + two * g2 * weight_k * fs_weight / phfreq(nu,iq)
642 : : end do
643 : : end do
644 : : end do
645 : : end do ! ikb
646 : : end do ! ik_start
647 : : end do
648 : :
649 [ # # ]: 0 : ABI_FREE(intp_gatm)
650 [ # # ]: 0 : ABI_FREE(gatm_full)
651 [ # # ]: 0 : ABI_FREE(g_req)
652 [ # # ]: 0 : ABI_FREE(gatm_real)
653 [ # # ]: 0 : ABI_FREE(gnu_real)
654 [ # # ]: 0 : ABI_FREE(eig_k)
655 [ # # # # ]: 0 : ABI_FREE(eig_kq)
656 : : end associate
657 : : end do
658 : :
659 : 0 : spin_factor = two / (gstore%nsppol * dtset%nspinor)
660 [ # # # # : 0 : phlambda = phlambda * spin_factor / edos_fermie
# # ]
661 : 0 : call xmpi_sum(phlambda, gstore%comm, ierr)
662 : 0 : if (average_lambda_qpath_degenerate) call average_lambda_degenerate_modes(nqpath, natom3, gstore%nsppol, phfreq, phlambda)
663 : :
664 : 0 : my_rank = xmpi_comm_rank(gstore%comm)
665 [ # # ]: 0 : if (my_rank == master) then
666 : 0 : call wrtout(ab_out, " Phonon frequencies and lambda(q,nu) along the q-path:", pre_newlines=1)
667 [ # # ]: 0 : do iq=1,nqpath
668 : 0 : write(msg, "(a,i0,a,3es16.8)")" q-path point ", iq, ": ", qpoints(:,iq)
669 : 0 : call wrtout(ab_out, msg)
670 : 0 : select case (gstore%nsppol)
671 : : case (1)
672 : 0 : call wrtout(ab_out, " nu omega (meV) lambda")
673 [ # # ]: 0 : do nu=1,natom3
674 : 0 : write(msg, "(i5,2x,es16.8,2x,es16.8)")nu, phfreq(nu,iq) * Ha_meV, phlambda(nu,iq,1)
675 : 0 : call wrtout(ab_out, msg)
676 : : end do
677 : : case (2)
678 : 0 : call wrtout(ab_out, " nu omega (meV) lambda_spin1 lambda_spin2")
679 [ # # ]: 0 : do nu=1,natom3
680 : 0 : write(msg, "(i5,2x,es16.8,2x,es16.8,2x,es16.8)")nu, phfreq(nu,iq) * Ha_meV, phlambda(nu,iq,1:2)
681 : 0 : call wrtout(ab_out, msg)
682 : : end do
683 : : case default
684 [ # # # ]: 0 : ABI_ERROR("Printing lambda(q,nu) supports only nsppol = 1 or 2")
685 : : end select
686 : : end do
687 : : end if
688 : :
689 [ # # ]: 0 : ABI_FREE(displ_cart4)
690 [ # # ]: 0 : ABI_FREE(displ_red4)
691 [ # # ]: 0 : ABI_FREE(displ_red)
692 : 0 : call qpath%free()
693 : 0 : call cwtime_report(" Wannier interpolation of lambda(q,nu) along q-path", cpu, wall, gflops)
694 : :
695 [ # # # # : 0 : end subroutine get_lambda_qpath_wan
# # # # #
# # # #
# ]
696 : : !!***
697 : :
698 : : !----------------------------------------------------------------------
699 : :
700 : : !!****f* m_migdal_eliashberg/average_lambda_degenerate_modes
701 : : !! NAME
702 : : !! average_lambda_degenerate_modes
703 : : !!
704 : : !! FUNCTION
705 : : !! Replace lambda(q,nu) inside each degenerate phonon subspace with its
706 : : !! arithmetic average. Frequencies are assumed to be ordered by branch,
707 : : !! as returned by ifc%fourq. The sum over every subspace is preserved.
708 : : !!
709 : : !! SOURCE
710 : :
711 : 0 : subroutine average_lambda_degenerate_modes(nqpath, nmode, nsppol, phfreq, phlambda)
712 : :
713 : : !Arguments ------------------------------------
714 : : !scalars
715 : : integer,intent(in) :: nqpath, nmode, nsppol
716 : : !arrays
717 : : real(dp),intent(in) :: phfreq(nmode,nqpath)
718 : : real(dp),intent(inout) :: phlambda(nmode,nqpath,nsppol)
719 : :
720 : : !Local variables-------------------------------
721 : : !scalars
722 : : integer :: iq, spin, first_mode, last_mode, ndeg
723 : : real(dp) :: lambda_avg
724 : : !----------------------------------------------------------------------
725 : :
726 [ # # ]: 0 : do iq=1,nqpath
727 : : first_mode = 1
728 [ # # ]: 0 : do while (first_mode <= nmode)
729 : : last_mode = first_mode
730 [ # # ]: 0 : do while (last_mode < nmode)
731 [ # # ]: 0 : if (abs(phfreq(last_mode + 1,iq) - phfreq(first_mode,iq)) > lambda_qpath_degen_tol) exit
732 : 0 : last_mode = last_mode + 1
733 : : end do
734 : :
735 : 0 : ndeg = last_mode - first_mode + 1
736 [ # # ]: 0 : if (ndeg > 1) then
737 [ # # ]: 0 : do spin=1,nsppol
738 [ # # ]: 0 : lambda_avg = sum(phlambda(first_mode:last_mode,iq,spin)) / ndeg
739 [ # # ]: 0 : phlambda(first_mode:last_mode,iq,spin) = lambda_avg
740 : : end do
741 : : end if
742 : 0 : first_mode = last_mode + 1
743 : : end do
744 : : end do
745 : :
746 : 0 : end subroutine average_lambda_degenerate_modes
747 : : !!***
748 : :
749 : : !----------------------------------------------------------------------
750 : :
751 : : !!****f* m_migdal_eliashberg/get_a2fw
752 : : !! NAME
753 : : !! get_a2fw
754 : : !!
755 : : !! FUNCTION
756 : : !! Compute the raw, mode-resolved lambda(q,nu) on the phonon IBZ and
757 : : !! construct a^2F(omega) from its weighted sum. Phonon modes are not
758 : : !! averaged when they are degenerate.
759 : : !!
760 : : !! INPUTS
761 : : !! gstore: Electron-phonon matrix elements. Electronic k-points must cover
762 : : !! the full BZ; phonon q-points may cover either the BZ or the IBZ.
763 : : !! edos_fermie: Electronic DOS at the Fermi level.
764 : : !! nw: Number of frequencies.
765 : : !! wmesh: Frequency mesh.
766 : : !!
767 : : !! OUTPUT
768 : : !! a2fw(nw): Eliashberg function in the historical raw normalization.
769 : : !! phfreq_qibz: Phonon frequencies on the q-point IBZ.
770 : : !! phlambda_qibz: DOS-normalized mode- and spin-resolved lambda(q,nu).
771 : : !!
772 : : !! SOURCE
773 : :
774 : 0 : subroutine get_a2fw(gstore, edos_fermie, nw, wmesh, a2fw, phfreq_qibz, phlambda_qibz)
775 : :
776 : : !Arguments ------------------------------------
777 : : class(gstore_t),intent(inout) :: gstore
778 : : integer,intent(in) :: nw
779 : : real(dp),intent(in) :: edos_fermie
780 : : real(dp),intent(in) :: wmesh(nw)
781 : : real(dp),intent(out) :: a2fw(nw)
782 : : real(dp),allocatable,intent(out) :: phfreq_qibz(:,:), phlambda_qibz(:,:,:)
783 : :
784 : : !Local variables-------------------------------
785 : : integer :: my_is, my_ik, my_iq, my_ip, in_k, im_kq, ierr, iq_ibz, isym_q, trev_q, nb_k, nb_kq
786 : : integer :: natom, natom3, g0_q(3), nk_batch, ik_start, ik_stop, ikb, nkb, ipc, nwan
787 : : real(dp) :: g2_qnu, wqnu, weight_k, weight_q, cpu, wall, gflops, spin_factor
788 : : logical :: isirr_q
789 : : !arrays
790 : : integer :: units(2)
791 : : real(dp) :: qpt(3)
792 : 0 : real(dp),allocatable :: dbl_delta_q(:,:,:), g2_mnkp(:,:,:,:), deltaw_nuq(:), displ_cart_dum(:,:,:,:)
793 : 0 : real(dp),allocatable :: displ_red_local(:,:,:)
794 : 0 : complex(dp),allocatable :: intp_gatm(:,:,:,:), gatm_full(:,:,:,:), g_req(:,:,:,:), gnu(:,:)
795 : : !----------------------------------------------------------------------
796 : :
797 [ # # ]: 0 : units = [std_out, ab_out]
798 : :
799 : 0 : call cwtime(cpu, wall, gflops, "start")
800 : 0 : call wrtout(units, sjoin(" Computing a^2F(w) with ph_smear:", ftoa(gstore%dtset%ph_smear * Ha_meV), "(meV)"), pre_newlines=1)
801 : :
802 [ # # ]: 0 : ABI_CHECK(gstore%kzone == "bz", "get_a2fw requires kzone == `bz`")
803 [ # # # # : 0 : ABI_CHECK(any(gstore%qzone == ["bz ", "ibz"]), "get_a2fw requires qzone == `bz` or `ibz`")
# # ]
804 [ # # ]: 0 : ABI_CHECK(any(gstore%with_cplex == [0, 1]), "get_a2fw requires with_cplex=0 or 1")
805 [ # # ]: 0 : if (gstore%with_cplex == 0) then
806 [ # # ]: 0 : ABI_CHECK(gstore%has_wannier, "get_a2fw with_cplex=0 requires a Wannier interpolator")
807 : : end if
808 [ # # ]: 0 : ABI_CHECK(edos_fermie > zero, "The electronic DOS at the Fermi level must be positive")
809 : :
810 [ # # ]: 0 : ABI_CHECK(allocated(gstore%gqk), "gstore%gqk is not allocated")
811 [ # # ]: 0 : ABI_CHECK(size(gstore%gqk) == gstore%my_nspins, "gstore%gqk is inconsistent with my_nspins")
812 : :
813 [ # # # # : 0 : ABI_MALLOC(deltaw_nuq, (nw))
# # # # #
# # # ]
814 : 0 : natom = gstore%cryst%natom
815 : 0 : natom3 = 3 * natom
816 [ # # # # : 0 : ABI_MALLOC(phfreq_qibz, (natom3, gstore%nqibz))
# # # # #
# # # # #
# # # # ]
817 [ # # # # : 0 : ABI_CALLOC(phlambda_qibz, (natom3, gstore%nqibz, gstore%nsppol))
# # # # #
# # # # #
# # # # #
# # # # #
# # # # ]
818 [ # # # # : 0 : ABI_MALLOC(displ_cart_dum, (2, 3, gstore%cryst%natom, natom3))
# # # # #
# # # # #
# # ]
819 [ # # ]: 0 : do iq_ibz=1,gstore%nqibz
820 : 0 : call gstore%ifc%fourq(gstore%cryst, gstore%qibz(:,iq_ibz), phfreq_qibz(:,iq_ibz), displ_cart_dum)
821 : : end do
822 [ # # ]: 0 : ABI_FREE(displ_cart_dum)
823 [ # # ]: 0 : a2fw = zero
824 : :
825 : : ! Loop over collinear spins.
826 [ # # ]: 0 : do my_is=1,gstore%my_nspins
827 : 0 : associate (gqk => gstore%gqk(my_is), cryst => gstore%cryst)
828 : 0 : if (gstore%with_cplex == 1) then
829 [ # # ]: 0 : ABI_CHECK(allocated(gqk%my_g2), "my_g2 is not allocated")
830 : : end if
831 [ # # ]: 0 : ABI_CHECK(allocated(gqk%my_wnuq), "my_wnuq is not allocated")
832 : :
833 : 0 : nb_k = gqk%nb_k; nb_kq = gqk%nb_kq
834 [ # # ]: 0 : ABI_CHECK_IEQ(nb_k, nb_kq, "gqk_dbldelta_qpt does not support nb_k != nb_kq")
835 [ # # ]: 0 : do ipc=1,gqk%my_npert
836 [ # # ]: 0 : ABI_CHECK(gqk%my_pertcases(ipc) == gqk%my_pert_start + ipc - 1, "gqk%my_pertcases is inconsistent with my_pert_start")
837 : : end do
838 : :
839 : : ! Weights for delta(e_{m k+q}) delta(e_{n k}) for my list of k-points.
840 [ # # # # : 0 : ABI_MALLOC(dbl_delta_q, (nb_kq, nb_k, gqk%my_nk))
# # # # #
# # # # #
# # # # #
# ]
841 [ # # # # : 0 : ABI_MALLOC(g2_mnkp, (nb_kq, nb_k, gqk%my_nk, gqk%my_npert))
# # # # #
# # # # #
# # # # #
# # # #
# ]
842 : :
843 [ # # ]: 0 : if (gstore%with_cplex == 0) then
844 : 0 : nwan = gqk%wan%nwan
845 : :
846 : 0 : nk_batch = gqk%wan%eph_kbatch_size(gqk%my_nk, natom3)
847 [ # # # # ]: 0 : ABI_CHECK(nk_batch >= 1 .and. nk_batch <= gqk%my_nk, "Invalid Wannier e-ph k-batch size")
848 : :
849 [ # # # # : 0 : ABI_MALLOC(intp_gatm, (nwan, nwan, gqk%my_npert, nk_batch))
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
850 [ # # # # : 0 : ABI_MALLOC(gatm_full, (nwan, nwan, natom3, nk_batch))
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
851 [ # # # # : 0 : ABI_MALLOC(g_req, (gqk%wan%nr_e, nwan, nwan, gqk%my_npert))
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
852 [ # # # # : 0 : ABI_MALLOC(gnu, (nwan, nwan))
# # # # #
# # # # #
# # # # ]
853 [ # # # # : 0 : ABI_MALLOC(displ_red_local, (2, natom3, gqk%my_npert))
# # # # #
# # # # #
# # # # ]
854 : : end if
855 : :
856 : : ! Loop over my q-points.
857 [ # # ]: 0 : do my_iq=1,gqk%my_nq
858 : 0 : iq_ibz = gqk%my_q2ibz(1,my_iq)
859 : 0 : isym_q = gqk%my_q2ibz(2,my_iq)
860 : 0 : trev_q = gqk%my_q2ibz(6,my_iq)
861 [ # # ]: 0 : g0_q = gqk%my_q2ibz(3:5,my_iq)
862 [ # # # # : 0 : isirr_q = isym_q == 1 .and. trev_q == 0 .and. all(g0_q == 0)
# # # # ]
863 [ # # # # ]: 0 : if (gstore%qzone == "bz" .and. .not. isirr_q) cycle
864 : :
865 : : ! Compute all integration weights for the double delta.
866 : 0 : call gqk%dbldelta_qpt(my_iq, gstore, gstore%dtset%eph_intmeth, gstore%dtset%eph_fsmear, qpt, weight_q, dbl_delta_q)
867 : :
868 [ # # ]: 0 : if (gstore%with_cplex == 1) then
869 : : ! Copy data to improve memory access in the loops below.
870 [ # # ]: 0 : do my_ip=1,gqk%my_npert
871 [ # # # # : 0 : g2_mnkp(:,:,:,my_ip) = gqk%my_g2(my_ip,:,my_iq,:,:)
# # ]
872 : : end do
873 : : else
874 : : ! Convert the symmetry-consistent phonon eigenvectors already stored
875 : : ! in gstore to reduced coordinates for the locally owned modes.
876 [ # # ]: 0 : do my_ip=1,gqk%my_npert
877 : : call phdispl_cart2red_nmodes(natom, 1, cryst%gprimd, &
878 : : gqk%my_displ_cart(:,:,:,my_ip,my_iq), &
879 : 0 : displ_red_local(:,:,my_ip:my_ip))
880 : : end do
881 : :
882 : : ! Interpolate bounded k blocks and immediately form |g_mnnu(k,q)|^2.
883 : 0 : call gqk%wan%prepare_eph_q(qpt, g_req)
884 [ # # # # : 0 : do ik_start=1,gqk%my_nk,nk_batch
# # # # ]
885 : 0 : ik_stop = min(gqk%my_nk, ik_start + nk_batch - 1)
886 : 0 : nkb = ik_stop - ik_start + 1
887 : : call gqk%wan%interp_eph_manyk_from_q(cryst, nkb, gqk%my_kpts(:,ik_start:ik_stop), qpt, &
888 : 0 : g_req, intp_gatm(:,:,:,1:nkb))
889 : :
890 [ # # # # : 0 : gatm_full(:,:,:,1:nkb) = czero
# # # # ]
891 [ # # ]: 0 : do ikb=1,nkb
892 [ # # ]: 0 : do ipc=1,gqk%my_npert
893 [ # # # # ]: 0 : gatm_full(:,:,gqk%my_pertcases(ipc),ikb) = intp_gatm(:,:,ipc,ikb)
894 : : end do
895 : : end do
896 [ # # ]: 0 : if (gqk%pert_comm%nproc > 1) call xmpi_sum(gatm_full(:,:,:,1:nkb), gqk%pert_comm%value, ierr)
897 : :
898 [ # # ]: 0 : do ikb=1,nkb
899 : 0 : my_ik = ik_start + ikb - 1
900 [ # # ]: 0 : do my_ip=1,gqk%my_npert
901 : 0 : wqnu = gqk%my_wnuq(my_ip,my_iq)
902 [ # # ]: 0 : if (wqnu < EPHTK_WTOL) then
903 [ # # # # ]: 0 : g2_mnkp(:,:,my_ik,my_ip) = zero; cycle
904 : : end if
905 [ # # # # ]: 0 : gnu = czero
906 [ # # ]: 0 : do ipc=1,natom3
907 : : gnu = gnu + gatm_full(:,:,ipc,ikb) * &
908 [ # # # # : 0 : (displ_red_local(1,ipc,my_ip) + j_dpc * displ_red_local(2,ipc,my_ip))
# # ]
909 : : end do
910 [ # # # # ]: 0 : gnu = gnu / sqrt(two * wqnu)
911 [ # # # # ]: 0 : g2_mnkp(:,:,my_ik,my_ip) = real(gnu * conjg(gnu), kind=dp)
912 : : end do
913 : : end do
914 : : end do
915 : : end if
916 : :
917 : : ! Loop over my phonon modes.
918 [ # # ]: 0 : do my_ip=1,gqk%my_npert
919 : 0 : wqnu = gqk%my_wnuq(my_ip, my_iq)
920 [ # # ]: 0 : if (wqnu < EPHTK_WTOL) cycle
921 : : ! Loop over my k-points.
922 [ # # ]: 0 : do my_ik=1,gqk%my_nk
923 : 0 : weight_k = gqk%my_wtk(my_ik)
924 : :
925 : : ! Sum over m_kq and n_k and accumulate.
926 [ # # ]: 0 : do in_k=1,nb_k
927 [ # # ]: 0 : do im_kq=1,nb_kq
928 : 0 : g2_qnu = g2_mnkp(im_kq, in_k, my_ik, my_ip)
929 : : phlambda_qibz(gqk%my_pertcases(my_ip),iq_ibz,gqk%spin) = &
930 : : phlambda_qibz(gqk%my_pertcases(my_ip),iq_ibz,gqk%spin) + two * g2_qnu * weight_k * &
931 : 0 : dbl_delta_q(im_kq, in_k, my_ik) / wqnu
932 : : end do
933 : : end do
934 : : end do
935 : : end do
936 : : end do ! my_iq
937 : :
938 [ # # ]: 0 : ABI_FREE(dbl_delta_q)
939 [ # # ]: 0 : ABI_FREE(g2_mnkp)
940 [ # # # # ]: 0 : if (gstore%with_cplex == 0) then
941 [ # # ]: 0 : ABI_FREE(intp_gatm)
942 [ # # ]: 0 : ABI_FREE(gatm_full)
943 [ # # ]: 0 : ABI_FREE(g_req)
944 [ # # ]: 0 : ABI_FREE(gnu)
945 [ # # ]: 0 : ABI_FREE(displ_red_local)
946 : : end if
947 : : end associate
948 : : end do ! my_is
949 : :
950 [ # # ]: 0 : ABI_FREE(deltaw_nuq)
951 : :
952 : : ! Normalize lambda(q,nu), then construct the normalized Eliashberg function
953 : : ! with the phonon-IBZ weights. edos_fermie is the total DOS, including spin
954 : : ! degeneracy. The factor two already present in the mode-resolved accumulator
955 : : ! is the conventional prefactor in lambda(q,nu), so do not divide the DOS by
956 : : ! two here.
957 : 0 : spin_factor = two / (gstore%nsppol * gstore%dtset%nspinor)
958 [ # # # # : 0 : phlambda_qibz = phlambda_qibz * spin_factor / edos_fermie
# # ]
959 : 0 : call xmpi_sum(phlambda_qibz, gstore%comm, ierr)
960 [ # # ]: 0 : a2fw = zero
961 [ # # ]: 0 : do iq_ibz=1,gstore%nqibz
962 [ # # ]: 0 : do my_ip=1,natom3
963 : 0 : wqnu = phfreq_qibz(my_ip,iq_ibz)
964 [ # # ]: 0 : if (wqnu < EPHTK_WTOL) cycle
965 [ # # # # ]: 0 : deltaw_nuq = gaussian(wmesh - wqnu, gstore%dtset%ph_smear)
966 [ # # # # ]: 0 : a2fw = a2fw + half * gstore%wtq(iq_ibz) * wqnu * sum(phlambda_qibz(my_ip,iq_ibz,:)) * deltaw_nuq
967 : : end do
968 : : end do
969 [ # # ]: 0 : a2fw = a2fw * (edos_fermie / two)
970 : :
971 : 0 : call cwtime_report(" get_a2fw", cpu, wall, gflops)
972 : :
973 [ # # # # : 0 : end subroutine get_a2fw
# # # # #
# # # ]
974 : : !!***
975 : :
976 : : !----------------------------------------------------------------------
977 : :
978 : : !!****f* m_migdal_eliashberg/iso_solver_free
979 : : !! NAME
980 : : !! iso_solver_free
981 : : !!
982 : : !! FUNCTION
983 : : !! Free dynamic memory
984 : : !!
985 : : !! SOURCE
986 : :
987 : 0 : subroutine iso_solver_free(iso)
988 : :
989 : : !Arguments ------------------------------------
990 : : class(iso_solver_t),intent(inout) :: iso
991 : : !----------------------------------------------------------------------
992 : :
993 [ # # ]: 0 : ABI_SFREE(iso%delta_iw)
994 [ # # ]: 0 : ABI_SFREE(iso%zeta_iw)
995 [ # # ]: 0 : ABI_SFREE(iso%prev_delta_iw)
996 [ # # ]: 0 : ABI_SFREE(iso%prev_zeta_iw)
997 [ # # ]: 0 : ABI_SFREE(iso%delta_iw_mix)
998 : :
999 : 0 : end subroutine iso_solver_free
1000 : : !!***
1001 : :
1002 : : !----------------------------------------------------------------------
1003 : :
1004 : : !!****f* m_migdal_eliashberg/iso_solver_solve
1005 : : !! NAME
1006 : : !! iso_solver_solve
1007 : : !!
1008 : : !! FUNCTION
1009 : : !!
1010 : : !! INPUTS
1011 : : !!
1012 : : !! OUTPUT
1013 : : !!
1014 : : !! SOURCE
1015 : :
1016 : 0 : subroutine iso_solver_solve(iso, itemp, kt, niw, imag_w, lambda_ij)
1017 : :
1018 : : !Arguments ------------------------------------
1019 : : !scalars
1020 : : class(iso_solver_t),intent(inout) :: iso
1021 : : integer,intent(in) :: itemp, niw
1022 : : real(dp),intent(in) :: kt
1023 : : !arrays
1024 : : real(dp),intent(in) :: imag_w(niw), lambda_ij(2 * niw)
1025 : :
1026 : : !Local variables-------------------------------
1027 : : !scalars
1028 : : integer,parameter :: master = 0
1029 : : integer :: nproc, my_rank, iter, ii, jj, converged
1030 : : real(dp) :: rr
1031 : : !arrays
1032 : 0 : real(dp),allocatable :: prev_vals(:)
1033 : :
1034 : : !----------------------------------------------------------------------
1035 : :
1036 : : ABI_UNUSED(lambda_ij)
1037 : :
1038 : 0 : nproc = xmpi_comm_size(iso%comm); my_rank = xmpi_comm_rank(iso%comm)
1039 : :
1040 [ # # # # : 0 : ABI_REMALLOC(iso%delta_iw_mix, (niw, iso%max_nmix))
# # # # #
# # # # #
# # # # #
# ]
1041 : :
1042 [ # # ]: 0 : if (itemp == 1) then
1043 : : ! Init values from scratch
1044 [ # # # # : 0 : ABI_CALLOC(iso%zeta_iw, (niw))
# # # # #
# # # # #
# # ]
1045 [ # # # # : 0 : ABI_CALLOC(iso%prev_zeta_iw, (niw))
# # # # #
# # # ]
1046 [ # # # # : 0 : ABI_CALLOC(iso%delta_iw, (niw))
# # # # #
# # # ]
1047 [ # # # # : 0 : ABI_CALLOC(iso%prev_delta_iw, (niw))
# # # # #
# # # ]
1048 : : else
1049 : : ! Init values from previous temperature. TODO: May use spline
1050 : 0 : call alloc_copy(iso%zeta_iw, prev_vals)
1051 [ # # # # : 0 : ABI_RECALLOC(iso%zeta_iw, (niw))
# # # # #
# # # # #
# # # # ]
1052 [ # # ]: 0 : ABI_MOVE_ALLOC(prev_vals, iso%prev_zeta_iw)
1053 : 0 : call alloc_copy(iso%delta_iw, prev_vals)
1054 [ # # # # : 0 : ABI_RECALLOC(iso%delta_iw, (niw))
# # # # #
# # # #
# ]
1055 [ # # ]: 0 : ABI_MOVE_ALLOC(prev_vals, iso%prev_delta_iw)
1056 : : end if
1057 : :
1058 : 0 : converged = 0
1059 [ # # ]: 0 : iter_loop: do iter=1,iso%max_niter
1060 : :
1061 [ # # ]: 0 : do ii=1,niw
1062 : : !if (mod(ii, nproc) /= my_rank) cycle ! MPI parallelism inside comm
1063 [ # # ]: 0 : do jj=1,niw
1064 : 0 : rr = one / sqrt(imag_w(jj) ** 2 + iso%prev_delta_iw(jj) ** 2)
1065 : 0 : iso%zeta_iw(ii) = iso%zeta_iw(ii) + imag_w(jj) * rr !* lambda(ii - jj)
1066 : 0 : iso%delta_iw(ii) = iso%delta_iw(ii) + rr * iso%prev_delta_iw(jj) !* (lambda(ii - jj) - mustar)
1067 : : end do
1068 : 0 : iso%zeta_iw(ii) = one + pi * kt / imag_w(ii) * iso%zeta_iw(ii)
1069 : 0 : iso%delta_iw(ii) = pi * kt * iso%delta_iw(ii) / iso%zeta_iw(ii)
1070 : : end do ! ii
1071 : :
1072 : : if (my_rank == master) then
1073 : : ! Write SCF cycle to stdout.
1074 : : ! Check for convergence.
1075 : : converged = 0
1076 : : end if
1077 : :
1078 : : if (converged == 2) exit iter_loop
1079 : :
1080 : : ! TODO: Mixing
1081 [ # # # # : 0 : iso%prev_zeta_iw = iso%zeta_iw
# # # # #
# ]
1082 [ # # # # : 0 : iso%prev_delta_iw = iso%delta_iw
# # # # #
# ]
1083 : :
1084 : : end do iter_loop
1085 : :
1086 : : ! Pade' to go to real axis
1087 : : ! Compute Delta F
1088 : : ! Compute QP DOS
1089 : :
1090 : : ! Write results to netcdf file
1091 : : if (my_rank == master) then
1092 : : end if
1093 : :
1094 [ # # ]: 0 : end subroutine iso_solver_solve
1095 : : !!***
1096 : :
1097 : : !----------------------------------------------------------------------
1098 : :
1099 : : !!****f* m_migdal_eliashberg/matsubara_mesh
1100 : : !! NAME
1101 : : !! mastubara_mesh
1102 : : !!
1103 : : !! FUNCTION
1104 : : !!
1105 : : !! INPUTS
1106 : : !!
1107 : : !! OUTPUT
1108 : : !!
1109 : : !! SOURCE
1110 : :
1111 : 0 : subroutine matsubara_mesh(bosons_or_fermions, kt, wmax, niw, imag_w)
1112 : :
1113 : : !Arguments ------------------------------------
1114 : : !scalars
1115 : : character(len=*),intent(in) :: bosons_or_fermions
1116 : : real(dp),intent(in) :: kt, wmax
1117 : : integer,intent(out) :: niw
1118 : : real(dp),allocatable,intent(out) :: imag_w(:)
1119 : :
1120 : : !Local variables-------------------------------
1121 : : integer :: nn
1122 : : !----------------------------------------------------------------------
1123 : :
1124 : 0 : select case (bosons_or_fermions)
1125 : : case ("bosons")
1126 : : ! 2 n pi kT
1127 : 0 : niw = nint(wmax / (two * pi * kt)) + 1
1128 [ # # # # : 0 : ABI_MALLOC(imag_w, (niw))
# # # # #
# # # #
# ]
1129 [ # # ]: 0 : do nn=0,niw-1
1130 : 0 : imag_w(nn + 1) = two * nn * pi * kt
1131 : : end do
1132 : :
1133 : : case ("fermions")
1134 : : ! 2 (n + 1) pi kT
1135 : 0 : niw = nint(wmax / (two * pi * kt))
1136 [ # # # # : 0 : ABI_MALLOC(imag_w, (niw))
# # # # #
# # # #
# ]
1137 [ # # ]: 0 : do nn=0,niw-1
1138 : 0 : imag_w(nn + 1) = two * (nn + 1) * pi * kt
1139 : : end do
1140 : :
1141 : : case default
1142 [ # # # ]: 0 : ABI_ERROR(sjoin("Wrong values for bosons_or_fermions:", bosons_or_fermions))
1143 : : end select
1144 : :
1145 : 0 : end subroutine matsubara_mesh
1146 : : !!***
1147 : :
1148 : : !----------------------------------------------------------------------
1149 : :
1150 : : !!****f* m_migdal_eliashberg/get_lambda_iso_iw
1151 : : !! NAME
1152 : : !! get_lambda_iso_iw
1153 : : !!
1154 : : !! FUNCTION
1155 : : !! Compute isotropic lambda along the imaginary axis
1156 : : !!
1157 : : !! NOTES
1158 : : !! This routine currently requires with_cplex=1, i.e. squared e-ph matrix
1159 : : !! elements precomputed and stored in gqk%my_g2. The with_cplex=0 on-demand
1160 : : !! Wannier backend implemented in get_a2fw is not yet supported here.
1161 : : !!
1162 : : !! INPUTS
1163 : : !!
1164 : : !! OUTPUT
1165 : : !!
1166 : : !! SOURCE
1167 : :
1168 : : subroutine get_lambda_iso_iw(gstore, nw, imag_w, lambda)
1169 : :
1170 : : !Arguments ------------------------------------
1171 : : class(gstore_t),intent(inout) :: gstore
1172 : : integer,intent(in) :: nw
1173 : : real(dp),intent(in) :: imag_w(nw)
1174 : : real(dp),intent(out) :: lambda(nw)
1175 : :
1176 : : !Local variables-------------------------------
1177 : : integer :: my_is, my_ik, my_iq, my_ip, in_k, im_kq, ierr, nb_k, nb_kq
1178 : : real(dp) :: g2, wqnu, weight_k, weight_q
1179 : : !arrays
1180 : : real(dp) :: qpt(3)
1181 : : real(dp),allocatable :: dbl_delta_q(:,:,:), g2_pmnk(:,:,:,:)
1182 : : !----------------------------------------------------------------------
1183 : :
1184 : : ABI_CHECK(gstore%qzone == "bz", "get_lambda_iso_iw assumes qzone == `bz`")
1185 : : ABI_CHECK(gstore%with_cplex == 1, "get_lambda_iso_iw requires squared e-ph matrix elements in memory (with_cplex=1)")
1186 : : !if (gstore%check_cplex_qkzone_gmode(cplex1, "bz", kzone, gmode, kfilter) result(ierr)
1187 : :
1188 : : lambda = zero
1189 : : do my_is=1,gstore%my_nspins
1190 : : associate (gqk => gstore%gqk(my_is))
1191 : : ABI_CHECK(allocated(gqk%my_g2), "my_g2 is not allocated")
1192 : : ABI_CHECK(allocated(gqk%my_wnuq), "my_wnuq is not allocated")
1193 : :
1194 : : nb_k = gqk%nb_k; nb_kq = gqk%nb_kq
1195 : : ABI_CHECK_IEQ(nb_k, nb_kq, "gqk_dbldelta_qpt does not support nb_k != nb_kq")
1196 : :
1197 : : ! Weights for delta(e_{m k+q}) delta(e_{n k}) for my list of k-points.
1198 : : ABI_MALLOC(dbl_delta_q, (nb_kq, nb_k, gqk%my_nk))
1199 : : ABI_MALLOC(g2_pmnk, (gqk%my_npert, nb_kq, nb_k, gqk%my_nk))
1200 : :
1201 : : do my_iq=1,gqk%my_nq
1202 : : ! Compute integration weights for the double delta.
1203 : : call gqk%dbldelta_qpt(my_iq, gstore, gstore%dtset%eph_intmeth, gstore%dtset%eph_fsmear, qpt, weight_q, dbl_delta_q)
1204 : :
1205 : : ! Copy data to improve memory access in the loops below.
1206 : : g2_pmnk = gqk%my_g2(:,:,my_iq,:,:)
1207 : :
1208 : : do my_ik=1,gqk%my_nk
1209 : : weight_k = gqk%my_wtk(my_ik)
1210 : : do in_k=1,nb_k
1211 : : do im_kq=1,nb_kq
1212 : : do my_ip=1,gqk%my_npert
1213 : : g2 = g2_pmnk(my_ip, im_kq, in_k, my_ik)
1214 : : ! TODO: handle wqnu ~ 0
1215 : : wqnu = gqk%my_wnuq(my_ip, my_iq)
1216 : : lambda(:) = lambda(:) + &
1217 : : two * wqnu / (imag_w(:) ** 2 + wqnu ** 2) * g2 * weight_k * weight_q * dbl_delta_q(im_kq, in_k, my_ik)
1218 : : end do
1219 : : end do
1220 : : end do
1221 : : end do
1222 : : end do ! my_iq
1223 : :
1224 : : ABI_FREE(dbl_delta_q)
1225 : : ABI_FREE(g2_pmnk)
1226 : : end associate
1227 : : end do ! my_is
1228 : :
1229 : : ! Take into account collinear spin
1230 : : lambda = lambda * (two / (gstore%nsppol * gstore%dtset%nspinor))
1231 : : call xmpi_sum(lambda, gstore%comm, ierr)
1232 : :
1233 : : end subroutine get_lambda_iso_iw
1234 : : !!***
1235 : :
1236 [ # # # # : 0 : end module m_migdal_eliashberg
# # # # #
# # # # #
# # # # #
# # # ]
1237 : : !!***
|
