error in routine broyden(6): factorization

Dear all,
I’m a new learner to EPW. When I try to do epw.x calculations, I met an error in the step of solving isotropic Eliashberg equations on imaginary-axis. I have 2 questions:

  1. How to fix the error ’ error in routine broyden(6): factorization’ ?
  2. Why the estimated Allen-Dynes Tc is so large ? What should I do to get the correct value?
    Here are the input and output.

[b:1qq65uqf] input:[/b:1qq65uqf]
[code:1qq65uqf]–
&inputepw
prefix = ‘h3s’,
amass(1) = 1.008,
amass(2) = 32.06
outdir = ‘./’

ep_coupling = .true.
elph = .true.
kmaps = .false.
epbwrite = .true.
epbread = .false.

epwwrite = .true.
epwread = .false.

etf_mem = 1

nbndsub = 9
nbndskip = 0

wannierize = .true.
num_iter = 10000
dis_win_max = 33.0
dis_win_min = 0.0
proj(1)=‘random’
!proj(2) = ‘S:p’
proj(2) = ‘H:s’

iverbosity = 2

parallel_k = .true.
parallel_q = .false.

eps_acustic = 1.0
ephwrite = .true.

fsthick = 0.4
eptemp = 300
degaussw = 0.10
degaussq = 0.5
nqstep = 500

eliashberg = .true.

liso = .true.
limag = .true.
lpade = .true.
lacon = .true.

conv_thr_iaxis = 1.0d-2
conv_thr_racon = 1.0d-2

wscut = 0.5

nstemp = 3
tempsmin = 190.00
tempsmax = 240.00

nsiter = 500

muc = 0.16

dvscf_dir = ‘./save’

nk1 = 6
nk2 = 6
nk3 = 6

nq1 = 2
nq2 = 2
nq3 = 2

mp_mesh_k = .true.
nkf1 = 12
nkf2 = 12
nkf3 = 12

nqf1 = 12
nqf2 = 12
nqf3 = 12
/
4 cartesian
0.000000000 0.000000000 0.000000000
0.000000000 0.000000000 -0.500000000
0.000000000 -0.500000000 -0.500000000
-0.500000000 -0.500000000 -0.500000000[/code:1qq65uqf]

[b:1qq65uqf]output:[/b:1qq65uqf]
[code:1qq65uqf]
......
......
===================================================================
Solve isotropic Eliashberg equations
===================================================================

 Finish reading .freq file 

              Fermi level (eV) =     3.0765347476E+01
 DOS(states/spin/eV/Unit Cell) =     1.4401140222E-02
        Electron smearing (eV) =     1.0000000000E-01
             Fermi window (eV) =     4.0000000000E-01
 Nr irreducible k-points within the Fermi shell =         7 out of        84
       2 bands within the Fermi window


 Finish reading .egnv file 


 Max nr of q-points =        34


 Finish reading .ikmap files


 Start reading .ephmat files


 Finish reading .ephmat files 

 lambda_max =             0.0450411   lambda_k_max =             0.0349339

 Electron-phonon coupling strength =    0.0153917

 Estimated Allen-Dynes Tc = 1304136.7081148 K for muc =    0.16000

 Estimated BCS superconducting gap =     197.7921979 eV

 temp(  1) =    190.00000 K 

 Solve isotropic Eliashberg equations on imaginary-axis

 Total number of frequency points nsiw (      1 ) =      5

 iter =      1   error =     2.4687234455E+03   Znormi(1) =     1.0000096887E+00   Deltai(1) =    -7.9926272214E-02
 iter =      2   error =     3.1032568548E-05   Znormi(1) =     1.0000323260E+00   Deltai(1) =    -7.9923543453E-02
 Convergence was reached in nsiter =      2

 iaxis_imag   :      0.00s CPU      0.00s WALL (       1 calls)

 Pade approximant of isotropic Eliashberg equations from imaginary-axis to real-axis

 pade =      4   error =     9.6971175492E-04   Re[Znorm(1)] =     1.0000324395E+00   Re[Delta(1)] =    -7.9917933118E-02

 raxis_pade   :      0.00s CPU      0.06s WALL (       1 calls)


 Analytic continuation of isotropic Eliashberg equations from imaginary-axis to real-axis

 Total number of frequency points nsw =   1000

 iter =      1   error =     3.5303127519E+00   Re[Znorm(1)] =     1.0095012619E+00   Re[Delta(1)] =     3.0934545775E-02
 iter =      2   error =     4.4514422372E-02   Re[Znorm(1)] =     1.0107920397E+00   Re[Delta(1)] =     3.1000378597E-02
 iter =      3   error =     2.3226943841E-02   Re[Znorm(1)] =     1.0094155556E+00   Re[Delta(1)] =     3.1121693953E-02
 iter =      4   error =     1.4163699664E-03   Re[Znorm(1)] =     1.0094607209E+00   Re[Delta(1)] =     3.1118151292E-02
 Convergence was reached in nsiter =      4

 raxis_acon   :      0.46s CPU      0.50s WALL (       1 calls)

 itemp =   1   total cpu time :     0.6 secs

 temp(  2) =    215.00000 K 

 Solve isotropic Eliashberg equations on imaginary-axis

 Total number of frequency points nsiw (      2 ) =      4

 iter =      1   error =     3.6474898270E+00   Znormi(1) =     1.0104612381E+00   Deltai(1) =     2.9609774910E-02
 iter =      2   error =     1.8145447398E+01   Znormi(1) =     1.0143183603E+00   Deltai(1) =     1.5311049912E-03
 iter =      3   error =     1.1027266065E+01   Znormi(1) =     1.0143319963E+00   Deltai(1) =    -1.5404809425E-04
 iter =      4   error =     2.6237333860E+00   Znormi(1) =     1.0143320299E+00   Deltai(1) =     9.4249689740E-05
 iter =      5   error =     2.6906797559E+01   Znormi(1) =     1.0143320754E+00   Deltai(1) =    -3.6403709234E-06
 iter =      6   error =     5.2143624565E+01   Znormi(1) =     1.0143320754E+00   Deltai(1) =     7.1038939502E-08
 iter =      7   error =     7.2414067743E+00   Znormi(1) =     1.0143320754E+00   Deltai(1) =    -1.1438373280E-08
 iter =      8   error =     1.7121020893E+01   Znormi(1) =     1.0143320754E+00   Deltai(1) =     7.0686178910E-10
 iter =      9   error =     5.5918259276E+01   Znormi(1) =     1.0143320754E+00   Deltai(1) =    -1.2884922422E-11
 iter =     10   error =     1.1674916477E+02   Znormi(1) =     1.0143320754E+00   Deltai(1) =     1.1042401721E-13
 iter =     11   error =     3.2608593559E+02   Znormi(1) =     1.0143320754E+00   Deltai(1) =    -4.1740120020E-16
 iter =     12   error =     7.2644700849E+03   Znormi(1) =     1.0143320754E+00   Deltai(1) =     5.4313404221E-20
 iter =     13   error =     9.4758877708E+00   Znormi(1) =     1.0143320754E+00   Deltai(1) =     5.6544495803E-21
 iter =     14   error =     3.1271261124E+03   Znormi(1) =     1.0143320754E+00   Deltai(1) =     1.5729283948E-24
 iter =     15   error =     1.5790915528E+03   Znormi(1) =     1.0143320754E+00   Deltai(1) =    -1.4355529987E-27
 iter =     16   error =     4.6182498842E+03   Znormi(1) =     1.0143320754E+00   Deltai(1) =    -2.7129426882E-31
 iter =     17   error =     2.2621787534E+04   Znormi(1) =     1.0143320754E+00   Deltai(1) =    -1.1850080281E-35
 iter =     18   error =     3.7699754882E+04   Znormi(1) =     1.0143320754E+00   Deltai(1) =    -3.1010668206E-40
 iter =     19   error =     7.6081446085E+04   Znormi(1) =     1.0143320754E+00   Deltai(1) =    -4.3192551601E-45
 iter =     20   error =     2.5811095566E+02   Znormi(1) =     1.0143320754E+00   Deltai(1) =    -1.3195073361E-47
 iter =     21   error =     1.1064835446E+02   Znormi(1) =     1.0143320754E+00   Deltai(1) =     1.2010977025E-49
 iter =     22   error =     2.8294154637E+03   Znormi(1) =     1.0143320754E+00   Deltai(1) =    -4.2471506223E-53
 iter =     23   error =     5.0995761014E+03   Znormi(1) =     1.0143320754E+00   Deltai(1) =     8.3233123787E-57
 iter =     24   error =     4.0748297992E+02   Znormi(1) =     1.0143320754E+00   Deltai(1) =     2.0294251698E-59
 iter =     25   error =     2.3054296097E+02   Znormi(1) =     1.0143320754E+00   Deltai(1) =     8.7729024036E-62
 iter =     26   error =     1.5360597909E+04   Znormi(1) =     1.0143320754E+00   Deltai(1) =     5.6867203476E-66
 iter =     27   error =     4.0984645680E+03   Znormi(1) =     1.0143320754E+00   Deltai(1) =     1.3266819502E-69
 iter =     28   error =     9.0543008089E+02   Znormi(1) =     1.0143320754E+00   Deltai(1) =    -1.7051109766E-72
 iter =     29   error =     1.0694719641E+02   Znormi(1) =     1.0143320754E+00   Deltai(1) =    -1.4520313479E-74
 iter =     30   error =     3.6912791082E+02   Znormi(1) =     1.0143320754E+00   Deltai(1) =    -3.9212172986E-77
 iter =     31   error =     9.0500315665E+02   Znormi(1) =     1.0143320754E+00   Deltai(1) =    -4.3289705471E-80
 iter =     32   error =     1.7397491957E+03   Znormi(1) =     1.0143320754E+00   Deltai(1) =    -2.4887526622E-83
 iter =     33   error =     3.3143472843E+03   Znormi(1) =     1.0143320754E+00   Deltai(1) =    -7.5200225422E-87
 iter =     34   error =     1.4183729964E+03   Znormi(1) =     1.0143320754E+00   Deltai(1) =    -5.2789037394E-90
 iter =     35   error =     2.2066633003E+03   Znormi(1) =     1.0143320754E+00   Deltai(1) =     2.3832886622E-93
 iter =     36   error =     6.8285808331E+02   Znormi(1) =     1.0143320754E+00   Deltai(1) =    -3.5043862381E-96
 iter =     37   error =     5.5213745986E+02   Znormi(1) =     1.0143320754E+00   Deltai(1) =     6.3182051807E-99
 iter =     38   error =     1.7660910107E+02   Znormi(1) =     1.0143320754E+00   Deltai(1) =     3.4704250772-101
 iter =     39   error =     2.3623341737E+02   Znormi(1) =     1.0143320754E+00   Deltai(1) =    -8.5142226839-104
 iter =     40   error =     7.5477225197E+00   Znormi(1) =     1.0143320754E+00   Deltai(1) =    -2.0277536460-104
 iter =     41   error =     3.1199926982E+03   Znormi(1) =     1.0143320754E+00   Deltai(1) =     6.3583812063-108
 iter =     42   error =     1.4744014679E+03   Znormi(1) =     1.0143320754E+00   Deltai(1) =    -4.3198221105-111
 iter =     43   error =     3.1923011944E+03   Znormi(1) =     1.0143320754E+00   Deltai(1) =     1.3521215438-114
 iter =     44   error =     2.0201245912E+03   Znormi(1) =     1.0143320754E+00   Deltai(1) =    -6.7050212014-118

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Error in routine broyden (6):
factorization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 stopping ...

[/code:1qq65uqf]

Thanks and regards,
Xiaohan

Dear Xiaohan,

To answer your questions:

  1. You can lower the broyden_beta parameter (default value 0.7). However, in your case the code is crashing most likely because you are attempting to solve the Eliashberg equation at temperatures greater than the Tc (estimated Tc based on Allen-Dynes is 197 K). You should start solving the equations from lower temperatures and look at the temperature dependence of the superconducting gap [look for instance at Fig. 4a in PHYSICAL REVIEW B 87, 024505 (2013)].

  2. I would guess that your calculations may not be converged. You need to converge the e-ph coupling strength with respect to the k and q-meshes. The Allen-Dynes expression on Tc also depends on which value you use for the Coulomb parameter mu*.

Best,
Roxana

Dear Roxana,

Thank you very much for your reply.
The first problem is truly due to that the temperature is greater than the Tc.
And the second problem is more likely related to my wrong setting of nbndsub and nbndskip. I recalculated epw with the same k , q points and muc as the article ([url:3dc9xphq]https://arxiv.org/abs/1606.09477[/url:3dc9xphq] ). I set nbndsub=12 , nbndskip=4 . And the Tc is 118K, which is still smaller than that in the article(166K).
So I want to ask what other converge tests should I take?About which parameters?
Thanks again.

Best,
Xiaohan

Dear Xiaohan,

For example you can check if the wannierization is correct (e.g, band structure after wannier interpolation) or if the decay of the Hamiltonian, dynamical matrix e-ph matrix elements are good enough [see Computer Physics Communications 181 (2010) 21402148].

Best,
Roxana