## ABINIT/m_eliashberg_1d [ Modules ]

[ Top ] [ Modules ]

NAME

  m_eliashberg_1d


FUNCTION

  Solve the Eliashberg equations in the isotropic case


  Copyright (C) 2008-2018 ABINIT group (MVer)
GNU General Public License, see ~abinit/COPYING
or http://www.gnu.org/copyleft/gpl.txt .


PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24
25 #include "abi_common.h"
26
27 module m_eliashberg_1d
28
29  use defs_basis
30  use defs_elphon
31  use m_errors
32  use m_abicore
33  use m_io_tools
34  use m_abicore
35
36  use m_numeric_tools,   only : simpson_int
37
38  implicit none
39
40  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


PARENTS

      eliashberg_1d


CHILDREN

SOURCE

272 subroutine eli_app_m_1d (delta_1d,lambda_1d,nmatsu,z_1d)
273
274
275 !This section has been created automatically by the script Abilint (TD).
276 !Do not modify the following lines by hand.
277 #undef ABI_FUNC
278 #define ABI_FUNC 'eli_app_m_1d'
279 !End of the abilint section
280
281  implicit none
282
283 !Arguments ------------------------------------
284 !scalars
285  integer,intent(in) :: nmatsu
286 !arrays
287  real(dp),intent(in) :: lambda_1d(-nmatsu:nmatsu),z_1d(-nmatsu:nmatsu)
288  real(dp),intent(inout) :: delta_1d(-nmatsu:nmatsu)
289
290 !Local variables-------------------------------
291 !scalars
292  integer :: imatsu,jmatsu,miguelflag
293  real(dp) :: zfact
294 !arrays
295  real(dp) :: delta_tmp(-nmatsu:nmatsu),freqfact(-nmatsu:nmatsu)
296
297 ! *********************************************************************
298
299  miguelflag = 0
300
301
302  do imatsu=-nmatsu,nmatsu
303    freqfact(imatsu) = one / abs(two*imatsu+one)
304  end do
305
306  delta_tmp(:) = delta_1d(:)
307
308  if (miguelflag == 1) then
309    do imatsu=-nmatsu,nmatsu
310 !    zfact = pi*tc / z_1d(imatsu)
311      zfact = one / z_1d(imatsu)
312
313      do jmatsu=max(-nmatsu,-nmatsu+imatsu),min(nmatsu,nmatsu+imatsu)
314        delta_tmp(imatsu) = delta_tmp(imatsu) &
315 &       + delta_1d(jmatsu) &
316 &       * lambda_1d(imatsu-jmatsu) &
317 &       * freqfact(jmatsu)
318      end do
319      delta_tmp(imatsu) = delta_tmp(imatsu)*zfact
320    end do
321
322  else
323
324 !  i < 0
325    do imatsu=-nmatsu,-1
326
327 !    j < 0
328      do jmatsu=max(-nmatsu,-nmatsu+imatsu),-1
329        delta_tmp(imatsu) = delta_tmp(imatsu) &
330 &       + lambda_1d(imatsu-jmatsu)*delta_1d(jmatsu)*freqfact(jmatsu) &
331 &       - lambda_1d(imatsu-jmatsu)*delta_1d(imatsu)*freqfact(imatsu)
332      end do
333 !    j > 0
334      do jmatsu=0,min(nmatsu,nmatsu+imatsu)
335        delta_tmp(imatsu) = delta_tmp(imatsu) &
336 &       + lambda_1d(imatsu-jmatsu)*delta_1d(jmatsu)*freqfact(jmatsu) &
337 &       + lambda_1d(imatsu-jmatsu)*delta_1d(imatsu)*freqfact(imatsu)
338      end do
339
340    end do
341
342 !  i > 0
343    do imatsu=0,nmatsu
344
345 !    j < 0
346      do jmatsu=max(-nmatsu,-nmatsu+imatsu),-1
347        delta_tmp(imatsu) = delta_tmp(imatsu) &
348 &       + lambda_1d(imatsu-jmatsu)*delta_1d(jmatsu)*freqfact(jmatsu) &
349 &       + lambda_1d(imatsu-jmatsu)*delta_1d(imatsu)*freqfact(imatsu)
350      end do
351 !    j > 0
352      do jmatsu=0,min(nmatsu,nmatsu+imatsu)
353        delta_tmp(imatsu) = delta_tmp(imatsu) &
354 &       + lambda_1d(imatsu-jmatsu)*delta_1d(jmatsu)*freqfact(jmatsu) &
355 &       - lambda_1d(imatsu-jmatsu)*delta_1d(imatsu)*freqfact(imatsu)
356      end do
357
358    end do
359
360  end if
361
362  delta_1d(:) = delta_tmp(:)
363
364 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


PARENTS

      eliashberg_1d


CHILDREN

SOURCE

397 subroutine eli_diag_m_1d (delta_1d,lambda_1d,maxeigval,mustar,nmatsu,tc,z_1d)
398
399  use m_linalg_interfaces
400
401 !This section has been created automatically by the script Abilint (TD).
402 !Do not modify the following lines by hand.
403 #undef ABI_FUNC
404 #define ABI_FUNC 'eli_diag_m_1d'
405 !End of the abilint section
406
407  implicit none
408
409 !Arguments ------------------------------------
410 !scalars
411  integer,intent(in) :: nmatsu
412  real(dp),intent(in) :: mustar,tc
413  real(dp),intent(out) :: maxeigval
414 !arrays
415  real(dp),intent(in) :: lambda_1d(-nmatsu:nmatsu),z_1d(-nmatsu:nmatsu)
416  real(dp),intent(inout) :: delta_1d(-nmatsu:nmatsu)
417
418 !Local variables-------------------------------
419 !scalars
420  integer :: imatsu,info,jmatsu,kmatsu,lwork,tmiguel
421  real(dp) :: si,sj,sqtimat,sqtjmat
422 !arrays
423  real(dp) :: mtm_eig(2*nmatsu+1),symm_mtm(-nmatsu:nmatsu,-nmatsu:nmatsu)
424  real(dp) :: work(3*(2*nmatsu+1))
425
426 ! *********************************************************************
427
428  tmiguel = 0
429
430  if (tmiguel == 1) then
431    do imatsu=-nmatsu,nmatsu
432      do jmatsu=-nmatsu,nmatsu
433
434        symm_mtm(imatsu,jmatsu) = zero
435        do kmatsu=max(-nmatsu+imatsu,-nmatsu+jmatsu,-nmatsu),min(nmatsu+imatsu,nmatsu+jmatsu,nmatsu)
436          symm_mtm(imatsu,jmatsu) = symm_mtm(imatsu,jmatsu) &
437 &         +  lambda_1d(kmatsu-imatsu)*lambda_1d(kmatsu-jmatsu) &
438 &         /  ( z_1d(kmatsu)*z_1d(kmatsu) )
439        end do
440 !      symm_mtm(imatsu,jmatsu) = symm_mtm(imatsu,jmatsu) / ((two*imatsu+one)*(two*jmatsu+one))
441        symm_mtm(imatsu,jmatsu) = symm_mtm(imatsu,jmatsu) * pi * tc * pi * tc
442      end do
443    end do
444
445  else
446
447    symm_mtm(:,:) = -mustar
448
449    si = -one
450    do imatsu=-nmatsu,nmatsu
451      sqtimat = one / sqrt(two*abs(imatsu)+one)
452      if (imatsu == 0) si = one
453      sj = -one
454      do jmatsu=max(-nmatsu,-nmatsu+imatsu),min(nmatsu,nmatsu+imatsu)
455        sqtjmat = one / sqrt(two*abs(jmatsu)+one)
456        if (jmatsu == 0) sj = one
457
458        symm_mtm(imatsu,jmatsu) = symm_mtm(imatsu,jmatsu) &
459 &       + lambda_1d(imatsu-jmatsu)*sqtimat*sqtjmat
460
461        symm_mtm(imatsu,imatsu) = symm_mtm(imatsu,imatsu) &
462 &       - lambda_1d(imatsu-jmatsu)*si*sj*sqtimat*sqtimat
463      end do
464    end do
465
466  end if
467
468  lwork = 3*(2*nmatsu+1)
469  call DSYEV('V', 'U', 2*nmatsu+1, symm_mtm, 2*nmatsu+1, mtm_eig, work, lwork, info )
470
471  write(std_out,*) 'last eigenvalues = '
472  write(std_out,*) mtm_eig(2*nmatsu-9:2*nmatsu+1)
473
474  do imatsu=-nmatsu,nmatsu
475    delta_1d(imatsu) = symm_mtm(imatsu,nmatsu)*sqrt(two*abs(imatsu)+one)
476  end do
477
478  maxeigval = mtm_eig(2*nmatsu+1)
479
480 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


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

513 subroutine eli_lambda_1d (a2f_1d,elph_ds,lambda_1d,nmatsu,tc)
514
515
516 !This section has been created automatically by the script Abilint (TD).
517 !Do not modify the following lines by hand.
518 #undef ABI_FUNC
519 #define ABI_FUNC 'eli_lambda_1d'
520 !End of the abilint section
521
522  implicit none
523
524 !Arguments ------------------------------------
525 !scalars
526  integer,intent(in) :: nmatsu
527  real(dp),intent(in) :: tc
528  type(elph_type),intent(in) :: elph_ds
529 !arrays
530  real(dp),intent(in) :: a2f_1d(elph_ds%na2f)
531  real(dp),intent(out) :: lambda_1d(-nmatsu:nmatsu)
532
533 !Local variables-------------------------------
534 !scalars
535  integer :: imatsu,iomega
536  real(dp) :: nu_matsu,nu_matsu2,omega,domega
537 !arrays
538  real(dp) :: lambda_int(elph_ds%na2f),tmplambda(elph_ds%na2f)
539
540 ! *********************************************************************
541 !
542 !MG: the step should be calculated locally using nomega and the extrema of the spectrum.
543 !One should not rely on previous calls for the setup of elph_ds%domega
544 !I will remove elph_ds%domega since mka2f.F90 will become a method of gamma_t
545  domega =elph_ds%domega
546
547  do imatsu=-nmatsu,nmatsu
548    nu_matsu = (two*imatsu)*pi*tc
549    nu_matsu2 = nu_matsu*nu_matsu
550
551    tmplambda(:) = zero
552    omega=domega
553    do iomega=2,elph_ds%na2f
554      tmplambda(iomega) = a2f_1d(iomega) * two * omega / (nu_matsu2 + omega*omega)
555      omega=omega+domega
556    end do
557    call simpson_int(elph_ds%na2f,domega,tmplambda,lambda_int)
558
559    lambda_1d(imatsu) = lambda_int(elph_ds%na2f)
560  end do
561
562 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


PARENTS

CHILDREN

      eli_app_m_1d


SOURCE

591 subroutine eli_m_iter_1d (delta_1d,lambda_1d,maxeigval,nmatsu,z_1d)
592
593
594 !This section has been created automatically by the script Abilint (TD).
595 !Do not modify the following lines by hand.
596 #undef ABI_FUNC
597 #define ABI_FUNC 'eli_m_iter_1d'
598 !End of the abilint section
599
600  implicit none
601
602 !Arguments ------------------------------------
603 !scalars
604  integer,intent(in) :: nmatsu
605  real(dp),intent(out) :: maxeigval
606 !arrays
607  real(dp),intent(in) :: lambda_1d(-nmatsu:nmatsu),z_1d(-nmatsu:nmatsu)
608  real(dp),intent(inout) :: delta_1d(-nmatsu:nmatsu)
609
610 !Local variables-------------------------------
611 !scalars
612  integer :: iiterm,imatsu,nfilter,ngeteig
613  real(dp) :: dnewnorm,dnorm,fact
614 !arrays
615  real(dp) :: delta_old(-nmatsu:nmatsu)
616
617 ! *********************************************************************
618
619  nfilter = 10
620  ngeteig = 10
621
622
623 !
624 !1) apply M matrix enough times to filter out largest eigenvalue
625 !
626  do iiterm=1,nfilter
627    call eli_app_m_1d (delta_1d,lambda_1d,nmatsu,z_1d)
628
629 !  DEBUG
630 !  dnorm=zero
631 !  do imatsu=-nmatsu,nmatsu
632 !  dnorm = dnorm + delta_1d(imatsu)*delta_1d(imatsu)/(two*imatsu+one)
633 !  end do
634 !  dnorm = sqrt(dnorm)
635 !  write(std_out,*) 'eli_m_iter_1d : dnorm ', dnorm
636 !  ENDDEBUG
637  end do
638
639 !
640 !2) calculate norm
641 !
642  dnorm=zero
643  do imatsu=-nmatsu,nmatsu
644    dnorm = dnorm + delta_1d(imatsu)*delta_1d(imatsu)/abs(two*imatsu+one)
645  end do
646  dnorm = sqrt(dnorm)
647
648 !normalize delta_1d
649  delta_1d(:) = delta_1d(:) / dnorm
650
651  delta_old = delta_1d
652
653 !DEBUG
654 !dnewnorm=zero
655 !do imatsu=-nmatsu,nmatsu
656 !dnewnorm = dnewnorm + delta_1d(imatsu)*delta_1d(imatsu)/abs(two*imatsu+one)
657 !end do
658 !dnewnorm = sqrt(dnewnorm)
659 !write(std_out,*) 'eli_m_iter_1d : dnewnorm1 ', dnewnorm
660 !ENDDEBUG
661
662 !
663 !3) re-apply M matrix ngeteig times
664 !
665  do iiterm=1,ngeteig
666    call eli_app_m_1d (delta_1d,lambda_1d,nmatsu,z_1d)
667 !  DEBUG
668 !  dnewnorm=zero
669 !  do imatsu=-nmatsu,nmatsu
670 !  dnewnorm = dnewnorm + delta_1d(imatsu)*delta_1d(imatsu)/abs(two*imatsu+one)
671 !  end do
672 !  dnewnorm = sqrt(dnewnorm)
673 !  write(std_out,*) 'eli_m_iter_1d : dnewnorm ', dnewnorm
674 !
675 !  do imatsu=-nmatsu,nmatsu
676 !  write (112,*) imatsu,delta_1d(imatsu)/delta_old(imatsu)
677 !  end do
678 !  write (112,*)
679 !  delta_old = delta_1d
680 !  ENDDEBUG
681  end do
682
683 !
684 !4) calculate new norm and estimate eigenvalue
685 !
686  dnewnorm=zero
687  do imatsu=-nmatsu,nmatsu
688    dnewnorm = dnewnorm + delta_1d(imatsu)*delta_1d(imatsu)/abs(two*imatsu+one)
689  end do
690  dnewnorm = sqrt(dnewnorm)
691
692 !maxeigval = exp ( log(dnewnorm/dnorm) / ngeteig )
693  maxeigval = exp ( log(dnewnorm) / ngeteig )
694
695  write(std_out,*) 'eli_m_iter_1d : maxeigval =' , maxeigval
696 !fact = exp(-log(maxeigval) * (ngeteig+nfilter))
697  fact = exp(-log(maxeigval) * (ngeteig))
698  do imatsu=-nmatsu,nmatsu
699    delta_1d(imatsu) = delta_1d(imatsu) * fact
700  end do
701
702 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


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

734 subroutine eli_z_1d (lambda_1d,nmatsu,z_1d)
735
736
737 !This section has been created automatically by the script Abilint (TD).
738 !Do not modify the following lines by hand.
739 #undef ABI_FUNC
740 #define ABI_FUNC 'eli_z_1d'
741 !End of the abilint section
742
743  implicit none
744
745 !Arguments ------------------------------------
746 !scalars
747  integer,intent(in) :: nmatsu
748 !arrays
749  real(dp),intent(in) :: lambda_1d(-nmatsu:nmatsu)
750  real(dp),intent(out) :: z_1d(-nmatsu:nmatsu)
751
752 !Local variables-------------------------------
753 !scalars
754  integer :: imatsu,jmatsu
755
756 ! *********************************************************************
757
758
759  do imatsu=0,nmatsu
760
761    z_1d(imatsu) = zero
762 !  count $\mathrm{sign}(omega_{Matsubara})$
763    do jmatsu=-nmatsu+imatsu,-1
764      z_1d(imatsu) = z_1d(imatsu) - lambda_1d(imatsu-jmatsu)
765    end do
766    do jmatsu=0,nmatsu
767      z_1d(imatsu) = z_1d(imatsu) + lambda_1d(imatsu-jmatsu)
768    end do
769
770 !  NOTE: the pi*Tc factor in Z cancels the one in the Matsubara frequency.
771    z_1d(imatsu) = one + z_1d(imatsu) / (two*imatsu+one)
772    z_1d(-imatsu) = z_1d(imatsu)
773  end do
774
775 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

PARENTS

      elphon


CHILDREN

NOTES

  na2f = number of frequency points for alpha^2F function


SOURCE

 77 subroutine eliashberg_1d(a2f_1d,elph_ds,mustar)
78
79
80 !This section has been created automatically by the script Abilint (TD).
81 !Do not modify the following lines by hand.
82 #undef ABI_FUNC
83 #define ABI_FUNC 'eliashberg_1d'
84 !End of the abilint section
85
86  implicit none
87
88 !Arguments ------------------------------------
89 !scalars
90  real(dp),intent(in) :: mustar
91  type(elph_type),intent(in) :: elph_ds
92 !arrays
93  real(dp),intent(in) :: a2f_1d(elph_ds%na2f)
94
95 !Local variables-------------------------------
96   ! for diagonalization of gammma matrix
97   ! output variables for gtdyn9+dfpt_phfrq
98 !scalars
99  integer :: iiter,imatsu
100  integer :: maxiter,nmatsu,unit_del,unit_lam,unit_z
101  real(dp) :: maxeigval,omega_cutoff
102  real(dp) :: tc
103  character(len=fnlen) :: fname
104 !arrays
105  real(dp),allocatable :: delta_1d(:),lambda_1d(:),mm_1d(:,:),z_1d(:)
106
107 ! *********************************************************************
108
109  call wrtout(std_out,'Solving the 1-D Eliashberg equation (isotropic case)',"COLL")
110
111  if (elph_ds%nsppol /= 1) then
112    write(std_out,*) 'eliashberg_1d is not coded for nsppol > 1 yet'
113    return
114  end if
115
116 !maximum number of iterations to find T_c
117  maxiter=30
118
119 !Fix nmatsu. Should add test in iiter loop to check if
120 !omega_cutoff is respected
121  nmatsu = 50
122 !write(std_out,*) ' eliashberg_1d : nmatsu = ', nmatsu
123
124  ABI_ALLOCATE(lambda_1d,(-nmatsu:nmatsu))
125  ABI_ALLOCATE(z_1d,(-nmatsu:nmatsu))
126  ABI_ALLOCATE(delta_1d,(-nmatsu:nmatsu))
127  ABI_ALLOCATE(mm_1d,(-nmatsu:nmatsu,-nmatsu:nmatsu))
128
129  unit_lam=get_unit()
130  fname=trim(elph_ds%elph_base_name) // "_LAM"
131  open (UNIT=unit_lam,FILE=fname,STATUS='REPLACE')
132  unit_z=get_unit()
133  fname=trim(elph_ds%elph_base_name) // "_Z"
134  open (UNIT=unit_z,FILE=fname,STATUS='REPLACE')
135  unit_del=get_unit()
136  fname=trim(elph_ds%elph_base_name) // "_DEL"
137  open (UNIT=unit_del,FILE=fname,STATUS='REPLACE')
138
139 !
140 !1) use linearized Eliashberg equation to find Tc
141 !! \sum_j \mathbf{M}_{ij} \Delta_j = \zeta \cdot \Delta_i i,j = 1 .. n_{\mathrm{Matsubara}}$142 !$\zeta = 1$gives T$_c\beta = \frac{1}{\mathrm{T}}\omega_i = (2 i + 1) \pi \mathrm{T}$143 !! \mathbf{M}_{ij} = \frac{\pi}{\beta} \frac{\lambda (\omega_i - \omega_j)}{Z (\omega_i)}$
144 !! Z (\omega_i) = 1 + \frac{\pi}{\beta \omega_i} \sum_j \lambda(\omega_i - \omega_j) \mathrm{sgn}(\omega_j)$145 ! 146 147 !initial guess for T$_c\$ in Hartree (1Ha =3.067e5 K)
148  tc = 0.0001
149 !
150 !big iterative loop
151 !
152  do iiter=1,maxiter
153
154    omega_cutoff = (two*nmatsu+one) * pi * tc
155
156 !
157 !  calculate array of lambda values
158 !
159    call eli_lambda_1d (a2f_1d,elph_ds,lambda_1d,nmatsu,tc)
160    write (unit_lam,'(a)') '#'
161    write (unit_lam,'(a)') '# ABINIT package : lambda file'
162    write (unit_lam,'(a)') '#'
163    write (unit_lam,'(a,I10,a)') '# lambda_1d array containing 2*', nmatsu, '+1 Matsubara frequency points'
164    write (unit_lam,'(a,E16.6,a,E16.6)') '#  from ', -omega_cutoff, ' to ', omega_cutoff
165    write (unit_lam,'(a)') '#  lambda_1d is the frequency dependent coupling constant '
166    write (unit_lam,'(a)') '#  in the Eliashberg equations '
167    write (unit_lam,'(a)') '#'
168    do imatsu=-nmatsu,nmatsu
169      write (unit_lam,*) imatsu,lambda_1d(imatsu)
170    end do
171    write (unit_lam,*)
172
173 !
174 !  calculate array of z values
175 !
176    call eli_z_1d (lambda_1d,nmatsu,z_1d)
177    write (unit_z,'(a)') '#'
178    write (unit_z,'(a)') '# ABINIT package : Z file'
179    write (unit_z,'(a)') '#'
180    write (unit_z,'(a,I10,a)') '# z_1d array containing 2*', nmatsu, '+1 Matsubara frequency points'
181    write (unit_z,'(a,E16.6,a,E16.6)') '# from ', -omega_cutoff, ' to ', omega_cutoff
182    write (unit_z,'(a)') '# z_1d is the renormalization factor in the Eliashberg equations'
183    write (unit_z,'(a)') '#'
184    do imatsu=-nmatsu,nmatsu
185      write (unit_z,*) imatsu,z_1d(imatsu)
186    end do
187
188 !  !
189 !  ! apply M matrix until a maximal eigenvalue is found.
190 !  !
191 !  call eli_m_iter_1d (delta_1d,lambda_1d,maxeigval,nmatsu,z_1d)
192
193 !
194 !  diagonalize M brute forcefully
195 !
196    call eli_diag_m_1d(delta_1d,lambda_1d,maxeigval,mustar,nmatsu,tc,z_1d)
197
198    write (unit_del,'(a)') '#'
199    write (unit_del,'(a)') '# eliashberg_1d : delta_1d = '
200    write (unit_del,'(a)') '#'
201    write (unit_del,'(a,i6,a)') '# delta_1d array containing 2*', nmatsu, '+1 Matsubara frequency points'
202    write (unit_z,'(a,E16.6,a,E16.6)') '# from ', -omega_cutoff, ' to ', omega_cutoff
203    write (unit_z,'(a)') '# delta_1d is the gap function in the Eliashberg equations'
204    write (unit_z,'(a)') '#'
205    do imatsu=-nmatsu,nmatsu
206      write (unit_del,*) imatsu,delta_1d(imatsu)
207    end do
208    write (unit_del,*)
209
210 !  if eigenvalue is < 1 increase T
211 !  else if eigenvalue is > 1 decrease T
212 !  if eigenvalue ~= 1 stop
213 !
214    if (abs(maxeigval-one) < tol8) then
215      write(std_out,*) 'Eliashberg Tc found = ', tc, ' (Ha) = ', tc/kb_HaK, ' (K)'
216      exit
217    else if (maxeigval > 0.001_dp) then
218      tc = tc * maxeigval
219    else
220      write(std_out,*) 'maxeigval is very small'
221      tc = tc * 1000.0_dp
222    end if
223
224
225  end do
226 !end iiter do
227
228  if (abs(maxeigval-one) > tol8) then
229    write(std_out,*) 'eliashberg_1d : Tc not converged. ', maxeigval, ' /= 1'
230    write(std_out,*) 'Eliashberg Tc nonetheless = ', tc, ' (Ha) = ', tc/kb_HaK, ' (K)'
231  end if
232
233  ABI_DEALLOCATE(lambda_1d)
234  ABI_DEALLOCATE(z_1d)
235  ABI_DEALLOCATE(delta_1d)
236  ABI_DEALLOCATE(mm_1d)
237
238  close (UNIT=unit_z)
239  close (UNIT=unit_lam)
240  close (UNIT=unit_del)
241
242  write(std_out,*) ' eliashberg_1d : end '
243
244 end subroutine eliashberg_1d