TABLE OF CONTENTS

ABINIT/m_migdal_eliashberg [ Modules ]

[ Top ] [ Modules ]

NAME

 m_migdal_eliashberg

FUNCTION

COPYRIGHT

  Copyright (C) 2008-2024 ABINIT group (MG)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

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 !, tolower, itoa, ftoa, ktoa, ltoa, strcat
37  !use m_numeric_tools,  only : arth, get_diag
38  use m_copy,            only : alloc_copy
39  use defs_datatypes ,   only : ebands_t
40  !use m_kpts,           only : kpts_timrev_from_kptopt
41  use m_ebands,          only : edos_t, ebands_get_edos
42  use m_gstore,          only : gstore_t
43 
44  implicit none
45 
46  private
47 
48  public :: migdal_eliashberg_iso
49  !public :: migdal_eliashberg_aniso

m_migdal_eliashberg/iso_solver_free [ Functions ]

[ Top ] [ m_migdal_eliashberg ] [ Functions ]

NAME

 iso_solver_free

FUNCTION

INPUTS

OUTPUT

SOURCE

118 subroutine iso_solver_free(iso)
119 
120 !Arguments ------------------------------------
121  class(iso_solver_t),intent(inout) :: iso
122 !----------------------------------------------------------------------
123 
124  ABI_SFREE(iso%delta_iw)
125  ABI_SFREE(iso%zeta_iw)
126  ABI_SFREE(iso%prev_delta_iw)
127  ABI_SFREE(iso%prev_zeta_iw)
128  ABI_SFREE(iso%delta_iw_mix)
129 
130 end subroutine iso_solver_free

m_migdal_eliashberg/iso_solver_solve [ Functions ]

[ Top ] [ m_migdal_eliashberg ] [ Functions ]

NAME

 iso_solver_solve

FUNCTION

INPUTS

OUTPUT

SOURCE

147 subroutine iso_solver_solve(iso, itemp, kt, niw, imag_w, lambda_ij)
148 
149 !Arguments ------------------------------------
150 !scalars
151  class(iso_solver_t),intent(inout) :: iso
152  integer,intent(in) :: itemp, niw
153  real(dp),intent(in) :: kt
154 !arrays
155  real(dp),intent(in) :: imag_w(niw), lambda_ij(2 * niw)
156 
157 !Local variables-------------------------------
158 !scalars
159  integer,parameter :: master = 0
160  integer :: nproc, my_rank, iter, ii, jj, converged
161  real(dp) :: rr
162 !arrays
163  real(dp),allocatable :: prev_vals(:)
164 
165 !----------------------------------------------------------------------
166 
167  ABI_UNUSED(lambda_ij)
168 
169  nproc = xmpi_comm_size(iso%comm); my_rank = xmpi_comm_rank(iso%comm)
170 
171  ABI_REMALLOC(iso%delta_iw_mix, (niw, iso%max_nmix))
172 
173  if (itemp == 1) then
174    ! Init values from scratch
175    ABI_CALLOC(iso%zeta_iw, (niw))
176    ABI_CALLOC(iso%prev_zeta_iw, (niw))
177    ABI_CALLOC(iso%delta_iw, (niw))
178    ABI_CALLOC(iso%prev_delta_iw, (niw))
179  else
180    ! Init values from previous temperature. TODO: May use spline
181    call alloc_copy(iso%zeta_iw, prev_vals)
182    ABI_RECALLOC(iso%zeta_iw, (niw))
183    ABI_MOVE_ALLOC(prev_vals, iso%prev_zeta_iw)
184    call alloc_copy(iso%delta_iw, prev_vals)
185    ABI_RECALLOC(iso%delta_iw, (niw))
186    ABI_MOVE_ALLOC(prev_vals, iso%prev_delta_iw)
187  end if
188 
189  converged = 0
190 iter_loop: do iter=1,iso%max_niter
191 
192    do ii=1,niw
193      !if (mod(ii, nproc) /= my_rank) cycle ! MPI parallelism inside comm
194      do jj=1,niw
195        rr = one / sqrt(imag_w(jj) ** 2 + iso%prev_delta_iw(jj) ** 2)
196        iso%zeta_iw(ii) = iso%zeta_iw(ii) + imag_w(jj) * rr  !* lambda(ii - jj)
197        iso%delta_iw(ii) = iso%delta_iw(ii) + rr * iso%prev_delta_iw(jj) !* (lambda(ii - jj) - mustar)
198      end do
199       iso%zeta_iw(ii) = one + pi * kt / imag_w(ii) * iso%zeta_iw(ii)
200       iso%delta_iw(ii) = pi * kt * iso%delta_iw(ii) / iso%zeta_iw(ii)
201    end do ! ii
202 
203    if (my_rank == master) then
204      ! Write SCF cycle to stdout.
205      ! Check for convergence.
206      converged = 0
207    end if
208 
209    if (converged == 2) exit iter_loop
210 
211    ! TODO: Mixing
212    iso%prev_zeta_iw = iso%zeta_iw
213    iso%prev_delta_iw = iso%delta_iw
214 
215  end do iter_loop
216 
217  ! Pade' to go to real axis
218  ! Compute Delta F
219  ! Compute QP DOS
220 
221  ! Write results to netcdf file
222  if (my_rank == master) then
223  end if
224 
225 end subroutine iso_solver_solve

m_migdal_eliashberg/iso_solver_t [ Types ]

[ Top ] [ m_migdal_eliashberg ] [ Types ]

NAME

 iso_solver_t

FUNCTION

NOTES

SOURCE

64 type, public :: iso_solver_t
65 
66   integer :: ntemp
67 
68   integer :: max_niter = -1
69   ! Maximum number of iterations.
70 
71   integer :: max_nmix = 4
72 
73   integer :: ncid
74 
75   integer :: comm = xmpi_undefined
76 
77   !integer :: niw = -1
78   ! Number of Matsubara frequencies.
79 
80   !real(dp) :: kt = -one
81   ! K * T in Ha.
82 
83   real(dp) :: tolerance = -one
84 
85   real(dp),allocatable :: zeta_iw(:), delta_iw(:)
86   real(dp),allocatable :: prev_zeta_iw(:), prev_delta_iw(:)
87   real(dp),allocatable :: delta_iw_mix(:,:)
88 
89 contains
90 
91   procedure :: free => iso_solver_free
92   ! Free dynamic memory
93 
94   procedure :: solve => iso_solver_solve
95 
96 end type iso_solver_t
97 
98 !private :: iso_solver_new

m_migdal_eliashberg/matsubara_mesh [ Functions ]

[ Top ] [ m_migdal_eliashberg ] [ Functions ]

NAME

 mastubara_mesh

FUNCTION

INPUTS

OUTPUT

SOURCE

385 subroutine matsubara_mesh(bosons_or_fermions, kt, wmax, niw, imag_w)
386 
387 !Arguments ------------------------------------
388 !scalars
389  character(len=*),intent(in) :: bosons_or_fermions
390  real(dp),intent(in) :: kt, wmax
391  integer,intent(out) :: niw
392  real(dp),allocatable,intent(out) :: imag_w(:)
393 
394 !Local variables-------------------------------
395  integer :: nn
396 !----------------------------------------------------------------------
397 
398  select case (bosons_or_fermions)
399  case ("bosons")
400    ! 2 n pi kT
401    niw = nint(wmax / (two * pi * kt)) + 1
402    ABI_MALLOC(imag_w, (niw))
403    do nn=0,niw-1
404      imag_w(nn + 1) = two * nn * pi * kt
405    end do
406 
407  case ("fermions")
408    ! 2 (n + 1) pi kT
409    niw = nint(wmax / (two * pi * kt))
410    ABI_MALLOC(imag_w, (niw))
411    do nn=0,niw-1
412      imag_w(nn + 1) = two * (nn  + 1) * pi * kt
413    end do
414 
415  case default
416    ABI_ERROR(sjoin("Wrong bosons_or_fermions:", bosons_or_fermions))
417  end select
418 
419 end subroutine matsubara_mesh

m_migdal_eliashberg/migdal_eliashberg_iso [ Functions ]

[ Top ] [ m_migdal_eliashberg ] [ Functions ]

NAME

 migdal_eliashber_iso

FUNCTION

INPUTS

OUTPUT

SOURCE

242 subroutine migdal_eliashberg_iso(gstore, dtset, dtfil)
243 
244 !Arguments ------------------------------------
245 !scalars
246  type(dataset_type),intent(in) :: dtset
247  type(datafiles_type),intent(in) :: dtfil
248  type(gstore_t),target,intent(inout) :: gstore
249 
250 !Local variables-------------------------------
251 !scalars
252  integer,parameter :: master = 0
253  integer :: nproc, my_rank, ierr, itemp, ntemp, niw, ncid !, my_nshiftq, nsppol !, iq_glob, ik_glob, ii ! out_nkibz,
254  integer :: edos_intmeth
255  !integer :: spin, natom3, cnt !, band, ib, nb, my_ik, my_iq, my_is
256  !integer :: ik_ibz, ik_bz, ebands_timrev
257  !integer :: iq_bz, iq_ibz !, ikq_ibz, ikq_bz
258  !integer :: ncid, spin_ncid, ncerr, gstore_fform
259  integer :: phmesh_size !, iw
260  real(dp) :: kt, wmax, cpu, wall, gflops
261  real(dp) :: edos_step, edos_broad !, sigma, ecut, eshift, eig0nk
262  !character(len=5000) :: msg
263  class(crystal_t),pointer :: cryst
264  class(ebands_t),pointer :: ebands
265  !class(ifc_type),target,intent(in) :: ifc
266  !type(gqk_t),pointer :: gqk
267  type(iso_solver_t) :: iso
268  type(edos_t) :: edos
269 !arrays
270  real(dp),allocatable :: ktmesh(:), lambda_ij(:), imag_w(:), imag_2w(:), phmesh(:), a2fw(:)
271 
272 !----------------------------------------------------------------------
273 
274  nproc = xmpi_comm_size(gstore%comm); my_rank = xmpi_comm_rank(gstore%comm)
275 
276  call wrtout(std_out, " Solving isotropic Migdal-Eliashberg equations on the imaginary axis", pre_newlines=2)
277  call cwtime(cpu, wall, gflops, "start")
278 
279  cryst => gstore%cryst; ebands => gstore%ebands
280  !natom3 = 3 * cryst%natom; nsppol = ebands%nsppol
281 
282  ! Consistency check
283  ierr = 0
284  ABI_CHECK_NOSTOP(gstore%qzone == "bz", "qzone == 'bz' is required", ierr)
285  ABI_CHECK_NOSTOP(gstore%gqk(1)%cplex == 1, "cplex == 1 is required", ierr)
286  ABI_CHECK(ierr == 0, "Wrong gstore object for migdal_eliashberg_iso. See messages above")
287 
288  ! Compute electron DOS.
289  edos_intmeth = 2; if (dtset%prtdos /= 0) edos_intmeth = dtset%prtdos
290  edos_step = dtset%dosdeltae; edos_broad = dtset%tsmear
291  edos_step = 0.01 * eV_Ha; edos_broad = 0.3 * eV_Ha
292  edos = ebands_get_edos(ebands, cryst, edos_intmeth, edos_step, edos_broad, gstore%comm)
293 
294  !! Store DOS per spin channel
295  !n0(:) = edos%gef(1:edos%nsppol)
296  if (my_rank == master) then
297    call edos%print(unit=std_out)
298    !call edos%print(unit=ab_out)
299    !path = strcat(dtfil%filnam_ds(4), "_EDOS")
300    !call wrtout(ab_out, sjoin("- Writing electron DOS to file:", path, ch10))
301    !call edos%write(path)
302  end if
303 
304  ! Compute phonon frequency mesh.
305  call gstore%ifc%get_phmesh(dtset%ph_wstep, phmesh_size, phmesh)
306 
307  ! Compute and store my phonon quantities
308  call gstore%calc_my_phonons(store_phdispl=.False.)
309 
310  ! Compute Eliashberg function a2F(w)
311  ABI_MALLOC(a2fw, (phmesh_size))
312  call gstore%get_a2fw(dtset, phmesh_size, phmesh, a2fw)
313 
314  ncid = nctk_noid
315  if (my_rank == master) then
316    NCF_CHECK(nctk_open_create(ncid, strcat(dtfil%filnam_ds(4), "_ISOME.nc") , xmpi_comm_self))
317    !write(777, *)"# phmesh (meV), a2fw"
318    !do iw=1, phmesh_size
319    !  write(777, *) phmesh(iw) * Ha_meV, a2fw(iw) / (edos%gef(0) / two)
320    !end do
321    !close(777)
322  end if
323 
324  ABI_FREE(a2fw)
325  ABI_FREE(phmesh)
326  call edos%free()
327 
328  call dtset%get_ktmesh(ntemp, ktmesh)
329  !NVHPC and LLVM don't like using this constructor because allocatable arrays aren't set.
330 #if defined FC_NVHPC || defined FC_LLVM
331   iso%ntemp=ntemp
332   iso%max_niter=10
333   iso%tolerance=tol10
334   iso%ncid=ncid
335   iso%comm=gstore%comm
336 #else
337  iso = iso_solver_t(ntemp=ntemp, max_niter=10, tolerance=tol10, ncid=ncid, comm=gstore%comm)
338 #endif
339 
340  do itemp=1,ntemp
341    ! Generate Matsubara mesh for this T with cutoff wmax.
342    kt = ktmesh(itemp)
343    wmax = one
344    call matsubara_mesh("bosons", kt, wmax, niw, imag_w)
345 
346    ! Compute lambda(w_i - w_j)
347    ABI_MALLOC(lambda_ij, (2 * niw))
348    ABI_MALLOC(imag_2w, (2 * niw))
349 
350    !call wrtout(std_out, " Computing lambda_iso_iw...")
351    !call gstore%get_lambda_iso_iw(dtset, 2 * niw, imag_2w, lambda_ij)
352    ABI_FREE(imag_2w)
353 
354    !call iso%solve(itemp, kt, niw, imag_w, lambda_ij)
355    ABI_FREE(lambda_ij)
356    ABI_FREE(imag_w)
357  end do ! itemp
358 
359  ABI_FREE(ktmesh)
360  call iso%free()
361 
362  if (my_rank == master) then
363    NCF_CHECK(nf90_close(ncid))
364  end if
365 
366  call cwtime_report(" migdal_eliashberg_iso:", cpu, wall, gflops)
367 
368 end subroutine migdal_eliashberg_iso