Gap energy abruptly goes negative at high temperatures

[CHILLL]

Dear experts and users,

I calculated superconducting properties of a material in EPW 5.4 and 5.1, and there is the same problem I see.

The isotropic results

Code: Select all

     Electron-phonon coupling strength =    0.8613774
 
     Estimated Allen-Dynes Tc =    17.174912 K for muc =    0.13000
  
     Estimated w_log in Allen-Dynes Tc =    33.076210 meV
  
     Estimated BCS superconducting gap =     2.604834 meV

look quite alright, but when the code solves the anisotropic Eliashberg equations I get negative gap energy for higher temperatures.
For example, up to 7 K everything looks alright. From the outpput file I get:

Code: Select all

Min. value of superconducting gap =     0.000000 meV

But at 10 K the situation changes:

Code: Select all

Min. value of superconducting gap =    -4.879436 meV

As you can see, I get approx. -5 meV gap energy. The gap distribution function also looks very different from the one for lower temperatures.

Could you please tell me if there is something I can check ? My wannier functions are good, and the phonons are stable.
Here is my EPW input without the wannier block:

Code: Select all

--
&inputepw
  prefix      = 'structure',
  outdir      = './'
 
 !this is a restart calculation
  ep_coupling = .false.
  elph        = .false.
  !! kmaps       = .true.
  epbwrite    = .false.
  epbread     = .false.
  epwwrite    = .false.
  epwread     = .true.


  etf_mem     =  1

  nbndsub     =  24 !Number of Wannier functions to utilize
 !! nbndskip    =  0
 
  iverbosity  = 2

!  ngaussw = -99
  
  fsthick     = 0.4  ! eV
 !! eptemp      = 300  ! K
  degaussw    = 0.1  ! eV

  ephwrite    = .false.

  eliashberg  = .true.

  laniso = .true.
  limag = .true.
  !lpade = .true.

  nsiter   = 1500
  conv_thr_iaxis = 1.0d-5

  wscut = 0.5  ! eV

  temps = 7 10 ! Critical temp ~ 6K

   muc     = 0.13

  dvscf_dir   = '../phonon/save'

 !nscf (Should be commensurate, so watch out!)
  nk1         = 7
  nk2         = 8
  nk3         = 1
 !ph
  nq1         = 7
  nq2         = 8
  nq3         = 1

  mp_mesh_k = .true.
  nkf1 = 140
  nkf2 = 160
  nkf3 = 1

  nqf1 = 70
  nqf2 = 80
  nqf3 = 1 
/
  20 cartesian ! like in ".dyn0" file!
   0.000000000000000E+00   0.000000000000000E+00   0.000000000000000E+00
   0.000000000000000E+00   0.150306232485474E+00   0.000000000000000E+00
   0.000000000000000E+00   0.300612464970948E+00   0.000000000000000E+00
   0.000000000000000E+00   0.450918697456422E+00   0.000000000000000E+00
   0.000000000000000E+00  -0.601224929941896E+00   0.000000000000000E+00
   0.142857142857143E+00   0.000000000000000E+00   0.000000000000000E+00
   0.142857142857143E+00   0.150306232485474E+00   0.000000000000000E+00
   0.142857142857143E+00   0.300612464970948E+00   0.000000000000000E+00
   0.142857142857143E+00   0.450918697456422E+00   0.000000000000000E+00
   0.142857142857143E+00  -0.601224929941896E+00   0.000000000000000E+00
   0.285714285714286E+00   0.000000000000000E+00   0.000000000000000E+00
   0.285714285714286E+00   0.150306232485474E+00   0.000000000000000E+00
   0.285714285714286E+00   0.300612464970948E+00   0.000000000000000E+00
   0.285714285714286E+00   0.450918697456422E+00   0.000000000000000E+00
   0.285714285714286E+00  -0.601224929941896E+00   0.000000000000000E+00
   0.428571428571429E+00   0.000000000000000E+00   0.000000000000000E+00
   0.428571428571429E+00   0.150306232485474E+00   0.000000000000000E+00
   0.428571428571429E+00   0.300612464970948E+00   0.000000000000000E+00
   0.428571428571429E+00   0.450918697456422E+00   0.000000000000000E+00
   0.428571428571429E+00  -0.601224929941896E+00   0.000000000000000E+00

Thank you!
Best,
Mikhail

[CHILLL]

A little update to the question.
In order to understand the problem better, I computed the EPC projected on the phonon band structure. and there is something I am concerned about.
In the output file(see below) I see that lamda_tr and gamma_tr go negative. I think I have never seen it before. Could you please explain what exactly is the difference between lamda and lamda_tr ? And is it ok to have such negative values ?

Code: Select all

   ismear =     1 iq =       1 coord.:   0.00000  0.00000  0.00000 wt:   1.00000 Temp:    2.500K
     -------------------------------------------------------------------
     lambda___(   1 )=       0.000000   gamma___=       0.000000 meV   omega=      0.5760 meV
     lambda_tr(   1 )=       0.000000   gamma_tr=       0.000000 meV   omega=      0.5760 meV
     lambda___(   2 )=       0.000000   gamma___=       0.000000 meV   omega=      0.2633 meV
     lambda_tr(   2 )=       0.000000   gamma_tr=       0.000000 meV   omega=      0.2633 meV
     lambda___(   3 )=       0.000000   gamma___=       0.000000 meV   omega=      0.0000 meV
     lambda_tr(   3 )=       0.000000   gamma_tr=       0.000000 meV   omega=      0.0000 meV
     lambda___(   4 )=       0.000007   gamma___=       0.000042 meV   omega=     51.8795 meV
     lambda_tr(   4 )=      -0.000000   gamma_tr=      -0.000000 meV   omega=     51.8795 meV
     lambda___(   5 )=       0.000071   gamma___=       0.000481 meV   omega=     56.7642 meV
     lambda_tr(   5 )=       0.000001   gamma_tr=       0.000004 meV   omega=     56.7642 meV
     lambda___(   6 )=       0.000019   gamma___=       0.000143 meV   omega=     59.9750 meV
     lambda_tr(   6 )=      -0.000000   gamma_tr=      -0.000003 meV   omega=     59.9750 meV
     lambda___(   7 )=       0.001571   gamma___=       0.015966 meV   omega=     69.6555 meV
     lambda_tr(   7 )=      -0.000085   gamma_tr=      -0.000863 meV   omega=     69.6555 meV
     lambda___(   8 )=       0.000001   gamma___=       0.000023 meV   omega=     90.0875 meV
     lambda_tr(   8 )=      -0.000000   gamma_tr=      -0.000000 meV   omega=     90.0875 meV
     lambda___(   9 )=       0.000015   gamma___=       0.000262 meV   omega=     91.6168 meV
     lambda_tr(   9 )=      -0.000000   gamma_tr=      -0.000002 meV   omega=     91.6168 meV
     lambda___(  10 )=       0.000140   gamma___=       0.003877 meV   omega=    115.0273 meV
     lambda_tr(  10 )=      -0.000041   gamma_tr=      -0.001124 meV   omega=    115.0273 meV
     lambda___(  11 )=       3.664260   gamma___=     125.245782 meV   omega=    127.7233 meV
     lambda_tr(  11 )=       0.016005   gamma_tr=       0.547054 meV   omega=    127.7233 meV
     lambda___(  12 )=       0.000000   gamma___=       0.000006 meV   omega=    138.7679 meV
     lambda_tr(  12 )=       0.000000   gamma_tr=       0.000008 meV   omega=    138.7679 meV
     lambda___(  13 )=       0.000260   gamma___=       0.010737 meV   omega=    140.4893 meV
     lambda_tr(  13 )=      -0.000027   gamma_tr=      -0.001098 meV   omega=    140.4893 meV
     lambda___(  14 )=       0.019133   gamma___=       0.942512 meV   omega=    153.3304 meV
     lambda_tr(  14 )=      -0.030710   gamma_tr=      -1.512761 meV   omega=    153.3304 meV
     lambda___(  15 )=       0.147621   gamma___=       7.965514 meV   omega=    160.4776 meV
     lambda_tr(  15 )=       0.000251   gamma_tr=       0.013526 meV   omega=    160.4776 meV
     lambda___(  16 )=       0.001241   gamma___=       0.090714 meV   omega=    186.8010 meV
     lambda_tr(  16 )=      -0.000489   gamma_tr=      -0.035739 meV   omega=    186.8010 meV
     lambda___(  17 )=       0.000796   gamma___=       0.058191 meV   omega=    186.8331 meV
     lambda_tr(  17 )=      -0.000313   gamma_tr=      -0.022884 meV   omega=    186.8331 meV
     lambda___(  18 )=       0.013999   gamma___=       1.215107 meV   omega=    203.5321 meV
     lambda_tr(  18 )=       0.001747   gamma_tr=       0.151636 meV   omega=    203.5321 meV
     lambda___( tot )=       3.849134
     lambda_tr( tot )=      -0.013661
     -------------------------------------------------------------------

Thank you!
Mikhail

[hpaudya1]

Hi Mikhail,

Could you please provide the following two plots/data?
1) The gap distribution function at every temps you have.
2) Scattered plot of \Delta_nk vs E- EF at every temps you have. This is the plot for example in MgB2, 5th, and 6th column of "MgB2.imag_aniso_gap_FS_15.00" at 15 K temp.

Best,
Hari

[CHILLL]

Hello Hari,

Here are the pictures you requested:
https://ibb.co/XV4C5CC

https://ibb.co/Qp0XTfg


Best,
Mikhail

[hpaudya1]

Hi Mikhail,

Numerically \Delta can have any values, we do not have any restriction in the code, however, a negative value of \Delta is not physical. I think the \Delta becomes negative if the following term is negative

Code: Select all

lambda(k, k', n-n') - muc*

see eq. (22) of https://journals.aps.org/prb/abstract/1 ... .87.024505
I believe, for muc* = 0, there will be no issue. So I suggest you run some tests with muc* = 0 (or small values than what you are using now), and see what happens.

I hope it helps.

Best,
Hari

[CHILLL]

Hello Hari,

I did those tests already, and you are absolutely right: mustar 0 solves the problem.
And when I start increasing it, the problem appears again at a value around 0.1(likely less than that).
Unfortunately, results with mustar=0 are hardly suitable for a research, and there must be a reason why lambda(k, k', n-n') goes so low.
Is there anything else you could recommend for me to try ?

Best,
Mikhail

[hpaudya1]

Hi Makhail,

One thing I am thinking right now is about 'wscut'. You have used 0.5 eV, can you say how did you pick this value? Is your maximum phonon frequency 205 meV? If yes, wscut = 0.5 eV is very low.

You can do the following,
1) try increasing wscut (you can get an idea from the MgB2 calculation we have done, where we have used 1 eV which is around 10 times the max. phonon freq.)
Please not that, as you increase temp., the number of frequency grid points (nswi) decreases for a fixed value of wscut.

2) Your 6K-7K calculations look ok with wscut = 0.5eV. keep nsiw = XX in your input and run at 9 K (XX from your 1st temp, so that it will ignore wscut)

Hope it helps.
Best,
Hari

[CHILLL]

Hello Hari,

Your advice seems to be working, thanks a lot !
The only thing is that the higher the temperature I set in the calculation the more frequencies I need to make the gap positive.
Is this normal behavior ?
For all my previous materials, I used to just set a certain wscut (a very small one compared to what I need now), and the I just use it for all the temperatures. Is it not the optimal way then ?

Best,
Mikhail

[hpaudya1]

Hi Mikhail,

That's great.

The only thing is that the higher the temperature I set in the calculation the more frequencies I need to make the gap positive. Is this normal behavior ?

I think so. At the temperature close to Tc, the gap value will be low and may require a larger wscut.

For all my previous materials, I used to just set a certain wscut (a very small one compared to what I need now), and the I just use it for all the temperatures. Is it not the optimal way then ?

You need to do a convergence test for the wscut for the particular system you are interested in.

Happy EPWing,
Hari

[CHILLL]

Hello Hari,

At the moment, for 16 K temperature I am trying

Code: Select all

     Total number of frequency points nsiw(     1) =    550
     Cutoff frequency wscut =     4.7690

with phonons not higher in energy than 200 meV.
I am trying such high value because

Code: Select all

     Total number of frequency points nsiw(     1) =    450
     Cutoff frequency wscut =     3.9027

still produces negative gap.

The problem is, the calculations become very slow and I am running out of available time for a job on my cluster. So could you please answer two questions:
1) Is there a way to run anisotropic Eliashberg calculations with more processors than used for ephmat calculations ? When I try to increase it I get an error.
2) Can I tune some other parameter that could help with faster wscut convergence ? For example, could denser interpolation grids decrease the necessary wscut value ?

Thank you!
Best,
Mikhail