Erroneous phonon band structure of SrTiO3 using PAWs
Posted: Tue Aug 19, 2014 3:10 pm
Hi,
I am attempting to calculate the phonon band structure of cubic SrTiO3 with DFPT. Using the PAW potentials generated by AtomPAW2Abinit and the Perdew-Wang functional I obtain the erroneous result show in the attached figure phonon_band_structure_PAW.png. Although I have carried out the correction for the TO-LO splitting, the dispersion is still not smooth at Gamma. Moreover, the phonon band structure calculated using the PAW potentials do not show the negative frequencies in the Gamma and M points due to the cubic phase not being the ground state, for instance reported in Ph. Ghosez et al., AIP Conf. Proc. 535 (1), 102 (2000). Since the latter work was carried out using norm-conserving PSP I've redone the calculation using OPIUM generated PSP and the PBE functional. The result is shown in phonon_band_struct_PSP.png. In this case the LO-TO splitting is correctly accounted for, and the band structure is generally in good agreement with that reported in Ghosez et al. Yet, there is an extra instability around X in my calculation, which the reference shows for BaTiO3 and Ba0.5Sr0.5TiO3, but not for SrTiO3. I'm guessing this is due to the use of PBE functional in my calculation which introduces a relative error of +1% on the lattice parameter compared with the use of the LDA in the reference, which I found produces an error of -1% on the lattice parameters in case of SrTiO3.
Does the implementation of DFPT with PAWs still have some issues in the current version Abinit 7.8.2, in particular with regard to the TO-LO splitting? I hope someone can clarify these findings.
My script for the generation of the derivative database using PAWs (based on the input file from the second tutorial on response calculations) and the input file for the analysis of the ddb, to obtain the phonon band structure are displayed below.
Best wishes,
Jonas
******************************************************************
******************************************************************
# Crystalline SrTiO3 : computation of the phonon spectrum
ndtset 22
#Set 1 : ground state self-consistency
getwfk1 0 # Cancel default
kptopt1 1 # Automatic generation of k points, taking
# into account the symmetry
nqpt1 0 # Cancel default
tolwfr1 1.0d-22 # SCF stopping criterion (modify default)
rfphon1 0 # Cancel default
iscf1 17
nstep1 1600
prtwf1 1
#Q vectors for all datasets
#Complete set of symmetry-inequivalent qpt chosen to be commensurate
# with kpt mesh so that only one set of GS wave functions is needed.
#Generated automatically by running GS calculation with kptopt=1,
# nshift=0, shiftk=0 0 0 (to include gamma) and taking output kpt set
# file as qpt set. Set nstep=1 so only one iteration runs.
nqpt 1 # One qpt for each dataset (only 0 or 1 allowed)
# This is the default for all datasets and must
# be explicitly turned off for dataset 1.
qpt2 0.00000000E+00 0.00000000E+00 0.00000000E+00
qpt3 0.00000000E+00 0.00000000E+00 0.00000000E+00
qpt4 1.66666667E-01 0.00000000E+00 0.00000000E+00
qpt5 3.33333333E-01 0.00000000E+00 0.00000000E+00
qpt6 5.00000000E-01 0.00000000E+00 0.00000000E+00
qpt7 1.66666667E-01 1.66666667E-01 0.00000000E+00
qpt8 3.33333333E-01 1.66666667E-01 0.00000000E+00
qpt9 5.00000000E-01 1.66666667E-01 0.00000000E+00
qpt10 3.33333333E-01 3.33333333E-01 0.00000000E+00
qpt11 5.00000000E-01 3.33333333E-01 0.00000000E+00
qpt12 5.00000000E-01 5.00000000E-01 0.00000000E+00
qpt13 1.66666667E-01 1.66666667E-01 1.66666667E-01
qpt14 3.33333333E-01 1.66666667E-01 1.66666667E-01
qpt15 5.00000000E-01 1.66666667E-01 1.66666667E-01
qpt16 3.33333333E-01 3.33333333E-01 1.66666667E-01
qpt17 5.00000000E-01 3.33333333E-01 1.66666667E-01
qpt18 5.00000000E-01 5.00000000E-01 1.66666667E-01
qpt19 3.33333333E-01 3.33333333E-01 3.33333333E-01
qpt20 5.00000000E-01 3.33333333E-01 3.33333333E-01
qpt21 5.00000000E-01 5.00000000E-01 3.33333333E-01
qpt22 5.00000000E-01 5.00000000E-01 5.00000000E-01
#Set 2 : Response function calculation of d/dk wave function
iscf2 -3 # Need this non-self-consistent option for d/dk
kptopt2 2 # Modify default to use time-reversal symmetry
rfphon2 0 # Cancel default
rfelfd2 2 # Calculate d/dk wave function only
tolwfr2 1.0d-22 # Use wave function residual criterion instead
nstep2 1600
prtwf2 1
#Set 3 : Response function calculation of Q=0 phonons and electric field pert.
getddk3 2 # d/dk wave functions from last dataset
kptopt3 2 # Modify default to use time-reversal symmetry
rfelfd3 3 # Electric-field perturbation response only
#Sets 4-22 : Finite-wave-vector phonon calculations (defaults for all datasets)
getwfk 1 # Use GS wave functions from dataset1
kptopt 3 # Need full k-point set for finite-Q response
rfphon 1 # Do phonon response
rfatpol 1 5 # Treat displacements of all atoms
rfdir 1 1 1 # Do all directions (symmetry will be used)
tolvrs 1.0d-8 # This default is active for sets 3-22
prtwf 0
#######################################################################
#Common input variables
# Definition of the unit cell
#acell 3*7.27888848134371
acell 3*7.27898250385070
rprim # a cubic unit cell
1 0 0
0 1 0
0 0 1
# Definition of the atom types
ntypat 3
znucl 38 22 8 # strontium (38), titanium (22), oxygen (8)
# Definition of the atoms
natom 5
typat 1 2 3 3 3
xred
0.0 0.0 0.0 # Sr sits at cubic site
0.5 0.5 0.5 # Ti sits at body-center site
0.0 0.5 0.5 # O sits at face-centered sites
0.5 0.0 0.5
0.5 0.5 0.0
#Gives the number of band, explicitely (do not take the default)
nband 20
#Exchange-correlation functional
ixc 7 # LDA-Perdew-Wang
#Definition of the planewave basis set
ecut 60 # Maximal kinetic energy cut-off, in Hartree
pawecutdg 80
#Definition of the k-point grid
ngkpt 6 6 6
nshiftk 1 # Use one copy of grid only (default)
shiftk 0.5 0.5 0.5
paral_kgb 0
#Definition of the SCF procedure
iscf 7 # Self-consistent calculation, using algorithm 5
nstep 400 # Maximal number of SCF cycles
diemac 6.0 # Although this is not mandatory, it is worth to
# precondition the SCF cycle. The model dielectric
# function used as the standard preconditioner
# is described in the "dielng" input variable section.
# The dielectric constant of AlAs is smaller that the one of Si (=12).
# add to conserve old < 6.7.2 behavior for calculating forces at each SCF step
optforces 1
******************************************************************
******************************************************************
!Input file for the anaddb code. Analysis of the SrTiO3 DDB
!Flags
ifcflag 1 ! Interatomic force constant flag
!Wavevector grid number 1 (coarse grid, from DDB)
brav 1 ! Bravais Lattice : 1-S.C., 2-F.C., 3-B.C., 4-Hex.)
ngqpt 6 6 6 ! Monkhorst-Pack indices
nqshft 1 ! number of q-points in repeated basic q-cell
q1shft 3*0.0
!Effective charges
asr 1 ! Acoustic Sum Rule. 1 => imposed asymetrically
chneut 1 ! Charge neutrality requirement for effective charges.
!Interatomic force constant info
dipdip 1 ! Dipole-dipole interaction treatment
!Phonon band structure output for band2eps - See note near end for
! dealing with gamma LO-TO splitting issue.
eivec 4
!Wavevector list number 1 (Reduced coordinates and normalization factor)
nph1l 101 ! number of phonons in list 1
qph1l 0 0 0 1.0
0 0.025 0 1.0
0 0.05 0 1.0
0 0.075 0 1.0
0 0.1 0 1.0
0 0.125 0 1.0
0 0.15 0 1.0
0 0.175 0 1.0
0 0.2 0 1.0
0 0.225 0 1.0
0 0.25 0 1.0
0 0.275 0 1.0
0 0.3 0 1.0
0 0.325 0 1.0
0 0.35 0 1.0
0 0.375 0 1.0
0 0.4 0 1.0
0 0.425 0 1.0
0 0.45 0 1.0
0 0.475 0 1.0
0 0.5 0 1.0
0.025 0.5 0 1.0
0.05 0.5 0 1.0
0.075 0.5 0 1.0
0.1 0.5 0 1.0
0.125 0.5 0 1.0
0.15 0.5 0 1.0
0.175 0.5 0 1.0
0.2 0.5 0 1.0
0.225 0.5 0 1.0
0.25 0.5 0 1.0
0.275 0.5 0 1.0
0.3 0.5 0 1.0
0.325 0.5 0 1.0
0.35 0.5 0 1.0
0.375 0.5 0 1.0
0.4 0.5 0 1.0
0.425 0.5 0 1.0
0.45 0.5 0 1.0
0.475 0.5 0 1.0
0.5 0.5 0 1.0
.475 .475 0 1.0
.45 .45 0 1.0
.425 .425 0 1.0
.4 .4 0 1.0
.375 .375 0 1.0
.35 .35 0 1.0
.325 .325 0 1.0
.3 .3 0 1.0
.275 .275 0 1.0
.25 .25 0 1.0
.225 .225 0 1.0
.2 .2 0 1.0
.175 .175 0 1.0
.15 .15 0 1.0
.125 .125 0 1.0
.1 .1 0 1.0
.075 .075 0 1.0
.05 .05 0 1.0
.025 .025 0 1.0
0 0 0 1.0
0.025 0.025 0.025 1.0
0.05 0.05 0.05 1.0
0.075 0.075 0.075 1.0
0.1 0.1 0.1 1.0
0.125 0.125 0.125 1.0
0.15 0.15 0.15 1.0
0.175 0.175 0.175 1.0
0.2 0.2 0.2 1.0
0.225 0.225 0.225 1.0
0.25 0.25 0.25 1.0
0.275 0.275 0.275 1.0
0.3 0.3 0.3 1.0
0.325 0.325 0.325 1.0
0.35 0.35 0.35 1.0
0.375 0.375 0.375 1.0
0.4 0.4 0.4 1.0
0.425 0.425 0.425 1.0
0.45 0.45 0.45 1.0
0.475 0.475 0.475 1.0
0.5 0.5 0.5 1.0
.475 0.5 .475 1.0
.45 0.5 .45 1.0
.425 0.5 .425 1.0
.4 0.5 .4 1.0
.375 0.5 .375 1.0
.35 0.5 .35 1.0
.325 0.5 .325 1.0
.3 0.5 .3 1.0
.275 0.5 .275 1.0
.25 0.5 .25 1.0
.225 0.5 .225 1.0
.2 0.5 .2 1.0
.175 0.5 .175 1.0
.15 0.5 .15 1.0
.125 0.5 .125 1.0
.1 0.5 .1 1.0
.075 0.5 .075 1.0
.05 0.5 .05 1.0
.025 0.5 .025 1.0
0 0.5 0 1.0
!Wavevector list number 2 (Cartesian directions for non-analytic gamma phonons)
!The output for this calculation must be cut-and-pasted into the
! t59_out.freq file to be used as band2eps input to get proper LO-TO
! splitting at gamma. Note that gamma occurrs twice.
nph2l 1 ! number of directions in list 2
qph2l 1.0 0.0 0.0 0.0
# This line added when defaults were changed (v5.3) to keep the previous, old behaviour
symdynmat 0
prtdos 1
ng2qpt 40 40 40
I am attempting to calculate the phonon band structure of cubic SrTiO3 with DFPT. Using the PAW potentials generated by AtomPAW2Abinit and the Perdew-Wang functional I obtain the erroneous result show in the attached figure phonon_band_structure_PAW.png. Although I have carried out the correction for the TO-LO splitting, the dispersion is still not smooth at Gamma. Moreover, the phonon band structure calculated using the PAW potentials do not show the negative frequencies in the Gamma and M points due to the cubic phase not being the ground state, for instance reported in Ph. Ghosez et al., AIP Conf. Proc. 535 (1), 102 (2000). Since the latter work was carried out using norm-conserving PSP I've redone the calculation using OPIUM generated PSP and the PBE functional. The result is shown in phonon_band_struct_PSP.png. In this case the LO-TO splitting is correctly accounted for, and the band structure is generally in good agreement with that reported in Ghosez et al. Yet, there is an extra instability around X in my calculation, which the reference shows for BaTiO3 and Ba0.5Sr0.5TiO3, but not for SrTiO3. I'm guessing this is due to the use of PBE functional in my calculation which introduces a relative error of +1% on the lattice parameter compared with the use of the LDA in the reference, which I found produces an error of -1% on the lattice parameters in case of SrTiO3.
Does the implementation of DFPT with PAWs still have some issues in the current version Abinit 7.8.2, in particular with regard to the TO-LO splitting? I hope someone can clarify these findings.
My script for the generation of the derivative database using PAWs (based on the input file from the second tutorial on response calculations) and the input file for the analysis of the ddb, to obtain the phonon band structure are displayed below.
Best wishes,
Jonas
******************************************************************
******************************************************************
# Crystalline SrTiO3 : computation of the phonon spectrum
ndtset 22
#Set 1 : ground state self-consistency
getwfk1 0 # Cancel default
kptopt1 1 # Automatic generation of k points, taking
# into account the symmetry
nqpt1 0 # Cancel default
tolwfr1 1.0d-22 # SCF stopping criterion (modify default)
rfphon1 0 # Cancel default
iscf1 17
nstep1 1600
prtwf1 1
#Q vectors for all datasets
#Complete set of symmetry-inequivalent qpt chosen to be commensurate
# with kpt mesh so that only one set of GS wave functions is needed.
#Generated automatically by running GS calculation with kptopt=1,
# nshift=0, shiftk=0 0 0 (to include gamma) and taking output kpt set
# file as qpt set. Set nstep=1 so only one iteration runs.
nqpt 1 # One qpt for each dataset (only 0 or 1 allowed)
# This is the default for all datasets and must
# be explicitly turned off for dataset 1.
qpt2 0.00000000E+00 0.00000000E+00 0.00000000E+00
qpt3 0.00000000E+00 0.00000000E+00 0.00000000E+00
qpt4 1.66666667E-01 0.00000000E+00 0.00000000E+00
qpt5 3.33333333E-01 0.00000000E+00 0.00000000E+00
qpt6 5.00000000E-01 0.00000000E+00 0.00000000E+00
qpt7 1.66666667E-01 1.66666667E-01 0.00000000E+00
qpt8 3.33333333E-01 1.66666667E-01 0.00000000E+00
qpt9 5.00000000E-01 1.66666667E-01 0.00000000E+00
qpt10 3.33333333E-01 3.33333333E-01 0.00000000E+00
qpt11 5.00000000E-01 3.33333333E-01 0.00000000E+00
qpt12 5.00000000E-01 5.00000000E-01 0.00000000E+00
qpt13 1.66666667E-01 1.66666667E-01 1.66666667E-01
qpt14 3.33333333E-01 1.66666667E-01 1.66666667E-01
qpt15 5.00000000E-01 1.66666667E-01 1.66666667E-01
qpt16 3.33333333E-01 3.33333333E-01 1.66666667E-01
qpt17 5.00000000E-01 3.33333333E-01 1.66666667E-01
qpt18 5.00000000E-01 5.00000000E-01 1.66666667E-01
qpt19 3.33333333E-01 3.33333333E-01 3.33333333E-01
qpt20 5.00000000E-01 3.33333333E-01 3.33333333E-01
qpt21 5.00000000E-01 5.00000000E-01 3.33333333E-01
qpt22 5.00000000E-01 5.00000000E-01 5.00000000E-01
#Set 2 : Response function calculation of d/dk wave function
iscf2 -3 # Need this non-self-consistent option for d/dk
kptopt2 2 # Modify default to use time-reversal symmetry
rfphon2 0 # Cancel default
rfelfd2 2 # Calculate d/dk wave function only
tolwfr2 1.0d-22 # Use wave function residual criterion instead
nstep2 1600
prtwf2 1
#Set 3 : Response function calculation of Q=0 phonons and electric field pert.
getddk3 2 # d/dk wave functions from last dataset
kptopt3 2 # Modify default to use time-reversal symmetry
rfelfd3 3 # Electric-field perturbation response only
#Sets 4-22 : Finite-wave-vector phonon calculations (defaults for all datasets)
getwfk 1 # Use GS wave functions from dataset1
kptopt 3 # Need full k-point set for finite-Q response
rfphon 1 # Do phonon response
rfatpol 1 5 # Treat displacements of all atoms
rfdir 1 1 1 # Do all directions (symmetry will be used)
tolvrs 1.0d-8 # This default is active for sets 3-22
prtwf 0
#######################################################################
#Common input variables
# Definition of the unit cell
#acell 3*7.27888848134371
acell 3*7.27898250385070
rprim # a cubic unit cell
1 0 0
0 1 0
0 0 1
# Definition of the atom types
ntypat 3
znucl 38 22 8 # strontium (38), titanium (22), oxygen (8)
# Definition of the atoms
natom 5
typat 1 2 3 3 3
xred
0.0 0.0 0.0 # Sr sits at cubic site
0.5 0.5 0.5 # Ti sits at body-center site
0.0 0.5 0.5 # O sits at face-centered sites
0.5 0.0 0.5
0.5 0.5 0.0
#Gives the number of band, explicitely (do not take the default)
nband 20
#Exchange-correlation functional
ixc 7 # LDA-Perdew-Wang
#Definition of the planewave basis set
ecut 60 # Maximal kinetic energy cut-off, in Hartree
pawecutdg 80
#Definition of the k-point grid
ngkpt 6 6 6
nshiftk 1 # Use one copy of grid only (default)
shiftk 0.5 0.5 0.5
paral_kgb 0
#Definition of the SCF procedure
iscf 7 # Self-consistent calculation, using algorithm 5
nstep 400 # Maximal number of SCF cycles
diemac 6.0 # Although this is not mandatory, it is worth to
# precondition the SCF cycle. The model dielectric
# function used as the standard preconditioner
# is described in the "dielng" input variable section.
# The dielectric constant of AlAs is smaller that the one of Si (=12).
# add to conserve old < 6.7.2 behavior for calculating forces at each SCF step
optforces 1
******************************************************************
******************************************************************
!Input file for the anaddb code. Analysis of the SrTiO3 DDB
!Flags
ifcflag 1 ! Interatomic force constant flag
!Wavevector grid number 1 (coarse grid, from DDB)
brav 1 ! Bravais Lattice : 1-S.C., 2-F.C., 3-B.C., 4-Hex.)
ngqpt 6 6 6 ! Monkhorst-Pack indices
nqshft 1 ! number of q-points in repeated basic q-cell
q1shft 3*0.0
!Effective charges
asr 1 ! Acoustic Sum Rule. 1 => imposed asymetrically
chneut 1 ! Charge neutrality requirement for effective charges.
!Interatomic force constant info
dipdip 1 ! Dipole-dipole interaction treatment
!Phonon band structure output for band2eps - See note near end for
! dealing with gamma LO-TO splitting issue.
eivec 4
!Wavevector list number 1 (Reduced coordinates and normalization factor)
nph1l 101 ! number of phonons in list 1
qph1l 0 0 0 1.0
0 0.025 0 1.0
0 0.05 0 1.0
0 0.075 0 1.0
0 0.1 0 1.0
0 0.125 0 1.0
0 0.15 0 1.0
0 0.175 0 1.0
0 0.2 0 1.0
0 0.225 0 1.0
0 0.25 0 1.0
0 0.275 0 1.0
0 0.3 0 1.0
0 0.325 0 1.0
0 0.35 0 1.0
0 0.375 0 1.0
0 0.4 0 1.0
0 0.425 0 1.0
0 0.45 0 1.0
0 0.475 0 1.0
0 0.5 0 1.0
0.025 0.5 0 1.0
0.05 0.5 0 1.0
0.075 0.5 0 1.0
0.1 0.5 0 1.0
0.125 0.5 0 1.0
0.15 0.5 0 1.0
0.175 0.5 0 1.0
0.2 0.5 0 1.0
0.225 0.5 0 1.0
0.25 0.5 0 1.0
0.275 0.5 0 1.0
0.3 0.5 0 1.0
0.325 0.5 0 1.0
0.35 0.5 0 1.0
0.375 0.5 0 1.0
0.4 0.5 0 1.0
0.425 0.5 0 1.0
0.45 0.5 0 1.0
0.475 0.5 0 1.0
0.5 0.5 0 1.0
.475 .475 0 1.0
.45 .45 0 1.0
.425 .425 0 1.0
.4 .4 0 1.0
.375 .375 0 1.0
.35 .35 0 1.0
.325 .325 0 1.0
.3 .3 0 1.0
.275 .275 0 1.0
.25 .25 0 1.0
.225 .225 0 1.0
.2 .2 0 1.0
.175 .175 0 1.0
.15 .15 0 1.0
.125 .125 0 1.0
.1 .1 0 1.0
.075 .075 0 1.0
.05 .05 0 1.0
.025 .025 0 1.0
0 0 0 1.0
0.025 0.025 0.025 1.0
0.05 0.05 0.05 1.0
0.075 0.075 0.075 1.0
0.1 0.1 0.1 1.0
0.125 0.125 0.125 1.0
0.15 0.15 0.15 1.0
0.175 0.175 0.175 1.0
0.2 0.2 0.2 1.0
0.225 0.225 0.225 1.0
0.25 0.25 0.25 1.0
0.275 0.275 0.275 1.0
0.3 0.3 0.3 1.0
0.325 0.325 0.325 1.0
0.35 0.35 0.35 1.0
0.375 0.375 0.375 1.0
0.4 0.4 0.4 1.0
0.425 0.425 0.425 1.0
0.45 0.45 0.45 1.0
0.475 0.475 0.475 1.0
0.5 0.5 0.5 1.0
.475 0.5 .475 1.0
.45 0.5 .45 1.0
.425 0.5 .425 1.0
.4 0.5 .4 1.0
.375 0.5 .375 1.0
.35 0.5 .35 1.0
.325 0.5 .325 1.0
.3 0.5 .3 1.0
.275 0.5 .275 1.0
.25 0.5 .25 1.0
.225 0.5 .225 1.0
.2 0.5 .2 1.0
.175 0.5 .175 1.0
.15 0.5 .15 1.0
.125 0.5 .125 1.0
.1 0.5 .1 1.0
.075 0.5 .075 1.0
.05 0.5 .05 1.0
.025 0.5 .025 1.0
0 0.5 0 1.0
!Wavevector list number 2 (Cartesian directions for non-analytic gamma phonons)
!The output for this calculation must be cut-and-pasted into the
! t59_out.freq file to be used as band2eps input to get proper LO-TO
! splitting at gamma. Note that gamma occurrs twice.
nph2l 1 ! number of directions in list 2
qph2l 1.0 0.0 0.0 0.0
# This line added when defaults were changed (v5.3) to keep the previous, old behaviour
symdynmat 0
prtdos 1
ng2qpt 40 40 40