TABLE OF CONTENTS


ABINIT/m_eliashberg_1d [ Modules ]

[ Top ] [ 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