Graph problem of Anisotropic superconducting gaps of mgb2

General
1

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.

    [img:x1yk3mlp]https://gss0.baidu.com/-fo3dSag_xI4khGko9WTAnF6hhy/zhidao/pic/item/d1a20cf431adcbef30ee5411a2af2edda2cc9fca.jpg[/img:x1yk3mlp]

This is my epw.in file:
[code:x1yk3mlp]–
&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
[/code:x1yk3mlp]

Best wishes
Bruce Mu

2

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

3

Dear Samuel,

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

Best regards,
Bruce Mu

4

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

5

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:1sibm0yt] &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
[/code:1sibm0yt]
nscf.in:
[code:1sibm0yt] &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
... ... ... ...
[/code:1sibm0yt]
ph.in:
[code:1sibm0yt]–
&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
/
[/code:1sibm0yt]

Thank you!

Best regards,
Bruce Mu

6

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