ABINIT Discussion Forums

Dear abinit users,

I am trying to calculate the Tc of PbTaSe2, which is experimentally around 3.8K. But, when I try to calculate it, my value is nearly ten times smaller!

Here is my groundstate input file and anaddb file. Any help is greatly appreciated!

ndtset 32

pawxcdev 0

#
# DATASET 1 : make ground state wavefunctions and density
#
tolwfr1 1.0d-14
nline1 8 # This is to expedite the convergence of higher-lying bands
rfphon1 0 # for DS1 do _not_ do perturbation
nqpt1 0 # for DS1 do _not_ do perturbation
prtwf1 1 # need GS wavefunctions for further runs
getwfk1 0
kptopt1 3
# enforce calculation of forces at each SCF step
optforces 1

qpt2 0.00000000E+00 0.00000000E+00 0.00000000E+00
qpt3 1.66666667E-01 0.00000000E+00 0.00000000E+00
qpt4 3.33333333E-01 0.00000000E+00 0.00000000E+00
qpt5 5.00000000E-01 0.00000000E+00 0.00000000E+00
qpt6 3.33333333E-01 1.66666667E-01 0.00000000E+00
qpt7 5.00000000E-01 1.66666667E-01 0.00000000E+00
qpt8 -3.33333333E-01 3.33333333E-01 0.00000000E+00
qpt9 0.00000000E+00 0.00000000E+00 5.00000000E-01
qpt10 1.66666667E-01 0.00000000E+00 5.00000000E-01
qpt11 3.33333333E-01 0.00000000E+00 5.00000000E-01
qpt12 5.00000000E-01 0.00000000E+00 5.00000000E-01
qpt13 3.33333333E-01 1.66666667E-01 5.00000000E-01
qpt14 5.00000000E-01 1.66666667E-01 5.00000000E-01
qpt15 -3.33333333E-01 3.33333333E-01 5.00000000E-01

#
# DS5 get DDK perturbation (stored in GKK matrix element files)
#
tolwfr16 1.0d-14
qpt16 0.0 0.0 0.0
rfphon16 0
rfelfd16 2
prtwf16 0

# WF on full BZ
tolwfr17 1.0d-14
prtwf17 1
getwfk17 1
nstep17 1
nline17 1
nqpt17 0
rfphon17 0

# Compute the GKK correctly (no gauge problem)
# This step is very quick as NSCF+nstep 1
qpt18 0.0 0.0 0.0
prtgkk18 1 # print out gkk files directly
iscf18 -2
tolwfr18 1.0d-14
nstep18 1
nline18 1
getwfk18 17
get1den18 2
#continue for other qpoints

# DS12 get the ddk matrix elements on the full grid
prtgkk32 1
rfphon32 0 # no phonons here
rfelfd32 2 # Activate the calculation of the d/dk perturbation
qpt32 0.0 0.0 0.0 # This is a calculation at the Gamma point
iscf32 -3 # The d/dk perturbation must be treated non SC-ly
tolwfr32 1.0d-14 # Must use tolwfr for non-self-consistent calculations
nstep32 1
nline32 1
getwfk32 17
get1den32 5

#
# general data for all phonon calculations:
#
prepgkk 1 # flag to calculate all perturbations for el-phon calculations
use_nonscf_gkk 0 # enforce old default and scf-calculate all perturbations
nqpt 1 # 1 qpoint at a time
tolvrs 1.e-8 # tolerance on 2DTE convergence: potential^(1) is what we need
getwfk 1 # all other DS get wf from DS1
prtwf 0
rfatpol 1 4 # all atoms are perturbed
rfdir 1 1 1 # all directions of perturbation
rfphon 1

#
# Common data for GS and perturbation runs
#

# The kpoint grid is minimalistic to keep the calculation manageable.
ngkpt 12 12 4
kptopt 3

# As the kinetic energy cutoff
ecut 33
pawecutdg 66
pawovlp -1

# Use a centered grid for the kpoints
nshiftk 1
shiftk 0.0 0.0 0.0

#
# Common data
#

acell 6.5755478745E+00 6.5755478745E+00 1.8026772888E+01
rprim
0.5 sqrt(3/4) 0
-0.5 sqrt(3/4) 0
0 0 1

nband 20

# include metallic occupation function with a small smearing
occopt 7
tsmear 0.001
natom 4
typat 3 1 3 2
xred 2/3 2/3 5/6
2/3 2/3 1/2
2/3 2/3 1/6
1 1 0
nstep 800
ntypat 3
znucl 82 73 34


anaddb file:

# turn on calculation of the electron-phonon quantities
elphflag 1
telphint 0

kptrlatt 12 0 0 0 12 0 0 0 4

elphsmear 0.01

# Path in reciprocal space along which the phonon linewidths
# and band structure will be calculated
nqpath 7
qpath
0.0 0.0 0.0
2/3 1/3 0.0
1/2 0.0 0.0
0.0 0.0 0.0
2/3 1/3 1/2
1/2 0.0 1/2
0.0 0.0 1/2

# Coulomb pseudopotential parameter
mustar 0.1

# Minimalistic qpoint grid

ngqpt 6 6 2

# impose acoustic sum rule in a symmetric way
asr 2
dipdip 1

chneut 0

# bravais lattice necessary
brav 1

# qpt grid
nqshft 1
q1shft 0.0 0.0 0.0
#
ifcflag 1
ifcana 1
ifcout 0
# ifc for all atoms?
natifc 0
atifc 1 2 3

# print dielectric matrix with freq dependence
dieflag 0

# print out eigenvectors and symmetrize dyn matrix
eivec 1

#Wavevector list number 1 (Reduced coordinates and normalization factor)
nph1l 1
qph1l
0.00000000E+00 0.00000000E+00 0.00000000E+00 1

Hello,

we have seen similar effects with some compounds. I think it is due to a rough interpolation of the EPC at low q and omega, which makes the a2F function too large at low w. This unbalances the calculation of omega_log, which is the large lambda limiting maximal value of Tc.

Could you post omega_log, the a2F function, lambda values, and test with different Gaussian smearings as well as the tetrahedron integration?

Also, 12 12 4 is a bit coarse. Even if your system has a large unit cell, unfortunately, the k grid convergence can be quite slow