Branch data Line data Source code
1 : : !!****m* ABINIT/m_eliashberg_1d
2 : : !! NAME
3 : : !! m_eliashberg_1d
4 : : !!
5 : : !! FUNCTION
6 : : !! Solve the Eliashberg equations in the isotropic case
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2008-2026 ABINIT group (MVer)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! SOURCE
15 : :
16 : : #if defined HAVE_CONFIG_H
17 : : #include "config.h"
18 : : #endif
19 : :
20 : : #include "abi_common.h"
21 : :
22 : : module m_eliashberg_1d
23 : :
24 : : use defs_basis
25 : : use defs_elphon
26 : : use m_errors
27 : : use m_abicore
28 : : use m_io_tools
29 : : use m_abicore
30 : :
31 : : use m_numeric_tools, only : simpson_int
32 : :
33 : : implicit none
34 : :
35 : : private
36 : : !!***
37 : :
38 : : public :: eliashberg_1d
39 : : !!***
40 : :
41 : : contains
42 : : !!***
43 : :
44 : : !!****f* m_eliashberg_1d/eliashberg_1d
45 : : !!
46 : : !! NAME
47 : : !! eliashberg_1d
48 : : !!
49 : : !! FUNCTION
50 : : !! Solve the Eliashberg equations in the isotropic case
51 : : !! First the linearized case, which allows the estimation of Tc
52 : : !! then the full case which gives the gap as a function of temperature.
53 : : !!
54 : : !! INPUTS
55 : : !! a2f_1d = alpha^2F function averaged over the FS (only energy dependence)
56 : : !! elph_ds = datastructure with phonon matrix elements
57 : : !! gkk2 = gkk2 matrix elements on full FS grid for each phonon mode
58 : : !! natom = number of atoms
59 : : !!
60 : : !! OUTPUT
61 : : !!
62 : : !! NOTES
63 : : !! na2f = number of frequency points for alpha^2F function
64 : : !!
65 : : !! SOURCE
66 : :
67 : 15 : subroutine eliashberg_1d(a2f_1d,elph_ds,mustar)
68 : :
69 : : !Arguments ------------------------------------
70 : : !scalars
71 : : real(dp),intent(in) :: mustar
72 : : type(elph_type),intent(in) :: elph_ds
73 : : !arrays
74 : : real(dp),intent(in) :: a2f_1d(elph_ds%na2f)
75 : :
76 : : !Local variables-------------------------------
77 : : ! for diagonalization of gammma matrix
78 : : ! output variables for gtdyn9+dfpt_phfrq
79 : : !scalars
80 : : integer :: iiter,imatsu
81 : : integer :: maxiter,nmatsu,unit_del,unit_lam,unit_z
82 : : real(dp) :: maxeigval,omega_cutoff
83 : : real(dp) :: tc
84 : : character(len=fnlen) :: fname
85 : : !arrays
86 : 15 : real(dp),allocatable :: delta_1d(:),lambda_1d(:),mm_1d(:,:),z_1d(:)
87 : :
88 : : ! *********************************************************************
89 : :
90 : 15 : call wrtout(std_out,'Solving the 1-D Eliashberg equation (isotropic case)',"COLL")
91 : :
92 [ + + ]: 15 : if (elph_ds%nsppol /= 1) then
93 : 1 : write(std_out,*) 'eliashberg_1d is not coded for nsppol > 1 yet'
94 : : return
95 : : end if
96 : :
97 : : !maximum number of iterations to find T_c
98 : 14 : maxiter=30
99 : :
100 : : !Fix nmatsu. Should add test in iiter loop to check if
101 : : !omega_cutoff is respected
102 : 14 : nmatsu = 50
103 : : !write(std_out,*) ' eliashberg_1d : nmatsu = ', nmatsu
104 : :
105 [ - + ]: 14 : ABI_MALLOC(lambda_1d,(-nmatsu:nmatsu))
106 [ - + ]: 14 : ABI_MALLOC(z_1d,(-nmatsu:nmatsu))
107 [ - + ]: 14 : ABI_MALLOC(delta_1d,(-nmatsu:nmatsu))
108 [ - + ]: 14 : ABI_MALLOC(mm_1d,(-nmatsu:nmatsu,-nmatsu:nmatsu))
109 : :
110 : 14 : unit_lam=get_unit()
111 [ + - ]: 14 : fname=trim(elph_ds%elph_base_name) // "_LAM"
112 : 14 : open (UNIT=unit_lam,FILE=fname,STATUS='REPLACE')
113 : 14 : unit_z=get_unit()
114 [ + - ]: 14 : fname=trim(elph_ds%elph_base_name) // "_Z"
115 : 14 : open (UNIT=unit_z,FILE=fname,STATUS='REPLACE')
116 : 14 : unit_del=get_unit()
117 [ + - ]: 14 : fname=trim(elph_ds%elph_base_name) // "_DEL"
118 : 14 : open (UNIT=unit_del,FILE=fname,STATUS='REPLACE')
119 : :
120 : : !
121 : : !1) use linearized Eliashberg equation to find Tc
122 : : !! \sum_j \mathbf{M}_{ij} \Delta_j = \zeta \cdot \Delta_i $ $i,j = 1 .. n_{\mathrm{Matsubara}}$
123 : : !$\zeta = 1$ gives T$_c$ $\beta = \frac{1}{\mathrm{T}}$ $\omega_i = (2 i + 1) \pi \mathrm{T}$
124 : : !! \mathbf{M}_{ij} = \frac{\pi}{\beta} \frac{\lambda (\omega_i - \omega_j)}{Z (\omega_i)}$
125 : : !! Z (\omega_i) = 1 + \frac{\pi}{\beta \omega_i} \sum_j \lambda(\omega_i - \omega_j) \mathrm{sgn}(\omega_j)$
126 : : !
127 : :
128 : : !initial guess for T$_c$ in Hartree (1Ha =3.067e5 K)
129 : 14 : tc = 0.0001
130 : : !
131 : : !big iterative loop
132 : : !
133 [ + + ]: 352 : do iiter=1,maxiter
134 : :
135 : 345 : omega_cutoff = (two*nmatsu+one) * pi * tc
136 : :
137 : : !
138 : : ! calculate array of lambda values
139 : : !
140 : 345 : call eli_lambda_1d (a2f_1d,elph_ds,lambda_1d,nmatsu,tc)
141 : 345 : write (unit_lam,'(a)') '#'
142 : 345 : write (unit_lam,'(a)') '# ABINIT package : lambda file'
143 : 345 : write (unit_lam,'(a)') '#'
144 : 345 : write (unit_lam,'(a,I10,a)') '# lambda_1d array containing 2*', nmatsu, '+1 Matsubara frequency points'
145 : 345 : write (unit_lam,'(a,E16.6,a,E16.6)') '# from ', -omega_cutoff, ' to ', omega_cutoff
146 : 345 : write (unit_lam,'(a)') '# lambda_1d is the frequency dependent coupling constant '
147 : 345 : write (unit_lam,'(a)') '# in the Eliashberg equations '
148 : 345 : write (unit_lam,'(a)') '#'
149 [ + + ]: 35190 : do imatsu=-nmatsu,nmatsu
150 : 35190 : write (unit_lam,*) imatsu,lambda_1d(imatsu)
151 : : end do
152 : 345 : write (unit_lam,*)
153 : :
154 : : !
155 : : ! calculate array of z values
156 : : !
157 : 345 : call eli_z_1d (lambda_1d,nmatsu,z_1d)
158 : 345 : write (unit_z,'(a)') '#'
159 : 345 : write (unit_z,'(a)') '# ABINIT package : Z file'
160 : 345 : write (unit_z,'(a)') '#'
161 : 345 : write (unit_z,'(a,I10,a)') '# z_1d array containing 2*', nmatsu, '+1 Matsubara frequency points'
162 : 345 : write (unit_z,'(a,E16.6,a,E16.6)') '# from ', -omega_cutoff, ' to ', omega_cutoff
163 : 345 : write (unit_z,'(a)') '# z_1d is the renormalization factor in the Eliashberg equations'
164 : 345 : write (unit_z,'(a)') '#'
165 [ + + ]: 35190 : do imatsu=-nmatsu,nmatsu
166 : 35190 : write (unit_z,*) imatsu,z_1d(imatsu)
167 : : end do
168 : :
169 : : ! !
170 : : ! ! apply M matrix until a maximal eigenvalue is found.
171 : : ! !
172 : : ! call eli_m_iter_1d (delta_1d,lambda_1d,maxeigval,nmatsu,z_1d)
173 : :
174 : : !
175 : : ! diagonalize M brute forcefully
176 : : !
177 : 345 : call eli_diag_m_1d(delta_1d,lambda_1d,maxeigval,mustar,nmatsu,tc,z_1d)
178 : :
179 : 345 : write (unit_del,'(a)') '#'
180 : 345 : write (unit_del,'(a)') '# eliashberg_1d : delta_1d = '
181 : 345 : write (unit_del,'(a)') '#'
182 : 345 : write (unit_del,'(a,i6,a)') '# delta_1d array containing 2*', nmatsu, '+1 Matsubara frequency points'
183 : 345 : write (unit_z,'(a,E16.6,a,E16.6)') '# from ', -omega_cutoff, ' to ', omega_cutoff
184 : 345 : write (unit_z,'(a)') '# delta_1d is the gap function in the Eliashberg equations'
185 : 345 : write (unit_z,'(a)') '#'
186 [ + + ]: 35190 : do imatsu=-nmatsu,nmatsu
187 : 35190 : write (unit_del,*) imatsu,delta_1d(imatsu)
188 : : end do
189 : 345 : write (unit_del,*)
190 : :
191 : : ! if eigenvalue is < 1 increase T
192 : : ! else if eigenvalue is > 1 decrease T
193 : : ! if eigenvalue ~= 1 stop
194 : : !
195 [ + + ]: 697 : if (abs(maxeigval-one) < tol8) then
196 : 7 : write(std_out,*) 'Eliashberg Tc found = ', tc, ' (Ha) = ', tc/kb_HaK, ' (K)'
197 : 7 : exit
198 [ + - ]: 338 : else if (maxeigval > 0.001_dp) then
199 : 338 : tc = tc * maxeigval
200 : : else
201 : 0 : write(std_out,*) 'maxeigval is very small'
202 : 0 : tc = tc * 1000.0_dp
203 : : end if
204 : :
205 : :
206 : : end do
207 : : !end iiter do
208 : :
209 [ + + ]: 14 : if (abs(maxeigval-one) > tol8) then
210 : 7 : write(std_out,*) 'eliashberg_1d : Tc not converged. ', maxeigval, ' /= 1'
211 : 7 : write(std_out,*) 'Eliashberg Tc nonetheless = ', tc, ' (Ha) = ', tc/kb_HaK, ' (K)'
212 : : end if
213 : :
214 [ - + ]: 14 : ABI_FREE(lambda_1d)
215 [ - + ]: 14 : ABI_FREE(z_1d)
216 [ - + ]: 14 : ABI_FREE(delta_1d)
217 [ - + ]: 14 : ABI_FREE(mm_1d)
218 : :
219 : 14 : close (UNIT=unit_z)
220 : 14 : close (UNIT=unit_lam)
221 : 14 : close (UNIT=unit_del)
222 : :
223 : 14 : write(std_out,*) ' eliashberg_1d : end '
224 : :
225 : : end subroutine eliashberg_1d
226 : : !!***
227 : :
228 : : !!****f* m_eliashberg_1d/eli_app_m_1d
229 : : !!
230 : : !! NAME
231 : : !! eli_app_m_1d
232 : : !!
233 : : !! FUNCTION
234 : : !! Apply the linearized Eliashberg matrix once to the input vector.
235 : : !!
236 : : !! INPUTS
237 : : !! lambda_1d = coupling constant as a function of frequency
238 : : !! nmatsu = number of Matsubara frequencies
239 : : !! tc = guess for critical temperature
240 : : !! z_1d = renormalization Z as a function of frequency
241 : : !!
242 : : !! SIDE EFFECTS
243 : : !! delta_1d = imaginary gap function as a function of frequency changed
244 : : !!
245 : : !! SOURCE
246 : :
247 : :
248 : : subroutine eli_app_m_1d (delta_1d,lambda_1d,nmatsu,z_1d)
249 : :
250 : : !Arguments ------------------------------------
251 : : !scalars
252 : : integer,intent(in) :: nmatsu
253 : : !arrays
254 : : real(dp),intent(in) :: lambda_1d(-nmatsu:nmatsu),z_1d(-nmatsu:nmatsu)
255 : : real(dp),intent(inout) :: delta_1d(-nmatsu:nmatsu)
256 : :
257 : : !Local variables-------------------------------
258 : : !scalars
259 : : integer :: imatsu,jmatsu,miguelflag
260 : : real(dp) :: zfact
261 : : !arrays
262 : : real(dp) :: delta_tmp(-nmatsu:nmatsu),freqfact(-nmatsu:nmatsu)
263 : :
264 : : ! *********************************************************************
265 : :
266 : : miguelflag = 0
267 : :
268 : :
269 : : do imatsu=-nmatsu,nmatsu
270 : : freqfact(imatsu) = one / abs(two*imatsu+one)
271 : : end do
272 : :
273 : : delta_tmp(:) = delta_1d(:)
274 : :
275 : : if (miguelflag == 1) then
276 : : do imatsu=-nmatsu,nmatsu
277 : : ! zfact = pi*tc / z_1d(imatsu)
278 : : zfact = one / z_1d(imatsu)
279 : :
280 : : do jmatsu=max(-nmatsu,-nmatsu+imatsu),min(nmatsu,nmatsu+imatsu)
281 : : delta_tmp(imatsu) = delta_tmp(imatsu) &
282 : : & + delta_1d(jmatsu) &
283 : : & * lambda_1d(imatsu-jmatsu) &
284 : : & * freqfact(jmatsu)
285 : : end do
286 : : delta_tmp(imatsu) = delta_tmp(imatsu)*zfact
287 : : end do
288 : :
289 : : else
290 : :
291 : : ! i < 0
292 : : do imatsu=-nmatsu,-1
293 : :
294 : : ! j < 0
295 : : do jmatsu=max(-nmatsu,-nmatsu+imatsu),-1
296 : : delta_tmp(imatsu) = delta_tmp(imatsu) &
297 : : & + lambda_1d(imatsu-jmatsu)*delta_1d(jmatsu)*freqfact(jmatsu) &
298 : : & - lambda_1d(imatsu-jmatsu)*delta_1d(imatsu)*freqfact(imatsu)
299 : : end do
300 : : ! j > 0
301 : : do jmatsu=0,min(nmatsu,nmatsu+imatsu)
302 : : delta_tmp(imatsu) = delta_tmp(imatsu) &
303 : : & + lambda_1d(imatsu-jmatsu)*delta_1d(jmatsu)*freqfact(jmatsu) &
304 : : & + lambda_1d(imatsu-jmatsu)*delta_1d(imatsu)*freqfact(imatsu)
305 : : end do
306 : :
307 : : end do
308 : :
309 : : ! i > 0
310 : : do imatsu=0,nmatsu
311 : :
312 : : ! j < 0
313 : : do jmatsu=max(-nmatsu,-nmatsu+imatsu),-1
314 : : delta_tmp(imatsu) = delta_tmp(imatsu) &
315 : : & + lambda_1d(imatsu-jmatsu)*delta_1d(jmatsu)*freqfact(jmatsu) &
316 : : & + lambda_1d(imatsu-jmatsu)*delta_1d(imatsu)*freqfact(imatsu)
317 : : end do
318 : : ! j > 0
319 : : do jmatsu=0,min(nmatsu,nmatsu+imatsu)
320 : : delta_tmp(imatsu) = delta_tmp(imatsu) &
321 : : & + lambda_1d(imatsu-jmatsu)*delta_1d(jmatsu)*freqfact(jmatsu) &
322 : : & - lambda_1d(imatsu-jmatsu)*delta_1d(imatsu)*freqfact(imatsu)
323 : : end do
324 : :
325 : : end do
326 : :
327 : : end if
328 : :
329 : : delta_1d(:) = delta_tmp(:)
330 : :
331 : : end subroutine eli_app_m_1d
332 : : !!***
333 : :
334 : : !!****f* m_eliashberg_1d/eli_diag_m_1d
335 : : !!
336 : : !! NAME
337 : : !! eli_diag_m_1d
338 : : !!
339 : : !! FUNCTION
340 : : !! diagonalize M matrix. Heavy and should be avoided for production.
341 : : !! Actually, since M is not symmetrical, diagonalize M^{t} M and
342 : : !! get right-eigenvalues and vectors
343 : : !!
344 : : !! INPUTS
345 : : !! lambda_1d = coupling constant as a function of frequency
346 : : !! mustar = Coulomb potential parameter in Eliashberg equation
347 : : !! nmatsu = number of Matsubara frequencies
348 : : !! tc = guess for critical temperature
349 : : !! z_1d = renormalization Z as a function of frequency
350 : : !!
351 : : !! OUTPUT
352 : : !! maxeigval = estimation for maximum eigenvalue of M
353 : : !!
354 : : !! SIDE EFFECTS
355 : : !! delta_1d = imaginary gap function as a function of frequency
356 : : !!
357 : : !! SOURCE
358 : :
359 : 345 : subroutine eli_diag_m_1d (delta_1d,lambda_1d,maxeigval,mustar,nmatsu,tc,z_1d)
360 : :
361 : : use m_linalg_interfaces
362 : :
363 : : !Arguments ------------------------------------
364 : : !scalars
365 : : integer,intent(in) :: nmatsu
366 : : real(dp),intent(in) :: mustar,tc
367 : : real(dp),intent(out) :: maxeigval
368 : : !arrays
369 : : real(dp),intent(in) :: lambda_1d(-nmatsu:nmatsu),z_1d(-nmatsu:nmatsu)
370 : : real(dp),intent(inout) :: delta_1d(-nmatsu:nmatsu)
371 : :
372 : : !Local variables-------------------------------
373 : : !scalars
374 : : integer :: imatsu,info,jmatsu,kmatsu,lwork,tmiguel
375 : : real(dp) :: si,sj,sqtimat,sqtjmat
376 : : !arrays
377 : 690 : real(dp) :: mtm_eig(2*nmatsu+1),symm_mtm(-nmatsu:nmatsu,-nmatsu:nmatsu)
378 : 690 : real(dp) :: work(3*(2*nmatsu+1))
379 : :
380 : : ! *********************************************************************
381 : :
382 : 345 : tmiguel = 0
383 : :
384 : : if (tmiguel == 1) then
385 : : do imatsu=-nmatsu,nmatsu
386 : : do jmatsu=-nmatsu,nmatsu
387 : :
388 : : symm_mtm(imatsu,jmatsu) = zero
389 : : do kmatsu=max(-nmatsu+imatsu,-nmatsu+jmatsu,-nmatsu),min(nmatsu+imatsu,nmatsu+jmatsu,nmatsu)
390 : : symm_mtm(imatsu,jmatsu) = symm_mtm(imatsu,jmatsu) &
391 : : & + lambda_1d(kmatsu-imatsu)*lambda_1d(kmatsu-jmatsu) &
392 : : & / ( z_1d(kmatsu)*z_1d(kmatsu) )
393 : : end do
394 : : ! symm_mtm(imatsu,jmatsu) = symm_mtm(imatsu,jmatsu) / ((two*imatsu+one)*(two*jmatsu+one))
395 : : symm_mtm(imatsu,jmatsu) = symm_mtm(imatsu,jmatsu) * pi * tc * pi * tc
396 : : end do
397 : : end do
398 : :
399 : : else
400 : :
401 [ + + + + ]: 3554535 : symm_mtm(:,:) = -mustar
402 : :
403 : : si = -one
404 [ + + ]: 35190 : do imatsu=-nmatsu,nmatsu
405 : 34845 : sqtimat = one / sqrt(two*abs(imatsu)+one)
406 [ + + ]: 34845 : if (imatsu == 0) si = one
407 : 34845 : sj = -one
408 [ + + ]: 2674785 : do jmatsu=max(-nmatsu,-nmatsu+imatsu),min(nmatsu,nmatsu+imatsu)
409 : 2639595 : sqtjmat = one / sqrt(two*abs(jmatsu)+one)
410 [ + + ]: 2639595 : if (jmatsu == 0) sj = one
411 : :
412 : : symm_mtm(imatsu,jmatsu) = symm_mtm(imatsu,jmatsu) &
413 : 2639595 : & + lambda_1d(imatsu-jmatsu)*sqtimat*sqtjmat
414 : :
415 : : symm_mtm(imatsu,imatsu) = symm_mtm(imatsu,imatsu) &
416 : 2674440 : & - lambda_1d(imatsu-jmatsu)*si*sj*sqtimat*sqtimat
417 : : end do
418 : : end do
419 : :
420 : : end if
421 : :
422 : 345 : lwork = 3*(2*nmatsu+1)
423 : 345 : call DSYEV('V', 'U', 2*nmatsu+1, symm_mtm, 2*nmatsu+1, mtm_eig, work, lwork, info )
424 : :
425 : 345 : write(std_out,*) 'last eigenvalues = '
426 : 345 : write(std_out,*) mtm_eig(2*nmatsu-9:2*nmatsu+1)
427 : :
428 [ + + ]: 35190 : do imatsu=-nmatsu,nmatsu
429 : 35190 : delta_1d(imatsu) = symm_mtm(imatsu,nmatsu)*sqrt(two*abs(imatsu)+one)
430 : : end do
431 : :
432 : 345 : maxeigval = mtm_eig(2*nmatsu+1)
433 : :
434 : 345 : end subroutine eli_diag_m_1d
435 : : !!***
436 : :
437 : : !!****f* m_eliashberg_1d/eli_lambda_1d
438 : : !!
439 : : !! NAME
440 : : !! eli_lambda_1d
441 : : !!
442 : : !! FUNCTION
443 : : !! In the solving of the 1D (energy only) Eliashberg equations, calculate
444 : : !! the lambda, which is the e-p coupling strength. See Allen and Mitrovic
445 : : !! Solid State Physics vol 37 ed Ehrenreich Seitz and Turnbull, p.45 [[cite:Allen1983c]]
446 : : !!
447 : : !! INPUTS
448 : : !! a2f_1d = 1D alpha2F function
449 : : !! elph_ds = elphon dataset
450 : : !! nmatsu = number of Matsubara frequencies
451 : : !! tc = guess for critical temperature
452 : : !!
453 : : !! OUTPUT
454 : : !! lambda_1d = coupling constant as a function of frequency
455 : : !!
456 : : !! NOTES
457 : : !! lambda is used at points which are differences of Matsubara freqs,
458 : : !! and hence is tabulated on points going through 0.
459 : : !!
460 : : !! SOURCE
461 : :
462 : 345 : subroutine eli_lambda_1d (a2f_1d,elph_ds,lambda_1d,nmatsu,tc)
463 : :
464 : : !Arguments ------------------------------------
465 : : !scalars
466 : : integer,intent(in) :: nmatsu
467 : : real(dp),intent(in) :: tc
468 : : type(elph_type),intent(in) :: elph_ds
469 : : !arrays
470 : : real(dp),intent(in) :: a2f_1d(elph_ds%na2f)
471 : : real(dp),intent(out) :: lambda_1d(-nmatsu:nmatsu)
472 : :
473 : : !Local variables-------------------------------
474 : : !scalars
475 : : integer :: imatsu,iomega
476 : : real(dp) :: nu_matsu,nu_matsu2,omega,domega
477 : : !arrays
478 : 690 : real(dp) :: lambda_int(elph_ds%na2f),tmplambda(elph_ds%na2f)
479 : :
480 : : ! *********************************************************************
481 : : !
482 : : !MG: the step should be calculated locally using nomega and the extrema of the spectrum.
483 : : !One should not rely on previous calls for the setup of elph_ds%domega
484 : : !I will remove elph_ds%domega since mka2f.F90 will become a method of gamma_t
485 : 345 : domega =elph_ds%domega
486 : :
487 [ + + ]: 35190 : do imatsu=-nmatsu,nmatsu
488 : 34845 : nu_matsu = (two*imatsu)*pi*tc
489 : 34845 : nu_matsu2 = nu_matsu*nu_matsu
490 : :
491 [ + + ]: 13972845 : tmplambda(:) = zero
492 : 34845 : omega=domega
493 [ + + ]: 13938000 : do iomega=2,elph_ds%na2f
494 : 13903155 : tmplambda(iomega) = a2f_1d(iomega) * two * omega / (nu_matsu2 + omega*omega)
495 : 13938000 : omega=omega+domega
496 : : end do
497 : 34845 : call simpson_int(elph_ds%na2f,domega,tmplambda,lambda_int)
498 : :
499 : 35190 : lambda_1d(imatsu) = lambda_int(elph_ds%na2f)
500 : : end do
501 : :
502 : 345 : end subroutine eli_lambda_1d
503 : : !!***
504 : :
505 : : !!****f* m_eliashberg_1d/eli_m_iter_1d
506 : : !! NAME
507 : : !! eli_m_iter_1d
508 : : !!
509 : : !! FUNCTION
510 : : !! Find largest eigenvalue of M matrix, to deduce superconducting Tc
511 : : !!
512 : : !! INPUTS
513 : : !! lambda_1d = coupling constant as a function of frequency
514 : : !! nmatsu = number of Matsubara frequencies
515 : : !! tc = guess for critical temperature
516 : : !! z_1d = renormalization Z as a function of frequency
517 : : !!
518 : : !! OUTPUT
519 : : !! maxeigval = estimation for maximum eigenvalue of M
520 : : !!
521 : : !! SIDE EFFECTS
522 : : !! delta_1d = imaginary gap function as a function of frequency
523 : : !!
524 : : !! SOURCE
525 : :
526 : : subroutine eli_m_iter_1d (delta_1d,lambda_1d,maxeigval,nmatsu,z_1d)
527 : :
528 : : !Arguments ------------------------------------
529 : : !scalars
530 : : integer,intent(in) :: nmatsu
531 : : real(dp),intent(out) :: maxeigval
532 : : !arrays
533 : : real(dp),intent(in) :: lambda_1d(-nmatsu:nmatsu),z_1d(-nmatsu:nmatsu)
534 : : real(dp),intent(inout) :: delta_1d(-nmatsu:nmatsu)
535 : :
536 : : !Local variables-------------------------------
537 : : !scalars
538 : : integer :: iiterm,imatsu,nfilter,ngeteig
539 : : real(dp) :: dnewnorm,dnorm,fact
540 : : !arrays
541 : : real(dp) :: delta_old(-nmatsu:nmatsu)
542 : :
543 : : ! *********************************************************************
544 : :
545 : : nfilter = 10
546 : : ngeteig = 10
547 : :
548 : :
549 : : !
550 : : !1) apply M matrix enough times to filter out largest eigenvalue
551 : : !
552 : : do iiterm=1,nfilter
553 : : call eli_app_m_1d (delta_1d,lambda_1d,nmatsu,z_1d)
554 : :
555 : : ! DEBUG
556 : : ! dnorm=zero
557 : : ! do imatsu=-nmatsu,nmatsu
558 : : ! dnorm = dnorm + delta_1d(imatsu)*delta_1d(imatsu)/(two*imatsu+one)
559 : : ! end do
560 : : ! dnorm = sqrt(dnorm)
561 : : ! write(std_out,*) 'eli_m_iter_1d : dnorm ', dnorm
562 : : ! ENDDEBUG
563 : : end do
564 : :
565 : : !
566 : : !2) calculate norm
567 : : !
568 : : dnorm=zero
569 : : do imatsu=-nmatsu,nmatsu
570 : : dnorm = dnorm + delta_1d(imatsu)*delta_1d(imatsu)/abs(two*imatsu+one)
571 : : end do
572 : : dnorm = sqrt(dnorm)
573 : :
574 : : !normalize delta_1d
575 : : delta_1d(:) = delta_1d(:) / dnorm
576 : :
577 : : delta_old = delta_1d
578 : :
579 : : !DEBUG
580 : : !dnewnorm=zero
581 : : !do imatsu=-nmatsu,nmatsu
582 : : !dnewnorm = dnewnorm + delta_1d(imatsu)*delta_1d(imatsu)/abs(two*imatsu+one)
583 : : !end do
584 : : !dnewnorm = sqrt(dnewnorm)
585 : : !write(std_out,*) 'eli_m_iter_1d : dnewnorm1 ', dnewnorm
586 : : !ENDDEBUG
587 : :
588 : : !
589 : : !3) re-apply M matrix ngeteig times
590 : : !
591 : : do iiterm=1,ngeteig
592 : : call eli_app_m_1d (delta_1d,lambda_1d,nmatsu,z_1d)
593 : : ! DEBUG
594 : : ! dnewnorm=zero
595 : : ! do imatsu=-nmatsu,nmatsu
596 : : ! dnewnorm = dnewnorm + delta_1d(imatsu)*delta_1d(imatsu)/abs(two*imatsu+one)
597 : : ! end do
598 : : ! dnewnorm = sqrt(dnewnorm)
599 : : ! write(std_out,*) 'eli_m_iter_1d : dnewnorm ', dnewnorm
600 : : !
601 : : ! do imatsu=-nmatsu,nmatsu
602 : : ! write (112,*) imatsu,delta_1d(imatsu)/delta_old(imatsu)
603 : : ! end do
604 : : ! write (112,*)
605 : : ! delta_old = delta_1d
606 : : ! ENDDEBUG
607 : : end do
608 : :
609 : : !
610 : : !4) calculate new norm and estimate eigenvalue
611 : : !
612 : : dnewnorm=zero
613 : : do imatsu=-nmatsu,nmatsu
614 : : dnewnorm = dnewnorm + delta_1d(imatsu)*delta_1d(imatsu)/abs(two*imatsu+one)
615 : : end do
616 : : dnewnorm = sqrt(dnewnorm)
617 : :
618 : : !maxeigval = exp ( log(dnewnorm/dnorm) / ngeteig )
619 : : maxeigval = exp ( log(dnewnorm) / ngeteig )
620 : :
621 : : write(std_out,*) 'eli_m_iter_1d : maxeigval =' , maxeigval
622 : : !fact = exp(-log(maxeigval) * (ngeteig+nfilter))
623 : : fact = exp(-log(maxeigval) * (ngeteig))
624 : : do imatsu=-nmatsu,nmatsu
625 : : delta_1d(imatsu) = delta_1d(imatsu) * fact
626 : : end do
627 : :
628 : : end subroutine eli_m_iter_1d
629 : : !!***
630 : :
631 : : !!****f* m_eliashberg_1d/eli_z_1d
632 : : !!
633 : : !! NAME
634 : : !! eli_z_1d
635 : : !!
636 : : !! FUNCTION
637 : : !! In the solving of the 1D (energy only) Eliashberg equations, calculate
638 : : !! the Z function, which is the renormalization factor. See Allen and Mitrovic
639 : : !! Solid State Physics vol 37 ed Ehrenreich Seitz and Turnbull [[cite:Allen1983c]]
640 : : !!
641 : : !! INPUTS
642 : : !! lambda_1d = coupling constant as a function of frequency
643 : : !! nmatsu = number of Matsubara frequencies
644 : : !!
645 : : !! OUTPUT
646 : : !! z_1d = renormalizing Z as a function of frequency
647 : : !!
648 : : !! NOTES
649 : : !! Because Z only depends on lambda(n-n'), and lambda(omega)
650 : : !! is an even function, Z is symmetrical in n and -n
651 : : !! hence only calculate for n>0 and complete the rest
652 : : !!
653 : : !! SOURCE
654 : :
655 : 345 : subroutine eli_z_1d (lambda_1d,nmatsu,z_1d)
656 : :
657 : : !Arguments ------------------------------------
658 : : !scalars
659 : : integer,intent(in) :: nmatsu
660 : : !arrays
661 : : real(dp),intent(in) :: lambda_1d(-nmatsu:nmatsu)
662 : : real(dp),intent(out) :: z_1d(-nmatsu:nmatsu)
663 : :
664 : : !Local variables-------------------------------
665 : : !scalars
666 : : integer :: imatsu,jmatsu
667 : :
668 : : ! *********************************************************************
669 : :
670 : :
671 [ + + ]: 17940 : do imatsu=0,nmatsu
672 : :
673 : 17595 : z_1d(imatsu) = zero
674 : : ! count $\mathrm{sign}(omega_{Matsubara})$
675 [ + + ]: 457470 : do jmatsu=-nmatsu+imatsu,-1
676 : 457470 : z_1d(imatsu) = z_1d(imatsu) - lambda_1d(imatsu-jmatsu)
677 : : end do
678 [ + + ]: 914940 : do jmatsu=0,nmatsu
679 : 914940 : z_1d(imatsu) = z_1d(imatsu) + lambda_1d(imatsu-jmatsu)
680 : : end do
681 : :
682 : : ! NOTE: the pi*Tc factor in Z cancels the one in the Matsubara frequency.
683 : 17595 : z_1d(imatsu) = one + z_1d(imatsu) / (two*imatsu+one)
684 : 17940 : z_1d(-imatsu) = z_1d(imatsu)
685 : : end do
686 : :
687 : 345 : end subroutine eli_z_1d
688 : : !!***
689 : :
690 : : end module m_eliashberg_1d
691 : : !!***
|
