Graph problem of Anisotropic superconducting gaps of mgb2

[Mu Haimen]

Dear all
I have a question that quite confused me when I calculated MgB2's superconducting gaps following the case in EPW examples.
The first figure is from my calculation, and the second is from article DOI:10.1103/PhysRevB.87.024505.

The question is why my gaps are thin (planus) and long while the ones in the literature are more concentrated (or sharper) ? How can I set the proper input parameters to obtain the same result as the figure from the PRB article ?
k- q-mesh in the 2 calculations are the same, both using 6^3 coarse k-mesh and 6^3 coarse q-mesh, 60^3 fine k-mesh and 30^3 fine q-mesh, and k-mesh in scf step is 24^3.





This is my epw.in file:

Code: Select all

--
&inputepw
  prefix      = 'MgB2',
  amass(1)    = 24.305,
  amass(2)    = 10.811
  outdir      = './'

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

  epwwrite = .true.
  epwread  = .false.

  etf_mem     =  1 

  nbndsub     =  5,
  nbndskip    =  0

  wannierize  = .true.
  num_iter    = 500
  dis_froz_max= 8.8
  proj(1)     = 'B:pz'
  proj(2)     = 'f=0.5,1.0,0.5:s'
  proj(3)     = 'f=0.0,0.5,0.5:s'
  proj(4)     = 'f=0.5,0.5,0.5:s'

  iverbosity  = 2

  eps_acustic = 2.0    ! Lowest boundary for the 
  ephwrite    = .true. ! Writes .ephmat files used when wliasberg = .true.

  fsthick     = 0.4  ! eV
  eptemp      = 300  ! K
  degaussw    = 0.10 ! eV
  nsmear      = 1
  delta_smear = 0.04 ! eV

  degaussq     = 0.5 ! meV
  nqstep       = 500

  eliashberg  = .true.

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

  conv_thr_iaxis = 1.0d-5

  wscut = 1.0   ! eV   Upper limit over frequency integration/summation in the Elisashberg eq

  temps(1) = 7
  temps(2) = 10
  temps(3) = 15
  temps(4) = 20
  temps(5) = 25
  temps(6) = 30
  temps(7) = 35
  temps(8) = 40
  temps(9) = 45
  temps(10) = 47.5
  temps(11) = 50

  nsiter   = 500

  muc     = 0.16

  dvscf_dir   = '../phonons/save'
  
  nk1         = 6
  nk2         = 6
  nk3         = 6

  nq1         = 6
  nq2         = 6
  nq3         = 6

  mp_mesh_k = .true.
  nkf1 = 60
  nkf2 = 60
  nkf3 = 60

  nqf1 = 30
  nqf2 = 30
  nqf3 = 30
 /
  28 cartesian
  0.000000000   0.000000000   0.000000000  0.004629630
  0.000000000   0.000000000   0.145933920  0.009259259
  0.000000000   0.000000000   0.291867841  0.009259259
  0.000000000   0.000000000  -0.437801761  0.004629630
  0.000000000   0.192450090   0.000000000  0.027777778
  0.000000000   0.192450090   0.145933920  0.055555556
  0.000000000   0.192450090   0.291867841  0.055555556
  0.000000000   0.192450090  -0.437801761  0.027777778
  0.000000000   0.384900179   0.000000000  0.027777778
  0.000000000   0.384900179   0.145933920  0.055555556
  0.000000000   0.384900179   0.291867841  0.055555556
  0.000000000   0.384900179  -0.437801761  0.027777778
  0.000000000  -0.577350269   0.000000000  0.013888889
  0.000000000  -0.577350269   0.145933920  0.027777778
  0.000000000  -0.577350269   0.291867841  0.027777778
  0.000000000  -0.577350269  -0.437801761  0.013888889
  0.166666667   0.288675135   0.000000000  0.027777778
  0.166666667   0.288675135   0.145933920  0.055555556
  0.166666667   0.288675135   0.291867841  0.055555556
  0.166666667   0.288675135  -0.437801761  0.027777778
  0.166666667   0.481125224   0.000000000  0.055555556
  0.166666667   0.481125224   0.145933920  0.111111111
  0.166666667   0.481125224   0.291867841  0.111111111
  0.166666667   0.481125224  -0.437801761  0.055555556
  0.333333333   0.577350269   0.000000000  0.009259259
  0.333333333   0.577350269   0.145933920  0.018518519
  0.333333333   0.577350269   0.291867841  0.018518519
  0.333333333   0.577350269  -0.437801761  0.009259259

Best wishes
Bruce Mu

[sponce]

Dear Bruce Mu,

It seems indeed you are using the same parameters. However the code was optimized, bugs were corrected etc.. The two results are indeed different but look somewhat similar.

If you want to get to the bottom of this, you should try to converge the results (increase the fine grids, decrease degaussw etc).

From the original paper from Roxana, the 60/30 grids were not yet fully converged so I would start there.

Maybe she has additional suggestions as well.

Best wishes,
Samuel

[Mu Haimen]

Dear Samuel,

Thank you for your reply , I will try it following your tips.

Best regards,
Bruce Mu

[roxana]

Hi,

I can add that the calculations in PRB were done with Quantum Espresso v4.0.3 and EPW v3.0. One major change between these older versions of the codes and the newer versions (QE starting from v5.* and EPW starting from 4.*) was some change in symmetry in the phonon calculations.

You mentioned that you used the same k- and q-grids as in the PRB. Are you also using the same cutoff of 60 Ry as in the paper?

Best,
Roxana

[Mu Haimen]

Hi, Roxana,

Indeed I used 60Ry cutoff in the scf calculation. My QE version is v6.3 and EPW is v5.0.0. I don't known if there is other option can affect the shape of gaps.

Below are my scf-, nscf-step and phonons-step input files:
scf.in:

Code: Select all

 &control
    calculation='scf',
    prefix='MgB2',
    pseudo_dir = '../pp/',
    outdir='./',
    tprnfor = .true.,
    tstress = .true.,
    etot_conv_thr = 1.0d-5
    forc_conv_thr = 1.0d-4
 /
 &system
    ibrav = 4,
    celldm(1) = 5.8260252227888,
    celldm(3) = 1.1420694129095,
    nat=  3,
    ntyp = 2,
    ecutwfc = 60
    smearing = 'mp'
    occupations = 'smearing'
    degauss = 0.02
 /
 &electrons
    diagonalization = 'david'
    mixing_mode = 'plain'
    mixing_beta = 0.7
    conv_thr =  1.0d-9
 /
ATOMIC_SPECIES
 Mg  24.305  Mg.pz-n-vbc.UPF
 B   10.811  B.pz-vbc.UPF
ATOMIC_POSITIONS crystal
Mg       0.000000000   0.000000000   0.000000000
B        0.333333333   0.666666667   0.500000000
B        0.666666667   0.333333333   0.500000000
K_POINTS AUTOMATIC
24 24 24  0 0 0

nscf.in:

Code: Select all

 &control
    calculation='nscf',
    prefix='MgB2',
    pseudo_dir = '../pp/',
    outdir='./',
    tprnfor = .true.,
    tstress = .true.,
    etot_conv_thr = 1.0d-5
    forc_conv_thr = 1.0d-4
 /
 &system
    ibrav = 4,
    celldm(1) = 5.8260252227888,
    celldm(3) = 1.1420694129095,
    nat=  3,
    ntyp = 2,
    ecutwfc = 60
    smearing = 'mp'
    occupations = 'smearing'
    degauss = 0.02
 /
 &electrons
    diagonalization = 'david'
    mixing_mode = 'plain'
    mixing_beta = 0.7
    conv_thr =  1.0d-9
 /
ATOMIC_SPECIES
 Mg  24.305  Mg.pz-n-vbc.UPF
 B   10.811  B.pz-vbc.UPF
ATOMIC_POSITIONS crystal
Mg       0.000000000   0.000000000   0.000000000
B        0.333333333   0.666666667   0.500000000
B        0.666666667   0.333333333   0.500000000
K_POINTS crystal
216
  0.00000000  0.00000000  0.00000000  4.629630e-03
  0.00000000  0.00000000  0.16666667  4.629630e-03
  ... ... ... ...

ph.in:

Code: Select all

--
&inputph
  prefix   = 'MgB2',
  recover  = .true.
  fildyn   = 'MgB2.dyn',
  amass(1) = 24.305,
  amass(2) = 10.811,
  outdir   = './'
  ldisp    = .true.,
  trans    = .true.,
  fildvscf = 'dvscf',
  nq1=6,
  nq2=6,
  nq3=6,
  tr2_ph   =  1.0d-16
 /

Thank you!

Best regards,
Bruce Mu

[roxana]

Hi,

I found the reason behind the difference in the gap plots. Apparently, at some point I changed the post processing step when writing on file the superconducting gaps on FS and didn't take into account the distance of the states with respect to the Fermi level. I reverted to the original version of the code (EPW 3) and the results for MgB2 are now the same as in my PRB.

The changes have been included in the developer version on gitlab and will also be part of a new version of QE that will be released at the end of this week (QE 6.4.1).

Best,
Roxana