TABLE OF CONTENTS
- ABINIT/m_eliashberg_1d
- m_eliashberg_1d/eli_app_m_1d
- m_eliashberg_1d/eli_diag_m_1d
- m_eliashberg_1d/eli_lambda_1d
- m_eliashberg_1d/eli_m_iter_1d
- m_eliashberg_1d/eli_z_1d
- m_eliashberg_1d/eliashberg_1d
ABINIT/m_eliashberg_1d [ Modules ]
NAME
m_eliashberg_1d
FUNCTION
Solve the Eliashberg equations in the isotropic case
COPYRIGHT
Copyright (C) 2008-2022 ABINIT group (MVer) 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
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
m_eliashberg_1d/eli_app_m_1d [ Functions ]
[ Top ] [ m_eliashberg_1d ] [ Functions ]
NAME
eli_app_m_1d
FUNCTION
Apply the linearized Eliashberg matrix once to the input vector.
INPUTS
lambda_1d = coupling constant as a function of frequency nmatsu = number of Matsubara frequencies tc = guess for critical temperature z_1d = renormalization Z as a function of frequency
SIDE EFFECTS
delta_1d = imaginary gap function as a function of frequency changed
SOURCE
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
m_eliashberg_1d/eli_diag_m_1d [ Functions ]
[ Top ] [ m_eliashberg_1d ] [ Functions ]
NAME
eli_diag_m_1d
FUNCTION
diagonalize M matrix. Heavy and should be avoided for production. Actually, since M is not symmetrical, diagonalize M^{t} M and get right-eigenvalues and vectors
INPUTS
lambda_1d = coupling constant as a function of frequency mustar = Coulomb potential parameter in Eliashberg equation nmatsu = number of Matsubara frequencies tc = guess for critical temperature z_1d = renormalization Z as a function of frequency
OUTPUT
maxeigval = estimation for maximum eigenvalue of M
SIDE EFFECTS
delta_1d = imaginary gap function as a function of frequency
SOURCE
359 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 real(dp) :: mtm_eig(2*nmatsu+1),symm_mtm(-nmatsu:nmatsu,-nmatsu:nmatsu) 378 real(dp) :: work(3*(2*nmatsu+1)) 379 380 ! ********************************************************************* 381 382 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 symm_mtm(:,:) = -mustar 402 403 si = -one 404 do imatsu=-nmatsu,nmatsu 405 sqtimat = one / sqrt(two*abs(imatsu)+one) 406 if (imatsu == 0) si = one 407 sj = -one 408 do jmatsu=max(-nmatsu,-nmatsu+imatsu),min(nmatsu,nmatsu+imatsu) 409 sqtjmat = one / sqrt(two*abs(jmatsu)+one) 410 if (jmatsu == 0) sj = one 411 412 symm_mtm(imatsu,jmatsu) = symm_mtm(imatsu,jmatsu) & 413 & + lambda_1d(imatsu-jmatsu)*sqtimat*sqtjmat 414 415 symm_mtm(imatsu,imatsu) = symm_mtm(imatsu,imatsu) & 416 & - lambda_1d(imatsu-jmatsu)*si*sj*sqtimat*sqtimat 417 end do 418 end do 419 420 end if 421 422 lwork = 3*(2*nmatsu+1) 423 call DSYEV('V', 'U', 2*nmatsu+1, symm_mtm, 2*nmatsu+1, mtm_eig, work, lwork, info ) 424 425 write(std_out,*) 'last eigenvalues = ' 426 write(std_out,*) mtm_eig(2*nmatsu-9:2*nmatsu+1) 427 428 do imatsu=-nmatsu,nmatsu 429 delta_1d(imatsu) = symm_mtm(imatsu,nmatsu)*sqrt(two*abs(imatsu)+one) 430 end do 431 432 maxeigval = mtm_eig(2*nmatsu+1) 433 434 end subroutine eli_diag_m_1d
m_eliashberg_1d/eli_lambda_1d [ Functions ]
[ Top ] [ m_eliashberg_1d ] [ Functions ]
NAME
eli_lambda_1d
FUNCTION
In the solving of the 1D (energy only) Eliashberg equations, calculate the lambda, which is the e-p coupling strength. See Allen and Mitrovic Solid State Physics vol 37 ed Ehrenreich Seitz and Turnbull, p.45 [[cite:Allen1983c]]
INPUTS
a2f_1d = 1D alpha2F function elph_ds = elphon dataset nmatsu = number of Matsubara frequencies tc = guess for critical temperature
OUTPUT
lambda_1d = coupling constant as a function of frequency
NOTES
lambda is used at points which are differences of Matsubara freqs, and hence is tabulated on points going through 0.
SOURCE
462 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 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 domega =elph_ds%domega 486 487 do imatsu=-nmatsu,nmatsu 488 nu_matsu = (two*imatsu)*pi*tc 489 nu_matsu2 = nu_matsu*nu_matsu 490 491 tmplambda(:) = zero 492 omega=domega 493 do iomega=2,elph_ds%na2f 494 tmplambda(iomega) = a2f_1d(iomega) * two * omega / (nu_matsu2 + omega*omega) 495 omega=omega+domega 496 end do 497 call simpson_int(elph_ds%na2f,domega,tmplambda,lambda_int) 498 499 lambda_1d(imatsu) = lambda_int(elph_ds%na2f) 500 end do 501 502 end subroutine eli_lambda_1d
m_eliashberg_1d/eli_m_iter_1d [ Functions ]
[ Top ] [ m_eliashberg_1d ] [ Functions ]
NAME
eli_m_iter_1d
FUNCTION
Find largest eigenvalue of M matrix, to deduce superconducting Tc
INPUTS
lambda_1d = coupling constant as a function of frequency nmatsu = number of Matsubara frequencies tc = guess for critical temperature z_1d = renormalization Z as a function of frequency
OUTPUT
maxeigval = estimation for maximum eigenvalue of M
SIDE EFFECTS
delta_1d = imaginary gap function as a function of frequency
SOURCE
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
m_eliashberg_1d/eli_z_1d [ Functions ]
[ Top ] [ m_eliashberg_1d ] [ Functions ]
NAME
eli_z_1d
FUNCTION
In the solving of the 1D (energy only) Eliashberg equations, calculate the Z function, which is the renormalization factor. See Allen and Mitrovic Solid State Physics vol 37 ed Ehrenreich Seitz and Turnbull [[cite:Allen1983c]]
INPUTS
lambda_1d = coupling constant as a function of frequency nmatsu = number of Matsubara frequencies
OUTPUT
z_1d = renormalizing Z as a function of frequency
NOTES
Because Z only depends on lambda(n-n'), and lambda(omega) is an even function, Z is symmetrical in n and -n hence only calculate for n>0 and complete the rest
SOURCE
655 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 do imatsu=0,nmatsu 672 673 z_1d(imatsu) = zero 674 ! count $\mathrm{sign}(omega_{Matsubara})$ 675 do jmatsu=-nmatsu+imatsu,-1 676 z_1d(imatsu) = z_1d(imatsu) - lambda_1d(imatsu-jmatsu) 677 end do 678 do jmatsu=0,nmatsu 679 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 z_1d(imatsu) = one + z_1d(imatsu) / (two*imatsu+one) 684 z_1d(-imatsu) = z_1d(imatsu) 685 end do 686 687 end subroutine eli_z_1d
m_eliashberg_1d/eliashberg_1d [ Functions ]
[ Top ] [ m_eliashberg_1d ] [ Functions ]
NAME
eliashberg_1d
FUNCTION
Solve the Eliashberg equations in the isotropic case First the linearized case, which allows the estimation of Tc then the full case which gives the gap as a function of temperature.
INPUTS
a2f_1d = alpha^2F function averaged over the FS (only energy dependence) elph_ds = datastructure with phonon matrix elements gkk2 = gkk2 matrix elements on full FS grid for each phonon mode natom = number of atoms
OUTPUT
NOTES
na2f = number of frequency points for alpha^2F function
SOURCE
67 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 real(dp),allocatable :: delta_1d(:),lambda_1d(:),mm_1d(:,:),z_1d(:) 87 88 ! ********************************************************************* 89 90 call wrtout(std_out,'Solving the 1-D Eliashberg equation (isotropic case)',"COLL") 91 92 if (elph_ds%nsppol /= 1) then 93 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 maxiter=30 99 100 !Fix nmatsu. Should add test in iiter loop to check if 101 !omega_cutoff is respected 102 nmatsu = 50 103 !write(std_out,*) ' eliashberg_1d : nmatsu = ', nmatsu 104 105 ABI_MALLOC(lambda_1d,(-nmatsu:nmatsu)) 106 ABI_MALLOC(z_1d,(-nmatsu:nmatsu)) 107 ABI_MALLOC(delta_1d,(-nmatsu:nmatsu)) 108 ABI_MALLOC(mm_1d,(-nmatsu:nmatsu,-nmatsu:nmatsu)) 109 110 unit_lam=get_unit() 111 fname=trim(elph_ds%elph_base_name) // "_LAM" 112 open (UNIT=unit_lam,FILE=fname,STATUS='REPLACE') 113 unit_z=get_unit() 114 fname=trim(elph_ds%elph_base_name) // "_Z" 115 open (UNIT=unit_z,FILE=fname,STATUS='REPLACE') 116 unit_del=get_unit() 117 fname=trim(elph_ds%elph_base_name) // "_DEL" 118 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 tc = 0.0001 130 ! 131 !big iterative loop 132 ! 133 do iiter=1,maxiter 134 135 omega_cutoff = (two*nmatsu+one) * pi * tc 136 137 ! 138 ! calculate array of lambda values 139 ! 140 call eli_lambda_1d (a2f_1d,elph_ds,lambda_1d,nmatsu,tc) 141 write (unit_lam,'(a)') '#' 142 write (unit_lam,'(a)') '# ABINIT package : lambda file' 143 write (unit_lam,'(a)') '#' 144 write (unit_lam,'(a,I10,a)') '# lambda_1d array containing 2*', nmatsu, '+1 Matsubara frequency points' 145 write (unit_lam,'(a,E16.6,a,E16.6)') '# from ', -omega_cutoff, ' to ', omega_cutoff 146 write (unit_lam,'(a)') '# lambda_1d is the frequency dependent coupling constant ' 147 write (unit_lam,'(a)') '# in the Eliashberg equations ' 148 write (unit_lam,'(a)') '#' 149 do imatsu=-nmatsu,nmatsu 150 write (unit_lam,*) imatsu,lambda_1d(imatsu) 151 end do 152 write (unit_lam,*) 153 154 ! 155 ! calculate array of z values 156 ! 157 call eli_z_1d (lambda_1d,nmatsu,z_1d) 158 write (unit_z,'(a)') '#' 159 write (unit_z,'(a)') '# ABINIT package : Z file' 160 write (unit_z,'(a)') '#' 161 write (unit_z,'(a,I10,a)') '# z_1d array containing 2*', nmatsu, '+1 Matsubara frequency points' 162 write (unit_z,'(a,E16.6,a,E16.6)') '# from ', -omega_cutoff, ' to ', omega_cutoff 163 write (unit_z,'(a)') '# z_1d is the renormalization factor in the Eliashberg equations' 164 write (unit_z,'(a)') '#' 165 do imatsu=-nmatsu,nmatsu 166 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 call eli_diag_m_1d(delta_1d,lambda_1d,maxeigval,mustar,nmatsu,tc,z_1d) 178 179 write (unit_del,'(a)') '#' 180 write (unit_del,'(a)') '# eliashberg_1d : delta_1d = ' 181 write (unit_del,'(a)') '#' 182 write (unit_del,'(a,i6,a)') '# delta_1d array containing 2*', nmatsu, '+1 Matsubara frequency points' 183 write (unit_z,'(a,E16.6,a,E16.6)') '# from ', -omega_cutoff, ' to ', omega_cutoff 184 write (unit_z,'(a)') '# delta_1d is the gap function in the Eliashberg equations' 185 write (unit_z,'(a)') '#' 186 do imatsu=-nmatsu,nmatsu 187 write (unit_del,*) imatsu,delta_1d(imatsu) 188 end do 189 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 if (abs(maxeigval-one) < tol8) then 196 write(std_out,*) 'Eliashberg Tc found = ', tc, ' (Ha) = ', tc/kb_HaK, ' (K)' 197 exit 198 else if (maxeigval > 0.001_dp) then 199 tc = tc * maxeigval 200 else 201 write(std_out,*) 'maxeigval is very small' 202 tc = tc * 1000.0_dp 203 end if 204 205 206 end do 207 !end iiter do 208 209 if (abs(maxeigval-one) > tol8) then 210 write(std_out,*) 'eliashberg_1d : Tc not converged. ', maxeigval, ' /= 1' 211 write(std_out,*) 'Eliashberg Tc nonetheless = ', tc, ' (Ha) = ', tc/kb_HaK, ' (K)' 212 end if 213 214 ABI_FREE(lambda_1d) 215 ABI_FREE(z_1d) 216 ABI_FREE(delta_1d) 217 ABI_FREE(mm_1d) 218 219 close (UNIT=unit_z) 220 close (UNIT=unit_lam) 221 close (UNIT=unit_del) 222 223 write(std_out,*) ' eliashberg_1d : end ' 224 225 end subroutine eliashberg_1d