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