ABINIT Discussion Forums

Hello,

I have been doing some calculations on a system of NaH but I am having problems getting the LO-TO splitting that is typical of a polar system like this. This system has been studied before and it is known to have an LO-TO splitting at the gamma point if the effective charges are included in the calculation (Ke and Tanaka PRB 71 024117 (2005)).

I have ran calculations before on other system and been able to get LO-TO splitting, for instance on ZnTe (a simple zinc-blende structure) and the tutorial structure AlAs.

Here is my input file for the calculation of the DDB files (Sorry about the size):


****************************************************************************************************************
****************************************************************************************************************
# Crystalline NaH : computation of the phonon spectrum

ndtset 10
#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
tolvrs1 1.0d-21 # SCF stopping criterion (modify default)
rfphon1 0 # Cancel default
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 2.50000000E-01 0.00000000E+00 0.00000000E+00
qpt5 5.00000000E-01 0.00000000E+00 0.00000000E+00
qpt6 2.50000000E-01 2.50000000E-01 0.00000000E+00
qpt7 5.00000000E-01 2.50000000E-01 0.00000000E+00
qpt8 -2.50000000E-01 2.50000000E-01 0.00000000E+00
qpt9 5.00000000E-01 5.00000000E-01 0.00000000E+00
qpt10 -2.50000000E-01 5.00000000E-01 2.50000000E-01

#nqpt 3-nkpt+2 are copied from make_qpt from kpt in
#the make_qpt.out file.
#insert nqpt 2 = gamma because it is the electrical part
#of the phonon around the gamma point. Think forces
#due to charges moving

#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-25 # Use wave function residual criterion instead

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-11 : 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 2 # Treat displacements of all atoms
rfdir 1 1 1 # Do all directions (symmetry will be used)
tolvrs 1.0d-11 # This default is active for sets 3-10

prtden 0
prtwf 0

#######################################################################
#Common input variables

#Definition of the unit cell
acell 6.4410481435E+00 6.4410481435E+00 6.4410481435E+00 Bohr

rprim 7.0710678119E-01 7.0710678119E-01 0.0000000000E+00
0.0000000000E+00 7.0710678119E-01 7.0710678119E-01
7.0710678119E-01 0.0000000000E+00 7.0710678119E-01

xred 5.0000000000E-01 5.0000000000E-01 5.0000000000E-01
0.0000000000E+00 0.0000000000E+00 0.0000000000E+00

#Definition of the atom types
ntypat 2 # There are two types of atom
znucl 1 11 # The keyword "znucl" refers to the atomic number of the
# possible type(s) of atom. The pseudopotential(s)
# mentioned in the "files" file must correspond
# to the type(s) of atom. Here, type 1 is the Aluminum,
# type 2 is the Arsenic.

#Definition of the atoms
natom 2 # There are two atoms
typat 1 2 # The first is of type 1 (Al), the second is of type 2 (As).

#Gives the number of band, explicitely (do not take the default)
nband 10

#Definition of the planewave basis set

ecut 40 Ha # Maximal kinetic energy cut-off, in Hartree

#Definition of the k-point grid #DEFINED BY findk DIRECTORY
ngkpt 8 8 8

nshiftk 1

shiftk 3*0.0

#Definition of the SCF procedure
# iscf 7 # Self-consistent calculation, using algorithm 5
nstep 3000 # Maximal number of SCF cycles
diemac 12.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.

occopt 7
tsmear 0.01 eV

*********************************************************************************************
****************************************************************************************************************

I merge the DDB's out of this run using: mrgddb < merge.in

***************************************************************************
****************************************************************************************************************
merge.ddb.out
NaH phonons on 4 4 4 mesh
8
./phonstat/phonstato_DS3_DDB
./phonstat/phonstato_DS4_DDB
./phonstat/phonstato_DS5_DDB
./phonstat/phonstato_DS6_DDB
./phonstat/phonstato_DS7_DDB
./phonstat/phonstato_DS8_DDB
./phonstat/phonstato_DS9_DDB
./phonstat/phonstato_DS10_DDB
**************************************************************************
****************************************************************************************************************

Then I run a phonon band-structure calculation using anaddb on the merge.ddb.out calculated above.

where the input file for the anaddb calculation on the merged DDB is:

****************************************************************************************************************
*****************************************************************************************************************
!Input file for the anaddb code. Analysis of the NaH 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 4 4 4 ! 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

# This line added when defaults were changed (v5.3) to keep the previous, old behaviour
symdynmat 1

!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

!Wavevector list number 1 (Reduced coordinates and normalization factor)
nph1l 1001 ! number of phonons in list 1

qph1l
0.5000 0.5000 0.5000 1.0
0.4975 0.4975 0.4975 1.0
0.4950 0.4950 0.4950 1.0
0.4925 0.4925 0.4925 1.0
0.4900 0.4900 0.4900 1.0
.
.
.
*****************************************************************************************
*****************************************************************************************

From what I can tell I am not doing anything that would not give me a correct LO-TO splitting at Gamma. But unfortunately I get none. All of my calculations are converged in the initial run that calculates the DDB's (there are 0 WARNING's and 0 COMMENT's in the log file).

However, when I look at the output of this run that generates the original DDB's I get the following which tells me there is a clear splitting seen around the Gamma point.
-----------------------------------------------------------------------------------------------------------------------------------------------------------
Phonon wavevector (reduced coordinates) : 0.00000 0.00000 0.00000
Phonon energies in Hartree :
2.711776E-05 2.712757E-05 2.716443E-05 2.174407E-03 2.174407E-03
2.174407E-03
Phonon frequencies in cm-1 :
- 5.951661E+00 5.953813E+00 5.961904E+00 4.772271E+02 4.772271E+02
- 4.772271E+02

Phonon at Gamma, with non-analyticity in the
direction (cartesian coordinates) 1.00000 0.00000 0.00000
Phonon energies in Hartree :
2.711783E-05 2.712938E-05 2.716541E-05 2.174407E-03 2.174407E-03
4.041351E-03
Phonon frequencies in cm-1 :
- 5.951675E+00 5.954210E+00 5.962119E+00 4.772271E+02 4.772271E+02
- 8.869739E+02

Phonon at Gamma, with non-analyticity in the
direction (cartesian coordinates) 0.00000 1.00000 0.00000
Phonon energies in Hartree :
2.711850E-05 2.712851E-05 2.716561E-05 2.174407E-03 2.174407E-03
4.041351E-03
Phonon frequencies in cm-1 :
- 5.951823E+00 5.954019E+00 5.962163E+00 4.772271E+02 4.772271E+02
- 8.869741E+02

Phonon at Gamma, with non-analyticity in the
direction (cartesian coordinates) 0.00000 0.00000 1.00000
Phonon energies in Hartree :
2.711965E-05 2.712771E-05 2.716525E-05 2.174407E-03 2.174407E-03
4.041350E-03
Phonon frequencies in cm-1 :
- 5.952075E+00 5.953845E+00 5.962084E+00 4.772271E+02 4.772271E+02
- 8.869739E+02
---------------------------------------------------------------------------------------------------------------------------------------------------------

However, at the end of the day when I am looking at the output of anaddb where it has analyzed the merged DDB. I get:
-------------------------------------------------------------------------------------------------------------------------------------------------------
Phonon at Gamma, with non-analyticity in the
direction (cartesian coordinates) 1.00000 0.00000 0.00000
Phonon energies in Hartree :
0.000000E+00 0.000000E+00 0.000000E+00 2.173016E-03 2.173016E-03
2.173017E-03
Phonon frequencies in cm-1 :
- 0.000000E+00 0.000000E+00 0.000000E+00 4.769218E+02 4.769219E+02
- 4.769222E+02
-------------------------------------------------------------------------------------------------------------------------

Which clearly tells me that anaddb does not see any phonon splitting (which is what I see in my plot, attached).


I've attached my abinit input file and my anaddb input file for viewing. and

I simply do not understand why my final answer out of anaddb does not have LO-TO splitting at the Gamma point. I would greatly appreciate any advice on this!

Thank you!

Gamma is an independent calculation where you get directly the LO/TO splitting. You hacve to activate your d/dk calculation and then the response to the electric field
and you should be able to get the LO/TO... with no problem.

regards

-aldo.

Hello,

After playing with the multiple flags in my input files for both the anaddb and abinit calculations I have found that changing occopt from 7 to 1 gives the correct LO/TO splitting at the gamma point. This seems strange to me. Is anaddb assuming I have a metal simply because I use occopt 7 instead of 1? Most of my calculations converge faster if occopt is set to 7 (gaussian smearing) than if it is set to 1 (assumed semi-conductor (all occupations automatically stop at the top of the valence band)).

Oh well, at least I figured it out. I hope this helps other people who have this same problem in the future.

@ aldo. Hello and thank you for replying, it is appreciated. You are correct that the calculations of the phonon's at the gamma point are a seperate calculation from those that are not at the gamma point, and that the d/dk calculation must first be ran at the gamma point. However, both calculations, at the gamma point and those off the gamma point, can use the same initial self consistent calculation's wavefunctions (if they are set up correctly), which is what I have done here.

[quote="Redorb"]This seems strange to me. Is anaddb assuming I have a metal simply because I use occopt 7 instead of 1? [\quote]
yes - this is no longer the case in the next version of abinit - now if you ask for the dielectric stuff and if the perturbations wrt E-field are present things are calculated. Normally if you have metallic occupation there is no dielectric response.


Matthieu