Comment on "First-principles study of the influence of (110)-oriented strain on the ferroelectric properties of rutile TiO2" [arXiv:1106.2820]

0.13 MB·English

by  Keith Refson

Comment on “First-principles study of the influence of (110)-oriented strain on the
ferroelectric properties of rutile TiO ”
2
∗
Keith Refson and Barbara Montanari
STFC, Rutherford Appleton Laboratory, Didcot, Oxfordshire, OX11 0QX, UK
Pavlin D. Mitev and Kersti Hermansson
Materials Chemistry, ˚Angstr¨om Laboratory, P.O. Box 538, SE-752 21 Uppsala, Sweden
3 Nicholas M. Harrison
1 Department of Chemistry,Imperial College London,
0 South Kensington Campus, London, SW7 2AZ, UK†
2
In a recent article, Gru¨nebohm et al. [Phys. Rev. B 84 132105 (2011)] report that they fail to
n
a reproduce the A2u ferroelectric instability of TiO2 in the rutile structure calculated with density
J functional theory within the PBE-GGA approximation by Montanari et al. [Chem. Phys. Lett
364, 528 (2002)]. We demonstrate that this disagreement arises from an erroneous treatment of
3
Ti 3s and 3p semi-core electrons as core in their calculations. Fortuitously the effect of the frozen
2
semi-core pseudopotential cancels the phonon instability of thePBE exchange-correlation, and the
] combination yields phonon frequencies similar to the LDAharmonic values.
i Gru¨nebohmetal. alsoattemptedandfailedtoreproducethesoftacousticphononmodeinstability
c
s under (110) strain reported by Mitev et al. [Phys. Rev. B 81 134303 (2010)]. For this mode the
- combination of PBE-GGA and frozen semi-core yields a small imaginary frequency of 9.8i. The
l
r failure of Gru¨nebohm et al. to find this mode probably arose from numerical limitations of the
t geometryoptimizationapproachinthepresenceofashallowdoublewellpotential;theoptimization
m
method is not suitable for locating such instabilities.
.
t
a PACSnumbers: 63.20.D-,71.15.Mb,77.84.-s,77.22.-d,77.80.bn
m
-
d I. INTRODUCTION GGA does not predict a phonon instability of the rutile
n structure,thatthereis nosofteningofthe A mode un-
2u
o
der(110)strain,noristhereadestabilizationunder(110)
RutileTiO isanincipientferroelectricmaterial. Over
c 2
strain at (1/2,1/2,1/4)by a soft acoustic mode.
[ a 20 yearperioddensity functional theory has beenused
todevelopanunderstandingofitsdielectricresponseand It would be highly unsatisfactory for this discrepancy
1
lattice dynamical properties. It has also established the to remain unresolved in the literature; it would weaken
v
4 exquisite sensitivity ofthese properties to strain,volume the basis for future studies on this topic, and under-
6 changes and thus to approximations such as the treat- mineconfidenceinthereproducibilityofDFTplane-wave
5 ment of electronic exchange and correlation1. methods forthis type ofstudy. Gru¨nebohmet al. offera
5 In 2002 two of us demonstrated2 that within the speculation that the discrepancy arises from “the differ-
1. PerdewBurkeandEnzerhofgeneralizedgradientapprox- ent potentials used and the resulting difference in lattice
0 imation (PBE-GGA)3 the A mode is unstable at the parameters”. We tested that hypothesis, and confirm it
2u
3 equilibrium geometry. The resulting structural instabil- as the most likely origin of the discrepancy. We show
1 itycorrespondstoaferroelectricdistortionofthe crystal that the error lies in the choice of pseudopotentials used
v: in qualitative disagreementwith the observedproperties by Gru¨nebohm et al. specifically their treatment of the
i of rutile TiO 4. In the local density approximation Ti 3s and 3p semicore electrons as part of the frozen
X 2
(LDA)5 the equilibrium cell volume is slightly smaller pseudopotential core instead of as valence.
r andastablestructureisobtained,givingaquantitatively Textbook wisdom holds that in the case of titanate
a
correctdescription of the lattice dynamics and dielectric perovskitesthefrozensemicoreapproximationintroduces
properties. This motivated the choice of LDA functional significant errors in ferroelectric properties and that 3p
for our study of Ref. 6 on the influence of strain on the and possibly 3s states must be explicitly included8. In-
ferroelectricpropertiesandphononmodesofrutileTiO . deed most pseudopotential studies of TiO and titanate
2 2
The major conclusion of that work was the existence of perovskites published since the mid 1990s have treated
a low frequency acoustic phonon branch over a broad the Ti 3s and 3p electrons as valence, claiming that the
region of the Brillouin Zone around (1/2,1/2,1/4)which frozen semicore approximation is inaccurate. However
becamesoftundernegativepressureoranisotropicstrain. thereisnoreportcontainingquantitativeevidenceofthe
In article Ref. 7, Gru¨nebohm, Ederer,and Entel, pub- consequences for phonon frequencies.
lished another study on the same topic, in which they Ghosez et al.13,14 performed band-by band decompo-
attempted and failed to reproduce those central findings sition of the Born effective charges and showed that the
of our work. Specifically they report that that the PBE- Ti3porbitalscontributeavalueof-0.22etotheeffective
2
a large error in the atomic configuration energy for s-d
TABLE I. Pseudopotential configurations used. A single ul-
promotion in Ti atoms of 90 meV obtained within the
trasoft projector was used for each core or semi-core state.
frozen-semicore approximation. However Perron et al.19
Two ultrasoft projectors were used for each valence state ex-
showed that the effect of a frozen-core pseudopotential
cept for UF1where a single projector was used for Ti 4s and
on structural parameters of TiO rutile and anatase is
UF3 where a single norm-conserving projector was used for 2
Ti 4s. For all frozen core and semi-core states a pseudized less than 0.5% provided that a nonlinear core correction
core charge9 was used. term is included20.
The question of how many electrons to treat as va-
Label UF1 UF2 UF3
lence is independent of other aspects of the pseudopo-
Ti ref. 2s22p63s23p64s23d2 3s23p64s23d2 4s23d2
tential construction or formalism, and holds equally for
r 1.2/1.2/1.8 1.8 2.2
c norm-conserving, ultrasoft or PAW pseudopotentials18,
r /a 0.5 1.5 1.3
inner 0 as it concerns the inclusion or omission of the physics of
l /a d d s
loc 0 polarizability of the Ti 3s and 3p electrons.
O ref. 1s22s22p4 2s22p4 2s22p4
In summary, previous investigations are consistent
r /a 0.8 1.3 1.3
c 0 with a picture in which structural parameters are little
r /a 0.5 0.9 0.9
inner 0 affected by the treatment of 3s and 3p, relative energies
Cutoff 1300 750 750
of phases and cohesive energies are noticeably affected.
/eV
Frozensemicorepseudopotentialshavepoortransferabil-
ity resulting in significant errors in dielectric properties
andchemicalreactionenergetics. Asnoresultshavebeen
TABLE II. Calculated structural parameters of rutile TiO2 reported on the consequence for phonon frequencies and
mechanical stability, we present them here.
a/A c/A u
PW-LDA (UF1) 4.550 2.919 0.3039
PW-LDA (UF2) 4.549 2.919 0.3039
PW-LDA (UF3) 4.597 2.902 0.3023 II. GGA INSTABILITY OF RUTILE
LCAO-LDAa 4.548 2.944 0.305 STRUCTURE
FP-LAPW-LDAb 4.558 2.920 0.3039
PW-GGA (UF1) 4.642 2.963 0.3051 We performed calculations on bulk rutile TiO using
2
PW-GGA (UF2) 4.642 2.964 0.3051 theCASTEPcode21,withaseriesofpseudopotentialsto
PW-GGA (UF3) 4.684 2.953 0.3035 testtheeffectofthefrozen-coreapproximationonsucces-
PW-GGAc 4.641 2.966 0.305 siveshells. Thepseudopotentialswereoftheultrasoftva-
PW-GGAd 4.664 2.969 0.3047 riety generated using Vanderbilt’s method22. Small core
LCAO-GGAe 4.623 2.987 0.306 radii were chosen to generate accurate but “hard” pseu-
dopotentials (i.e. which require a large cutoff) for both
a Ref.1and10 elements to leave the frozen-core approximation as the
b Ref.11 dominant contribution to the error. Details of these po-
c Ref.2
tentials are listed in table I. Three sets of pseudopoten-
d Ref.7
e Ref.12 tials are designated UF1, UF2, and UF3 denoting that
theTi1s,2sand2p,and3sand3pshellswerefrozenre-
spectively. The UF3 set corresponds to the “large core”
charge on the Ti atoms in SrTiO and -0.43e in the case set used by Gru¨nebohm et al. and the UF2 set corre-
3
of BaTiO perovskites. The same authors also argue15 sponds to the “small core” set used by Montanari and
3
that small contributions to this effective charge are re- Harrison and by previous authors. The UF1 set gives a
sponsible for the destabilization of the cubic phase by “nearly all electron” calculation in which only the Ti 1s
the softening of an optic phonon and a phase transition orbitals are in the core, and Ti 2s and 2p and O 1s are
to the rhombohedral phase. Tegner et al.16 performed explicitly treated as valence. The computational cost of
a comprehensive investigation on the effect of inclusion treating core states as valence is surprisingly reasonable
of Ti 3s and 3p on cohesive energies, phase stability and because the sharply-peakedcore orbitals are represented
electronicdensityofstatesofmetallicTi,concludingthat mostly by the pseudized augmentation functions of the
frozensemicoreTipseudopotentialsresultinpoortrans- Vanderbilt scheme. A very fine FFT grid is required
ferability andthatthe relativestability ofthe ω andhcp only for the electron density to represent the augmenta-
phases is reversed. Deskins17 found an 0.3 eV difference tion charge, not for the soft orbital functions. This “all
in the cohesive energy of TiO and again that semicore electron pseudopotential” approach has previously used
2
Ti pseudopotentials exhibit worse transferability in the successfully for calculations on iron at TPa pressures23.
case of changes of oxidation state, with errors of up to In each case a geometry optimization was performed
0.1 eV in dissociationreaction energies of TiCl and ab- to optimize both lattice parameters and the internal co-
n
sorption energetics of a variety of molecules and atoms ordinate, followed by a finite displacement calculation
on rutile TiO surfaces. Holzwarth et al.18 highlighted of the dynamical matrix. This was diagonalized to yield
2
3
TABLE III. Calculated and experimental Γ-point phonon frequencies of bulkrutile TiO (cm−1). (TO frequencies only.)
2
Mode LDA-UF1 LDA-UF2 LDAa LDA-LAPWb PBE-UF1 PBE-UF2 PBEa PBE-UF3c PBE-UF3 Expta
B 99.2 101.7 104.0 111.5 31.0i 28.4i 79.2 96.1 91.2 113
2u
A 125.1 128.7 154.4 153.0 89.6i 106.6i 86.3i 143.7 115.7 142
2u
E 135.2 138.8 191.4 156.4 69.8i 62.5i 124.0 137.7 111.2 189
u
B 136.9 137.1 137.0 138.2 127.1 150.5 154.2 131.2 144.4 142
2g
E 381.5 382.4 383.9 384.8 353.8 355.1 353.5 383.2 366.4 388
u
B 390.9 392.4 393.0 398.2 362.3 356.4 357.5 400.5 380.3 406
2u
A 417.6 419.5 421.7 416.0 408.5 417.9 423.6 407.6 405.6
2g
−
E 463.1 463.5 463.2 463.7 429.1 428.9 429.2 472.2 464.3 455
g
E 486.7 487.3 488.4 487.0 469.6 469.7 488.4 480.0 491.4 494
u
A 609.9 610.4 611.6 612.0 567.8 567.3 565.9 626.9 611.1 610
1g
B 817.3 818.2 824.7 816.1 770.8 771.4 774.3 815.2 794.6 825
1g
a Ref.2
b Ref.11. Thesearepseudopotential-LAPW, notall-electron
FP-LAPWcalculations
c CalculatedusingtheoptimizedlatticeparametersofUF2
thegamma-pointphononfrequencies. TheBrillouinzone
was sampled using a Monkhorst Pack24 grid of dimen- TABLE IV. Frequencies of the lowest mode at
q=(0.5,0.5,0.25) calculated using using a √2 √2 4
sions 4×4×6 and plane-wave cutoffs for each USP set supercell (cm−1). × ×
are listed in table I. These choicesyield a convergenceof
betterthan1meV/atomintotalenergy,0.00025Ainlat- UF1 UF2 UF3
ticeparameter. Thegeometrywasrelaxeduntiltheforce PBE 77.1i 72.6i 9.8i
residualoneachatomwaslessthan0.005eV/A.Formost LDA 50.0 55.1 74.8
phonons the frequencies were converged to better than
0.1 cm−1; for the most sensitive phonon modes (which
become imaginary) the maximum error was 3 cm−1.
UF2 configuration contains all of the electronic freedom
The calculated structural parameters are shown in ta-
requiredforaccuratelatticedynamics. Thisconfirmsthe
ble II. LDA lattice parameters for UF1 and UF2 are
correctnessofthe“small-core”pseudopotentialresultsof
within 0.2% of those calculated using the reference all-
Ref.2whereitwasdemonstratedthatall-electronLCAO
electron techniques1,11. Freezing the Ti semicore states
methodsgavethesameA instabilityastheplane-wave
slightly increases the lattice parameters but by less than 2u
pseudopotential calculations. We note one discrepancy
1%, consistent with the values reported by Gru¨nebohm
withthe priorresults ofRef. 2, namely thatthe B and
et al. and by previous studies19,20. 2u
E modes also exhibit an instability in addition to the
u
Phonon frequencies are presented in table III. This
A instability previously described.
2u
shows clearly that (a) the rutile structure is predicted
to suffer from a phonon instability within PBE and (b)
that freezing the semicore states has the effect of elimi-
nating the instability and yields all real frequencies. As III. MODE SOFTENING AT Q=(1/2,1/2,1/4)
theatomicreferenceconfigurationisthesameasusedby
Gru¨nebohmetal. itishighlylikelythatthisexplainsthe Turning to the softacousticmode frequencies,we per-
failure to reproduce the instability in their calculation. formedsupercellphononcalculationsusinga√2 √2 4
× ×
The change in frequencies when the UF3 pseudopoten- supercellanddiagonalizedatthecommensuratewavevec-
tialsareusedatthe lattice parametersofthe UF2 calcu- tor q=(0.5,0.5,0.25). Frequencies of the lowest acoustic
lation is much smaller than the UF2-UF3 change. This modes are presented in table IV. Clearly the UF3 pseu-
provesthatthediscrepancybetweenGru¨nebohm’sresult dopotentialmakesasubstantialerrorinthismode,which
and ours is a direct consequence of the omission of elec- isstiffenedwhenthe effectofsemicorepolarizationisex-
tronic polarizability and hybridisation of the Ti 3s and cluded. As with the Γ -point phonons, there is little dif-
3p states and not an indirect consequence of the volume ference between the unfrozen core UF2 compared with
change. the nearly all-electron UF1 frequencies. The UF1 and
The accuracy of the pseudopotentials is confirmed by UF2 LDA results are in goodagreementwith our earlier
the LDA results which are in excellent agreement with calculations in Ref. 6. However our current calculations
previous work including FP-LAPW calculations. The predictthatwithinthePBEapproximationthisacoustic
similarity of the UF2 results to the nearly all-electron modehasimaginaryfrequencyforallthreepseudopoten-
UF1valuesforbothLDAandPBEdemonstratesthatthe tialconfigurationsconsidered. AswiththeGammapoint
4
results,wewouldexpecttheUF3casetomostcloselyre- search for the soft displacement using this method was
produce the results of Ref. 7. unsuccessful. This is absence of evidence for the soft
Gru¨nebohm et al. state that “we could not repro- mode, not evidence of absence.
duce the destabilization of the system along an acous- Inourexperiencetheonlyreliablewaytoexploreasoft
tic phonon mode under (110) strain found in Ref. 825, mode using geometry optimization is to first perform a
for which an energy gain of several meV has been pre- phonon calculation, and to identify the soft eigenvector
dicted, evenif we use a 2.√2 2.√2 4 supercell,which by diagonalization of a dynamical matrix. This eigen-
is commensurable with the d×isplacem×ent pattern of this vector can be used as an initial perturbation from the
mode”. symmetric supercell. The geometry optimization must
use very tight force and energy tolerances,but is usually
However the geometry optimization method as em-
abletoproceedwithoutinterferencefromthestiffmodes.
ployed by Gru¨nebohm et al. is unsuitable for finding
the soft displacement pattern. Our experience is that
quasi-Newtontypeoptimizationmethodsusuallyfailany
IV. CONCLUSION
attempt to identify a shallow soft mode, or at best con-
verge extremely slowly. Discovery of a single set of
We have demonstrated that the finding of two of us
atomic displacements which slightly lowers the energy
that rutile TiO is mechanically unstable to a ferroelec-
in a high-dimensional search space of otherwise stiff di- 2
tric A soft mode in the PBE-GGA approximation is
rections is a highly computationally challenging proce- 2u
robust. Gru¨nebohm et al.’s contrary finding can be at-
dure. Severalfactors contribute to the difficulty. First is
tributed to an error in the PAW pseudopotentials they
the high dimensionality of the search space, 576 for the
used, namely the frozen-core approximation applied to
√2 √2 4cellemployed(whichisactually4timeslarger
× × the Ti 3s and 3p electrons. Our convergent series of
thanstrictlynecessarytorepresentthissofteigenvector).
pseudopotentials to a highly accurate “nearly all elec-
Within this space the desired eigenvector is represented
tron”limit demonstratesthat the effectof freezingTi 3s
by a single direction.
and 3p is to stabilise the PBE rutile structure of TiO .
Thesecondfactoristheshapeanddepthoftheanhar- 2
In fact this error almost cancels the PBE instability to
monic potential, which has a well depth of only 2.5 meV
yield Γ point frequencies very similar to the LDA result
(seeFig. 5ofRef.6)andatinygradientalongthemode.
butthecancellationislesscompleteelsewhereintheBril-
Forsmalldisplacementthisyieldsatypicalforceoneach
loun Zone.
Ti of 0.01 eV/A, with a maximum of 0.03 eV/A. Conse-
We further argue that the geometry optimisa-
quently the configuration at the initial points along the
tion methods with convergence tolerances reported by
optimizationwouldbedeemedalreadyconvergedaccord-
Gru¨nebohm et al. are not capable of identifying shallow
ing to the toleranceof0.01eV/Ausedby Gru¨nebohmet soft modes in the Brillouin zone. Our result6 of a soft
al. In our experience a tolerance of 0.001 eV/A or even
mode at q=(1/2,1/2,1/4)is also robustly confirmed.
smaller is essential in such circumstances.
Thus the situation is that in the incipient ferroelectric
The third factor is the difficulty of finding a suitable rutile-TiO the LDA and PBE-GGA approximations to
2
starting configuration for the optimization. In a super- density functional theory yield qualitatively distinct de-
cell of perfect rutile the gradient in every direction is scriptions. Within PBE-GGA, rutile TiO is predicted
2
zero, so an initial perturbation is required to break the to be ferroelectric with an equilibrium geometry that is
supercellsymmetry. Arandomizeddisplacementpattern unstable withrespectto displacementalonga numberof
must also excite the stiff modes and yield a much larger phonon eigenvectors4. Within the LDA, rutile-TiO is
2
increase in the energy than the decrease being sought. stable and a reasonable description of its dielectric re-
The resulting geometry optimization almost always fails sponse and lattice dynamics is obtained.
to find the soft mode displacement. We tested this for
the soft acoustic mode in rutile TiO , using the UF2
2
pseudopotential, a force tolerance of 0.001 eV/A and as ACKNOWLEDGMENTS
expectedtheoptimizationreturnedtothe symmetricsu-
percell configuration. We thank the EPSRC for support, of the UK na-
FromtheUF3phononfrequencyresultsintableIVwe tional supercomputing facility (HECToR) under grant
wouldexpectanevenshallowerpotentialandeventinier EP/F036809/1. Other calculations used STFC’s e-
gradient as a result of the additional stiffening resulting Science facility. K.R. would also like to thank Professor
fromthe neglectofsemi-corepolarization. Itistherefore Chris J. Pickard for introducing him to the concept of
doubly unsurprising that Gru¨nebohm et al.’s attempt to all-electron pseudopotentials.
∗ [email protected] † STFC Daresbury Laboratory, Cheshire, WA4 4AD, UK
5
1 B. Montanari and N. M. Harrison, 12 J. Muscat, N. Harrison, and G. Thornton,
Journal of Physics, Condensed Matter 16, 273 (2004). Physical Review B 59, 2320 (1999).
2 B. Montanari and N. Harrison, 13 P.Ghosez,X.Gonze,P.Lambin, andJ.-P.P.Michenaud,
Chemical Physics Letters 364, 528 (2002). Physical Review B 51, 6765 (1995).
3 J. P. Perdew, K. Burke, and M. Ernzerhof, 14 P. Ghosez, J. P. Michenaud, and X. Gonze,
Physical ReviewLetters 78, 1396 (1997). Physical Review B 58, 6224 (1998).
4 The shallow double well potential obtained using PBE- 15 P. Ghosez, X. Gonze, and J.-P. Michenaud,
GGA yields a prediction of a mechanical instability in Ferroelectrics 194, 39 (1997), arXiv:9706294 [cond-mat].
the classical zero-temperature “athermal” limit. However 16 B. E. Tegner and G. J. Ackland,
this ignores nuclear zero-point quantum effects and ther- Computational Materials Science 52, 2 (2012).
mal motion. It is possible that a fully quantum statistical 17 N.AaronDeskins,Chemical Physics Letters471, 75 (2009).
mechanical PBE calculation which included these effects 18 N.A.W.Holzwarth,G.E.Matthews,R.B.Dunning, and
would predict an undistorted rutile structure even at low A. R.Tackett, Physical Review B 55, 2005 (1997).
temperature. 19 H. Perron, C. Domain, J. Roques,
5 J. P. Perdew and A. Zunger, R. Drot, E. Simoni, and H. Catalette,
Physical ReviewB 23, 5048 (1981). Theoretical Chemistry Accounts 117, 565 (2007).
6 P. D. Mitev, K. Hermansson, B. Montanari, and K. Ref- 20 A much larger discrepancy between small- and large-core
son, Physical Review B 81, 134303 (2010). pseudopotentialsinearlycalculationswasshowntobedue
7 A. Gru¨nebohm, C. Ederer, and P. Entel, to the omission of the “non-linear core correction” term9,
Physical ReviewB 84, 132105 (2011). essential for the correct treatment of exchange and corre-
8 p. 265 in R. M. Martin, lation in thepresence of core-valence overlap in Ti12.
Electronic Structure: Basic Theory and Practical Methods, 21 S. J. Clark, M. D. Segall, C. J. Pickard, P. J. P. J. Has-
1st ed. (Cambridge UniversityPress, 2004). nip, M. I. J. I. J. Probert, K. Refson, and M. C. Payne,
9 S. Louie, S. Froyen, and M. Cohen, Zeitschrift fuerKristallographie 220, 567 (2005).
Physical ReviewB 26, 1738 (1982). 22 D. Vanderbilt,Physical Review B 41, 7892 (1990).
10 J. Muscat, V. Swamy, and N. M. Harrison, 23 C. J. Pickard and R. J. Needs,
Physical ReviewB 65, 224112 (2002). Journal of Physics: Condensed matter 21, 452205 (2009).
11 G.Cangiani,AbinitiostudyofthepropertiesofTiO2rutile 24 H. J. Monkhorst and J. D. Pack,
and anatase polytypes, Phd, Ecole Polytechnique Federale Physical Review B 13, 5188 (1976).
de Lausanne(2003). 25 OurRef. 6.