TABLE OF CONTENTS
- ABINIT/eli_app_m_1d
- ABINIT/eli_diag_m_1d
- ABINIT/eli_lambda_1d
- ABINIT/eli_m_iter_1d
- ABINIT/eli_z_1d
- ABINIT/eliashberg_1d
ABINIT/eli_app_m_1d [ Functions ]
NAME
eli_app_m_1d
FUNCTION
Apply the linearized Eliashberg matrix once to the input vector.
COPYRIGHT
Copyright (C) 2004-2018 ABINIT group (MVer) This file is distributed under the terms of the GNU General Public Licence, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
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
PARENTS
eliashberg_1d
CHILDREN
SOURCE
254 subroutine eli_app_m_1d (delta_1d,lambda_1d,nmatsu,z_1d) 255 256 use defs_basis 257 use defs_elphon 258 use m_profiling_abi 259 260 !This section has been created automatically by the script Abilint (TD). 261 !Do not modify the following lines by hand. 262 #undef ABI_FUNC 263 #define ABI_FUNC 'eli_app_m_1d' 264 !End of the abilint section 265 266 implicit none 267 268 !Arguments ------------------------------------ 269 !scalars 270 integer,intent(in) :: nmatsu 271 !arrays 272 real(dp),intent(in) :: lambda_1d(-nmatsu:nmatsu),z_1d(-nmatsu:nmatsu) 273 real(dp),intent(inout) :: delta_1d(-nmatsu:nmatsu) 274 275 !Local variables------------------------------- 276 !scalars 277 integer :: imatsu,jmatsu,miguelflag 278 real(dp) :: zfact 279 !arrays 280 real(dp) :: delta_tmp(-nmatsu:nmatsu),freqfact(-nmatsu:nmatsu) 281 282 ! ********************************************************************* 283 284 miguelflag = 0 285 286 287 do imatsu=-nmatsu,nmatsu 288 freqfact(imatsu) = one / abs(two*imatsu+one) 289 end do 290 291 delta_tmp(:) = delta_1d(:) 292 293 if (miguelflag == 1) then 294 do imatsu=-nmatsu,nmatsu 295 ! zfact = pi*tc / z_1d(imatsu) 296 zfact = one / z_1d(imatsu) 297 298 do jmatsu=max(-nmatsu,-nmatsu+imatsu),min(nmatsu,nmatsu+imatsu) 299 delta_tmp(imatsu) = delta_tmp(imatsu) & 300 & + delta_1d(jmatsu) & 301 & * lambda_1d(imatsu-jmatsu) & 302 & * freqfact(jmatsu) 303 end do 304 delta_tmp(imatsu) = delta_tmp(imatsu)*zfact 305 end do 306 307 else 308 309 ! i < 0 310 do imatsu=-nmatsu,-1 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 ! i > 0 328 do imatsu=0,nmatsu 329 330 ! j < 0 331 do jmatsu=max(-nmatsu,-nmatsu+imatsu),-1 332 delta_tmp(imatsu) = delta_tmp(imatsu) & 333 & + lambda_1d(imatsu-jmatsu)*delta_1d(jmatsu)*freqfact(jmatsu) & 334 & + lambda_1d(imatsu-jmatsu)*delta_1d(imatsu)*freqfact(imatsu) 335 end do 336 ! j > 0 337 do jmatsu=0,min(nmatsu,nmatsu+imatsu) 338 delta_tmp(imatsu) = delta_tmp(imatsu) & 339 & + lambda_1d(imatsu-jmatsu)*delta_1d(jmatsu)*freqfact(jmatsu) & 340 & - lambda_1d(imatsu-jmatsu)*delta_1d(imatsu)*freqfact(imatsu) 341 end do 342 343 end do 344 345 end if 346 347 delta_1d(:) = delta_tmp(:) 348 349 end subroutine eli_app_m_1d
ABINIT/eli_diag_m_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
COPYRIGHT
Copyright (C) 2004-2018 ABINIT group (MVer) This file is distributed under the terms of the GNU General Public Licence, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
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
PARENTS
eliashberg_1d
CHILDREN
NOTES
SOURCE
391 subroutine eli_diag_m_1d (delta_1d,lambda_1d,maxeigval,mustar,nmatsu,tc,z_1d) 392 393 use defs_basis 394 use defs_elphon 395 use m_profiling_abi 396 use m_linalg_interfaces 397 398 !This section has been created automatically by the script Abilint (TD). 399 !Do not modify the following lines by hand. 400 #undef ABI_FUNC 401 #define ABI_FUNC 'eli_diag_m_1d' 402 !End of the abilint section 403 404 implicit none 405 406 !Arguments ------------------------------------ 407 !scalars 408 integer,intent(in) :: nmatsu 409 real(dp),intent(in) :: mustar,tc 410 real(dp),intent(out) :: maxeigval 411 !arrays 412 real(dp),intent(in) :: lambda_1d(-nmatsu:nmatsu),z_1d(-nmatsu:nmatsu) 413 real(dp),intent(inout) :: delta_1d(-nmatsu:nmatsu) 414 415 !Local variables------------------------------- 416 !scalars 417 integer :: imatsu,info,jmatsu,kmatsu,lwork,tmiguel 418 real(dp) :: si,sj,sqtimat,sqtjmat 419 !arrays 420 real(dp) :: mtm_eig(2*nmatsu+1),symm_mtm(-nmatsu:nmatsu,-nmatsu:nmatsu) 421 real(dp) :: work(3*(2*nmatsu+1)) 422 423 ! ********************************************************************* 424 425 tmiguel = 0 426 427 if (tmiguel == 1) then 428 do imatsu=-nmatsu,nmatsu 429 do jmatsu=-nmatsu,nmatsu 430 431 symm_mtm(imatsu,jmatsu) = zero 432 do kmatsu=max(-nmatsu+imatsu,-nmatsu+jmatsu,-nmatsu),min(nmatsu+imatsu,nmatsu+jmatsu,nmatsu) 433 symm_mtm(imatsu,jmatsu) = symm_mtm(imatsu,jmatsu) & 434 & + lambda_1d(kmatsu-imatsu)*lambda_1d(kmatsu-jmatsu) & 435 & / ( z_1d(kmatsu)*z_1d(kmatsu) ) 436 end do 437 ! symm_mtm(imatsu,jmatsu) = symm_mtm(imatsu,jmatsu) / ((two*imatsu+one)*(two*jmatsu+one)) 438 symm_mtm(imatsu,jmatsu) = symm_mtm(imatsu,jmatsu) * pi * tc * pi * tc 439 end do 440 end do 441 442 else 443 444 symm_mtm(:,:) = -mustar 445 446 si = -one 447 do imatsu=-nmatsu,nmatsu 448 sqtimat = one / sqrt(two*abs(imatsu)+one) 449 if (imatsu == 0) si = one 450 sj = -one 451 do jmatsu=max(-nmatsu,-nmatsu+imatsu),min(nmatsu,nmatsu+imatsu) 452 sqtjmat = one / sqrt(two*abs(jmatsu)+one) 453 if (jmatsu == 0) sj = one 454 455 symm_mtm(imatsu,jmatsu) = symm_mtm(imatsu,jmatsu) & 456 & + lambda_1d(imatsu-jmatsu)*sqtimat*sqtjmat 457 458 symm_mtm(imatsu,imatsu) = symm_mtm(imatsu,imatsu) & 459 & - lambda_1d(imatsu-jmatsu)*si*sj*sqtimat*sqtimat 460 end do 461 end do 462 463 end if 464 465 lwork = 3*(2*nmatsu+1) 466 call DSYEV('V', 'U', 2*nmatsu+1, symm_mtm, 2*nmatsu+1, mtm_eig, work, lwork, info ) 467 468 write(std_out,*) 'last eigenvalues = ' 469 write(std_out,*) mtm_eig(2*nmatsu-9:2*nmatsu+1) 470 471 do imatsu=-nmatsu,nmatsu 472 delta_1d(imatsu) = symm_mtm(imatsu,nmatsu)*sqrt(two*abs(imatsu)+one) 473 end do 474 475 maxeigval = mtm_eig(2*nmatsu+1) 476 477 end subroutine eli_diag_m_1d
ABINIT/eli_lambda_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
COPYRIGHT
Copyright (C) 2004-2018 ABINIT group (MVer) This file is distributed under the terms of the GNU General Public Licence, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
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
SIDE EFFECTS
PARENTS
eliashberg_1d
CHILDREN
NOTES
lambda is used at points which are differences of Matsubara freqs, and hence is tabulated on points going through 0.
SOURCE
519 subroutine eli_lambda_1d (a2f_1d,elph_ds,lambda_1d,nmatsu,tc) 520 521 use defs_basis 522 use defs_elphon 523 use m_profiling_abi 524 525 use m_numeric_tools, only : simpson_int 526 527 !This section has been created automatically by the script Abilint (TD). 528 !Do not modify the following lines by hand. 529 #undef ABI_FUNC 530 #define ABI_FUNC 'eli_lambda_1d' 531 !End of the abilint section 532 533 implicit none 534 535 !Arguments ------------------------------------ 536 !scalars 537 integer,intent(in) :: nmatsu 538 real(dp),intent(in) :: tc 539 type(elph_type),intent(in) :: elph_ds 540 !arrays 541 real(dp),intent(in) :: a2f_1d(elph_ds%na2f) 542 real(dp),intent(out) :: lambda_1d(-nmatsu:nmatsu) 543 544 !Local variables------------------------------- 545 !scalars 546 integer :: imatsu,iomega 547 real(dp) :: nu_matsu,nu_matsu2,omega,domega 548 !arrays 549 real(dp) :: lambda_int(elph_ds%na2f),tmplambda(elph_ds%na2f) 550 551 ! ********************************************************************* 552 ! 553 !MG: the step should be calculated locally using nomega and the extrema of the spectrum. 554 !One should not rely on previous calls for the setup of elph_ds%domega 555 !I will remove elph_ds%domega since mka2f.F90 will become a method of gamma_t 556 domega =elph_ds%domega 557 558 do imatsu=-nmatsu,nmatsu 559 nu_matsu = (two*imatsu)*pi*tc 560 nu_matsu2 = nu_matsu*nu_matsu 561 562 tmplambda(:) = zero 563 omega=domega 564 do iomega=2,elph_ds%na2f 565 tmplambda(iomega) = a2f_1d(iomega) * two * omega / (nu_matsu2 + omega*omega) 566 omega=omega+domega 567 end do 568 call simpson_int(elph_ds%na2f,domega,tmplambda,lambda_int) 569 570 lambda_1d(imatsu) = lambda_int(elph_ds%na2f) 571 end do 572 573 end subroutine eli_lambda_1d
ABINIT/eli_m_iter_1d [ Functions ]
NAME
eli_m_iter_1d
FUNCTION
Find largest eigenvalue of M matrix, to deduce superconducting Tc
COPYRIGHT
Copyright (C) 2004-2018 ABINIT group (MVer) This file is distributed under the terms of the GNU General Public Licence, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
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
CHILDREN
eli_app_m_1d
NOTES
PARENTS
CHILDREN
SOURCE
614 subroutine eli_m_iter_1d (delta_1d,lambda_1d,maxeigval,nmatsu,z_1d) 615 616 use defs_basis 617 use defs_elphon 618 use m_profiling_abi 619 620 !This section has been created automatically by the script Abilint (TD). 621 !Do not modify the following lines by hand. 622 #undef ABI_FUNC 623 #define ABI_FUNC 'eli_m_iter_1d' 624 use interfaces_77_ddb, except_this_one => eli_m_iter_1d 625 !End of the abilint section 626 627 implicit none 628 629 !Arguments ------------------------------------ 630 !scalars 631 integer,intent(in) :: nmatsu 632 real(dp),intent(out) :: maxeigval 633 !arrays 634 real(dp),intent(in) :: lambda_1d(-nmatsu:nmatsu),z_1d(-nmatsu:nmatsu) 635 real(dp),intent(inout) :: delta_1d(-nmatsu:nmatsu) 636 637 !Local variables------------------------------- 638 !scalars 639 integer :: iiterm,imatsu,nfilter,ngeteig 640 real(dp) :: dnewnorm,dnorm,fact 641 !arrays 642 real(dp) :: delta_old(-nmatsu:nmatsu) 643 644 ! ********************************************************************* 645 646 nfilter = 10 647 ngeteig = 10 648 649 650 ! 651 !1) apply M matrix enough times to filter out largest eigenvalue 652 ! 653 do iiterm=1,nfilter 654 call eli_app_m_1d (delta_1d,lambda_1d,nmatsu,z_1d) 655 656 ! DEBUG 657 ! dnorm=zero 658 ! do imatsu=-nmatsu,nmatsu 659 ! dnorm = dnorm + delta_1d(imatsu)*delta_1d(imatsu)/(two*imatsu+one) 660 ! end do 661 ! dnorm = sqrt(dnorm) 662 ! write(std_out,*) 'eli_m_iter_1d : dnorm ', dnorm 663 ! ENDDEBUG 664 end do 665 666 ! 667 !2) calculate norm 668 ! 669 dnorm=zero 670 do imatsu=-nmatsu,nmatsu 671 dnorm = dnorm + delta_1d(imatsu)*delta_1d(imatsu)/abs(two*imatsu+one) 672 end do 673 dnorm = sqrt(dnorm) 674 675 !normalize delta_1d 676 delta_1d(:) = delta_1d(:) / dnorm 677 678 delta_old = delta_1d 679 680 !DEBUG 681 !dnewnorm=zero 682 !do imatsu=-nmatsu,nmatsu 683 !dnewnorm = dnewnorm + delta_1d(imatsu)*delta_1d(imatsu)/abs(two*imatsu+one) 684 !end do 685 !dnewnorm = sqrt(dnewnorm) 686 !write(std_out,*) 'eli_m_iter_1d : dnewnorm1 ', dnewnorm 687 !ENDDEBUG 688 689 ! 690 !3) re-apply M matrix ngeteig times 691 ! 692 do iiterm=1,ngeteig 693 call eli_app_m_1d (delta_1d,lambda_1d,nmatsu,z_1d) 694 ! DEBUG 695 ! dnewnorm=zero 696 ! do imatsu=-nmatsu,nmatsu 697 ! dnewnorm = dnewnorm + delta_1d(imatsu)*delta_1d(imatsu)/abs(two*imatsu+one) 698 ! end do 699 ! dnewnorm = sqrt(dnewnorm) 700 ! write(std_out,*) 'eli_m_iter_1d : dnewnorm ', dnewnorm 701 ! 702 ! do imatsu=-nmatsu,nmatsu 703 ! write (112,*) imatsu,delta_1d(imatsu)/delta_old(imatsu) 704 ! end do 705 ! write (112,*) 706 ! delta_old = delta_1d 707 ! ENDDEBUG 708 end do 709 710 ! 711 !4) calculate new norm and estimate eigenvalue 712 ! 713 dnewnorm=zero 714 do imatsu=-nmatsu,nmatsu 715 dnewnorm = dnewnorm + delta_1d(imatsu)*delta_1d(imatsu)/abs(two*imatsu+one) 716 end do 717 dnewnorm = sqrt(dnewnorm) 718 719 !maxeigval = exp ( log(dnewnorm/dnorm) / ngeteig ) 720 maxeigval = exp ( log(dnewnorm) / ngeteig ) 721 722 write(std_out,*) 'eli_m_iter_1d : maxeigval =' , maxeigval 723 !fact = exp(-log(maxeigval) * (ngeteig+nfilter)) 724 fact = exp(-log(maxeigval) * (ngeteig)) 725 do imatsu=-nmatsu,nmatsu 726 delta_1d(imatsu) = delta_1d(imatsu) * fact 727 end do 728 729 end subroutine eli_m_iter_1d
ABINIT/eli_z_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
COPYRIGHT
Copyright (C) 2004-2018 ABINIT group (MVer) This file is distributed under the terms of the GNU General Public Licence, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
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
PARENTS
eliashberg_1d
CHILDREN
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
768 subroutine eli_z_1d (lambda_1d,nmatsu,z_1d) 769 770 use defs_basis 771 use defs_elphon 772 use m_profiling_abi 773 774 !This section has been created automatically by the script Abilint (TD). 775 !Do not modify the following lines by hand. 776 #undef ABI_FUNC 777 #define ABI_FUNC 'eli_z_1d' 778 !End of the abilint section 779 780 implicit none 781 782 !Arguments ------------------------------------ 783 !scalars 784 integer,intent(in) :: nmatsu 785 !arrays 786 real(dp),intent(in) :: lambda_1d(-nmatsu:nmatsu) 787 real(dp),intent(out) :: z_1d(-nmatsu:nmatsu) 788 789 !Local variables------------------------------- 790 !scalars 791 integer :: imatsu,jmatsu 792 793 ! ********************************************************************* 794 795 796 do imatsu=0,nmatsu 797 798 z_1d(imatsu) = zero 799 ! count $\mathrm{sign}(omega_{Matsubara})$ 800 do jmatsu=-nmatsu+imatsu,-1 801 z_1d(imatsu) = z_1d(imatsu) - lambda_1d(imatsu-jmatsu) 802 end do 803 do jmatsu=0,nmatsu 804 z_1d(imatsu) = z_1d(imatsu) + lambda_1d(imatsu-jmatsu) 805 end do 806 807 ! NOTE: the pi*Tc factor in Z cancels the one in the Matsubara frequency. 808 z_1d(imatsu) = one + z_1d(imatsu) / (two*imatsu+one) 809 z_1d(-imatsu) = z_1d(imatsu) 810 end do 811 812 end subroutine eli_z_1d
ABINIT/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.
COPYRIGHT
Copyright (C) 2004-2018 ABINIT group (MVer) This file is distributed under the terms of the GNU General Public Licence, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
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
PARENTS
elphon
CHILDREN
NOTES
na2f = number of frequency points for alpha^2F function
SOURCE
37 #if defined HAVE_CONFIG_H 38 #include "config.h" 39 #endif 40 41 #include "abi_common.h" 42 43 44 subroutine eliashberg_1d(a2f_1d,elph_ds,mustar) 45 46 use defs_basis 47 use defs_elphon 48 use m_io_tools 49 use m_profiling_abi 50 51 !This section has been created automatically by the script Abilint (TD). 52 !Do not modify the following lines by hand. 53 #undef ABI_FUNC 54 #define ABI_FUNC 'eliashberg_1d' 55 use interfaces_14_hidewrite 56 use interfaces_77_ddb, except_this_one => eliashberg_1d 57 !End of the abilint section 58 59 implicit none 60 61 !Arguments ------------------------------------ 62 ! needed for phonon interpolation 63 !scalars 64 real(dp),intent(in) :: mustar 65 type(elph_type),intent(in) :: elph_ds 66 !arrays 67 real(dp),intent(in) :: a2f_1d(elph_ds%na2f) 68 69 !Local variables------------------------------- 70 ! for diagonalization of gammma matrix 71 ! output variables for gtdyn9+dfpt_phfrq 72 !scalars 73 integer :: iiter,imatsu 74 integer :: maxiter,nmatsu,unit_del,unit_lam,unit_z 75 real(dp) :: maxeigval,omega_cutoff 76 real(dp) :: tc 77 character(len=fnlen) :: fname 78 !arrays 79 real(dp),allocatable :: delta_1d(:),lambda_1d(:),mm_1d(:,:),z_1d(:) 80 81 ! ********************************************************************* 82 83 call wrtout(std_out,'Solving the 1-D Eliashberg equation (isotropic case)',"COLL") 84 85 if (elph_ds%nsppol /= 1) then 86 write(std_out,*) 'eliashberg_1d is not coded for nsppol > 1 yet' 87 return 88 end if 89 90 !maximum number of iterations to find T_c 91 maxiter=30 92 93 !Fix nmatsu. Should add test in iiter loop to check if 94 !omega_cutoff is respected 95 nmatsu = 50 96 !write(std_out,*) ' eliashberg_1d : nmatsu = ', nmatsu 97 98 ABI_ALLOCATE(lambda_1d,(-nmatsu:nmatsu)) 99 ABI_ALLOCATE(z_1d,(-nmatsu:nmatsu)) 100 ABI_ALLOCATE(delta_1d,(-nmatsu:nmatsu)) 101 ABI_ALLOCATE(mm_1d,(-nmatsu:nmatsu,-nmatsu:nmatsu)) 102 103 unit_lam=get_unit() 104 fname=trim(elph_ds%elph_base_name) // "_LAM" 105 open (UNIT=unit_lam,FILE=fname,STATUS='REPLACE') 106 unit_z=get_unit() 107 fname=trim(elph_ds%elph_base_name) // "_Z" 108 open (UNIT=unit_z,FILE=fname,STATUS='REPLACE') 109 unit_del=get_unit() 110 fname=trim(elph_ds%elph_base_name) // "_DEL" 111 open (UNIT=unit_del,FILE=fname,STATUS='REPLACE') 112 113 ! 114 !1) use linearized Eliashberg equation to find Tc 115 !! \sum_j \mathbf{M}_{ij} \Delta_j = \zeta \cdot \Delta_i $ $i,j = 1 .. n_{\mathrm{Matsubara}}$ 116 !$\zeta = 1$ gives T$_c$ $\beta = \frac{1}{\mathrm{T}}$ $\omega_i = (2 i + 1) \pi \mathrm{T}$ 117 !! \mathbf{M}_{ij} = \frac{\pi}{\beta} \frac{\lambda (\omega_i - \omega_j)}{Z (\omega_i)}$ 118 !! Z (\omega_i) = 1 + \frac{\pi}{\beta \omega_i} \sum_j \lambda(\omega_i - \omega_j) \mathrm{sgn}(\omega_j)$ 119 ! 120 121 !initial guess for T$_c$ in Hartree (1Ha =3.067e5 K) 122 tc = 0.0001 123 ! 124 !big iterative loop 125 ! 126 do iiter=1,maxiter 127 128 omega_cutoff = (two*nmatsu+one) * pi * tc 129 130 ! 131 ! calculate array of lambda values 132 ! 133 call eli_lambda_1d (a2f_1d,elph_ds,lambda_1d,nmatsu,tc) 134 write (unit_lam,'(a)') '#' 135 write (unit_lam,'(a)') '# ABINIT package : lambda file' 136 write (unit_lam,'(a)') '#' 137 write (unit_lam,'(a,I10,a)') '# lambda_1d array containing 2*', nmatsu, '+1 Matsubara frequency points' 138 write (unit_lam,'(a,E16.6,a,E16.6)') '# from ', -omega_cutoff, ' to ', omega_cutoff 139 write (unit_lam,'(a)') '# lambda_1d is the frequency dependent coupling constant ' 140 write (unit_lam,'(a)') '# in the Eliashberg equations ' 141 write (unit_lam,'(a)') '#' 142 do imatsu=-nmatsu,nmatsu 143 write (unit_lam,*) imatsu,lambda_1d(imatsu) 144 end do 145 write (unit_lam,*) 146 147 ! 148 ! calculate array of z values 149 ! 150 call eli_z_1d (lambda_1d,nmatsu,z_1d) 151 write (unit_z,'(a)') '#' 152 write (unit_z,'(a)') '# ABINIT package : Z file' 153 write (unit_z,'(a)') '#' 154 write (unit_z,'(a,I10,a)') '# z_1d array containing 2*', nmatsu, '+1 Matsubara frequency points' 155 write (unit_z,'(a,E16.6,a,E16.6)') '# from ', -omega_cutoff, ' to ', omega_cutoff 156 write (unit_z,'(a)') '# z_1d is the renormalization factor in the Eliashberg equations' 157 write (unit_z,'(a)') '#' 158 do imatsu=-nmatsu,nmatsu 159 write (unit_z,*) imatsu,z_1d(imatsu) 160 end do 161 162 ! ! 163 ! ! apply M matrix until a maximal eigenvalue is found. 164 ! ! 165 ! call eli_m_iter_1d (delta_1d,lambda_1d,maxeigval,nmatsu,z_1d) 166 167 ! 168 ! diagonalize M brute forcefully 169 ! 170 call eli_diag_m_1d(delta_1d,lambda_1d,maxeigval,mustar,nmatsu,tc,z_1d) 171 172 write (unit_del,'(a)') '#' 173 write (unit_del,'(a)') '# eliashberg_1d : delta_1d = ' 174 write (unit_del,'(a)') '#' 175 write (unit_del,'(a,i6,a)') '# delta_1d array containing 2*', nmatsu, '+1 Matsubara frequency points' 176 write (unit_z,'(a,E16.6,a,E16.6)') '# from ', -omega_cutoff, ' to ', omega_cutoff 177 write (unit_z,'(a)') '# delta_1d is the gap function in the Eliashberg equations' 178 write (unit_z,'(a)') '#' 179 do imatsu=-nmatsu,nmatsu 180 write (unit_del,*) imatsu,delta_1d(imatsu) 181 end do 182 write (unit_del,*) 183 184 ! if eigenvalue is < 1 increase T 185 ! else if eigenvalue is > 1 decrease T 186 ! if eigenvalue ~= 1 stop 187 ! 188 if (abs(maxeigval-one) < tol8) then 189 write(std_out,*) 'Eliashberg Tc found = ', tc, ' (Ha) = ', tc/kb_HaK, ' (K)' 190 exit 191 else if (maxeigval > 0.001_dp) then 192 tc = tc * maxeigval 193 else 194 write(std_out,*) 'maxeigval is very small' 195 tc = tc * 1000.0_dp 196 end if 197 198 199 end do 200 !end iiter do 201 202 if (abs(maxeigval-one) > tol8) then 203 write(std_out,*) 'eliashberg_1d : Tc not converged. ', maxeigval, ' /= 1' 204 write(std_out,*) 'Eliashberg Tc nonetheless = ', tc, ' (Ha) = ', tc/kb_HaK, ' (K)' 205 end if 206 207 ABI_DEALLOCATE(lambda_1d) 208 ABI_DEALLOCATE(z_1d) 209 ABI_DEALLOCATE(delta_1d) 210 ABI_DEALLOCATE(mm_1d) 211 212 close (UNIT=unit_z) 213 close (UNIT=unit_lam) 214 close (UNIT=unit_del) 215 216 write(std_out,*) ' eliashberg_1d : end ' 217 218 219 end subroutine eliashberg_1d