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 : !!***
|
