Authors: Adrien Richard, Furio Corà
Categories: Article
Source: The Journal of Physical Chemistry. C, Nanomaterials and Interfaces
Interactions on the Polymorphic Stability of Crystalline Oxides
The accurate determination of relative phase stabilities
using
DFT methods is a significant challenge when some of these can vary
by only a few kJ/mol. Here, we demonstrate that for a selection of
oxides (TiO2, MnO2, and ZnO) the inclusion of
dispersion interactions, accomplished using the DFT-D3 correction
scheme, allows for the correct ordering and an improved calculation
of the energy differences between polymorphic phases. The energetic
correction provided is of the same order of magnitude as the energy
difference between phases. D3-corrected hybrid functionals systematically
yield results closest to experiment. We propose that the inclusion
of dispersion interactions makes a significant contribution to the
relative energetics of polymorphic phases, especially those with different
densities, and should therefore be included for calculations of relative
energies using DFT methods.
Density functional theory (DFT) has become a regular contributor in modern solid state chemistry, as it gives valuable insight into atomic and electronic properties of solids as well as their functional behavior, complementary to experimental results.^1,2^
The prediction of novel materials and their design requires the calculation of thermodynamic stability to a high level of precision, a need that is accentuated in materials with rich polymorphism, where differences in phase energies can be only a few kJ/mol.^3^ This has been a long-standing issue with DFT, as the choice of functional has a strong influence on the calculated formation energies and relative phase stabilities. No functional to date has produced consistent and reliable phase ordering for materials with different polymorphic forms.^4,5^ Hybrid exchange functionals (HF-DFT functionals), the current state-of-the-art for accuracy, resolve the self-interaction error (SIE) of local DFT functionals, but do not include dispersion forces.^6^ These are typically considered to be negligible in magnitude compared with the total cohesive energy of ceramic materials, dominated by Coulomb and exchange interactions. However, despite dispersion interactions representing only a small fraction of the total binding energy in solids, their energetic contribution is of the same order of magnitude as the difference in cohesive energy between different polymorphic phases.
For reliable determination of these relative phase stabilities,
higher levels of theories can be employed such as second-order Møller–Plesset
perturbation theory (MP2),^7^ Quantum Monte
Carlo (QMC),^8^ Random Phase Approximation
(RPA),^9^ or even Configuration Interaction
(CI).^10^ However, these methods come at
a significantly higher computational cost, which effectively limits
their application to systems possessing only a few atoms. These methods
also show conflicting results, with QMC-based studies predicting the
correct energetic ordering only between certain TiO2 polymorphs,
for example.^8,9^ In recent literature, SCAN (Strongly
Constrained and Appropriately Normed functional) and other meta-GGA
functionals^11,12^ have gained attention in relation
to their successful prediction of relative phase stability in compounds,
such as MnO2.^13^ However, they
are still affected by the SIE. This can be palliated by the inclusion
of empirical parameters, such as the Hubbard U term,
but this naturally leads to a degree of empiricism and poor transferability
of results.^14^
Prediction of the relative
stability of TiO2 polymorphs
is a long-standing issue for DFT calculations; this has been attributed
to DFT’s inability to resolve the SIE arising from systems
with localized d/f-electrons, despite the formal d^0^ electronic configuration of the fully oxidized Ti^4+^ cation,^4,15−19^ which makes the addition of a U term
less justifiable. The order of stability found in the experiment,^20^ i.e., rutile < brookite < anatase is not
reproduced correctly by standard GGA DFT.^21,22^ To resolve this issue, researchers have typically relied on two
main the inclusion of the Hubbard U term
or some post-SCF dispersion correction.^4,18,23−30^
The Hubbard U term^15,31^ introduces
on-site Coulomb and exchange terms in parametric form to account for
noninteger or double occupation of a subset of states such as highly
localized d- and f-electron shells.^32^ A
number of studies justify its use in the calculation of the relative
polymorphic stability of TiO2 through the correction of
the unphysical delocalization of Ti 3d electron states.
Indeed, most studies using DFT+U report an ordering
of polymorph internal energies which matches experimental findings.
However, all studies recognize that the use of DFT+U is only a temporary solution, as not only does the inclusion of
the U term tend to worsen the accuracy of other calculated
properties, but it is also a system-specific solution.^14,33−36^
Over the past decade, DFT methods and published literature have allowed for a quantitative understanding of dispersion interactions, with the DFT-D semiempirical schemes proposed by Grimme being among the most widely applied.^37^ These treatments have shown to accurately describe the structural properties of layered materials and have proven to be especially important for calculations involving molecular crystals or low-dimensional systems, as the London dispersion interaction allows for a correct description of the intermolecular forces involved, thus producing equilibrium structures and energetic values significantly closer to experiment.^38−41^ The DFT-D schemes use damped interatomic potentials corresponding to the instantaneous dipole-induced dipole interaction so that the dispersion corrected energy becomes
where EDFT is
the usual mean-field DFT energy and Edisp is an empirical dispersion correction that is defined in the latest
DFT-D implementation developed for solids at the time of study, DFT-D3,
as the sum of two- and three-body energies, i.e., Edisp = E^2^ + E^3^. The most important two-body term is given by
The first sum spans all atomic pairs
in the system, while Cn^AB^ represents
the averaged
(isotropic) nth-order dispersion coefficient for
atom pair AB, and rAB is their internuclear distance. fd,n(rAB) denotes the damping function applied to the energy correction.
We can easily see from the representation of this sum that as rAB^n^ tends to 0, the entire sum
will tend toward infinity, which will in turn theoretically infinitely
increase the contribution of the two-body term toward the total dispersion
correction energy EDFT-D3. E^3^ has a general expression of the form C6/r^6^ and C9(3 cos θ~a~ cos θ~b~ cos θ~c~ + 1)/(rABrBCrCA)^3^ and is further defined
in ref (40). The magnitude
of the dispersion correction depends on interatomic distances and
hence on the density of the studied materials.
The effect of
dispersion on the formation energy of phases has
been calculated and discussed in previous work,^6,42^ notably
on the phase stabilities of various cesium halides.^43−45^ Similar to
TiO2, the phase stabilities of alkali-metal halides have
been a long-standing problem for DFT and their nature as simple, benchmark
ionic structures spurred researchers to investigate the effects of
dispersion as a computationally cheap solution to the problem. Although
the importance of dispersion has been discussed for clays,^46^ metal–organic frameworks, and zeolites,^42,47^ the main focus there was on the adsorption of molecules rather than
polymorphic stability. While dispersion interactions are known to
be of key importance in the structural chemistry of molecular and
low-dimensional crystals, they are often overlooked when considering
bulk solids with strong ionic or covalent bonding. Few papers touch
upon the inclusion of dispersion interactions in already established
functionals as a beneficial tool toward the correct calculation of
the energetic ordering of TiO2 polymorphs.^39,48−50^ There is no follow-up research on other compounds.
With the goal of identifying a low computational overhead method
able to reproduce the polymorphic energy ordering of TiO2, as found in experiment, we present here a study of the relative
stability of TiO2 polymorphs spanning a variety of functionals,
which include the GGA functional PBE^51^ and
the hybrid HF-DFT functionals B3LYP,^52,53^ HSE06,^54^ and PBE0,^55^ along
with all their DFT-D3^40^ equivalents. The
DFT-D3 correction scheme was chosen as the representative of dispersion
correction methods in solids; it was shown to provide better accuracy
than DFT-D2 due to the latter’s overestimation of dispersion
interactions and arguably the absence of a three-body term.^6^ The study was extended to MnO2 and
ZnO, both with a richness of crystalline polymorphs, to probe for
wider applicability.
All first-principles calculations were performed using the CRYSTAL17^56^ code, in which crystalline orbitals are expanded as a linear combination of atom-centered Gaussian basis sets. For each material investigated, calculations were performed using PBE, B3LYP, HSE06, and PBE0, with and without Grimme’s D3 correction. The basis sets chosen to represent the different atomic species were all of triple valence plus polarization quality selected from the CRYSTAL Basis Set database and consistently used throughout. The basis sets used are O-8-411d1,^57^ Ti-86-411(d31),^57,58^ Mn-86-411d41G,^59^ and Zn-86-411d31G.^60^
All structures studied have been fully optimized in the space group (SG) indicated by experimental studies. CRYSTAL17^56^ default tolerances have been used for the selection of integrals, SCF convergence and geometry optimization. Reciprocal space integration has been performed via a Monkhorst–Pack mesh^61^ using an 8 × 8 x 8 grid of points for all phases.
Our calculations focus on internal energies only and
do not include
contributions from zero-point energy (ZPE) and vibrational entropy.
All previous studies that estimated the effect of ZPE and vibrational
entropy on the polymorphic phase ordering of TiO2 show
that the impact of these effects are negligible and are unable to
modify the predicted relative stability between rutile and anatase.^9,12,24,30,33−35,50,62^ ZPE and vibrational entropy effects
were thus also neglected for MnO2 and ZnO.
For the
comparison of ferromagnetic and antiferromagnetic orders
in MnO2 polymorphs, calculations employed a broken-symmetry
approach and were initiated from an ionic solution, where each Mn^4+^ ion was assigned a spin-up (α) or spin-down (β)
state as required by the magnetic order studied.
We have investigated nine polymorphs of
TiO2 whose unit cells are illustrated in Figure 1. Table 1 displays a summary of information on the
Ti local environment, space group, and references for the initial
structures used in our geometry optimizations. In most phases, Ti
is 6-coordinated and in an octahedral (O~h) environment. TiO6~ octahedra connect in the lattices
by a varied extent of corner and edge sharing, giving rise to tetragonal,
orthorhombic, and monoclinic unit cells. The exceptions to the 6-coordination
of Ti are baddeleyite and cotunnite, which have Ti in 7- and 9-fold
coordinated environments, respectively. Rutile is known to be the
ground state polymorph of TiO2 from experiment.^20^ The other considered polymorphs are formed at
different experimental conditions according to the TiO2 phase diagram.^33^
Figure 1 Unit cells of the different studied TiO
2polymorphs
In this section, we focus exclusively
on energetics. Structural parameters such as lattice parameters and
bond distances are reported in the SI and
only discussed in a global comparison of results. A collection of
relative energies for all 9 TiO2 polymorphs from experimental
and computational literature, as well as values calculated in the
present work using all standard and D3-corrected functionals is provided
in Table 2. All energies
refer to that of rutile, the experimental ground state of TiO2 at ambient pressure.^20^ Different
functionals provide energy estimates that differ widely not only in
magnitude but also in sign, a situation that is not supportive of
predictive and transferable applications. Reliable thermochemical
measurements are available for only a subset of the TiO2 polymorphs; these are rutile, anatase and brookite^20^ i.e. the three polymorphs observed at ambient pressure.
There are several calorimetry experiments providing relative stability
data between these polymorphs (Table 2). Here, we consider the latest and most accurate measurements
provided by Ranade et al.^20^ Of all the
previous computational literature data considered, only the PBE+U study of ref (35) identifies rutile as the ground state, although it misrepresents
the relative energy of brookite. RPA results instead find rutile to
be stable over anatase, but with an energy difference of only +3.2
meV/f.u. compared to the +27.1 meV/f.u. from experiment. Results of
our calculations indicate that the inclusion of dispersion through
the DFT-D3 post-SCF method substantially modifies the relative stability.
Indeed, while for all uncorrected functionals anatase and brookite
are stable over rutile, the inverse is true upon inclusion of dispersion.
In Figure 2, we examine
in diagrammatic form the relative stability of the rutile, anatase,
brookite and TiO2-B phases, grouping calculated results
into standard and D3-corrected functionals. All the “standard”
functionals incorrectly predict rutile as the least stable of these
four polymorphs. On the other hand, all hybrid HF-DFT D3-corrected
functionals (B3LYP-D3, HSE06-D3, and PBE0-D3) correctly reproduce
the experimental order of phase stability.
Figure 2 Stability of the anatase, brookite, and TiO
2-B polymorphic phases of TiO2relative to the experimental ground state rutile^20^ (ΔE in meV/f.u. TiO2) calculated with the different “standard” and dispersion-corrected functionals. Bar numbers 1–4 refer to “standard” functionals (1: B3LYP, HSE06, PBE, 4: PBE0), 5–8 refer to their dispersion-corrected counterparts (5: B3LYP-D3, HSE06-D3, PBE-D3, PBE0-D3), and 9 refers to experiment where available.^20^
From a quantitative point of view, the HSE06-D3
functional provides
the most accurate estimate of relative energy compared to experiment.^20^ Indeed, the HSE06-D3 errors are of only +2 and
−6 meV/f.u. for brookite and anatase, respectively. Before
corrections, the corresponding HSE06 errors are −28 and −85
meV/f.u.. The above comparison indicates that dispersion contributes
as much as +87 and +42 meV/f.u to the energy of anatase and brookite
relative to rutile, i.e., its quantitative influence is of the same
order of magnitude as the energy difference between polymorphs and,
as such, shows that dispersion interactions cannot be neglected in
the prediction of the relative stability of TiO2 polymorphs.
It is encouraging to observe that DFT-D3-corrected functionals predict
the correct energy ordering for, at least, the rutile, anatase, brookite,
and TiO2-B polymorphs. Similar quantitative effects of
dispersion are observed for the HSE06, PBE, and PBE0 functionals,
while for B3LYP the effect of dispersion on relative energies is even
stronger. While inclusion of dispersion resolves the energetic order
of the phases considered, it does still yield incorrect results for
columbite. This polymorph is more stable than rutile before accounting
for dispersion and inclusion of dispersion enhances the energy difference
due to the higher density of columbite relative to rutile (Table 1).
In Table 2, we also include the relative stability provided by single-point energy calculations using the D3 functionals on geometries previously optimized with the equivalent standard functionals. These data are indicated with the acronym SEP-D3. The results in Table 2 indicate that the change in relative energy arises mostly from the D3 dispersion correction, while structural changes resulting from the inclusion of dispersion play only a minor role.
Many more ceramics are known to have
a similar richness of polymorphic phases as TiO2. Here,
we extend our study to MnO2 and ZnO compositions to verify
whether the trends observed in TiO2 have more general validity.
MnO2 shows polymorphism from various packings of MnO6 octahedra and is of interest for a variety of energy and
environmental applications.^75^ There has
been some amount of work on the phase ordering of MnO2 polymorphs,
although most papers focus on a single form and all use “standard”
DFT methods.^76−81^ Estimating the effect of dispersion in MnO2 is complicated
by the additional contribution of the magnetic order, a recognized
shortcoming of local DFT functionals due to the SIE, which is usually
corrected through DFT+U. However, HF-DFT functionals
also fail to accurately reproduce the energetic ordering of MnO2 polymorphs from experiment, despite resolving the SIE inherent
to GGA DFT. Previous studies have attributed this failure to several
interdependent factors such as the artificial underbinding of O^2–^ ligands and/or an inadequate description
of exchange and correlation in the MnO2 polymorphs.^81^
Contrarily to TiO2, no research
has applied methods beyond DFT to resolving the energetic ordering
of MnO2 polymorphs, due to dispersion being considered
of minor importance compared to the clarification of magnetic order.
However, the SCAN functional is shown to be the method of choice to
yield the correct ordering of polymorphic stability, as it is the
only one to predict β-MnO2, the experimental ground
state,^82^ as more stable than ramsdellite
(R) MnO2.^13,83^ The only calorimetric measurements
available are for the β-/R-MnO2 comparison.^82^ PBEsol, unlike other GGA functionals, shows
the correct order of phase stability.^13,83^ References (13) and (83) also show that the inclusion
of the Hubbard U term can cause errors in hybridization
between O 2p and Mn 3d valence orbitals,
resulting in unfavorable distortions of the MnO6 octahedra.
We study six polymorphs of MnO2 shown in Figure 3: rutile-structured pyrolusite (β-MnO2), the hollandite
α form, the intergrowth γ form, the orthorhombic ramsdellite
(R) form, the spinel λ form, and the layered δ form. Mn
is 6-coordinated in all polymorphs. Table 3 displays information relating to the local
environment of Mn, atomic density, space group, and references for
the initial structures used in our geometry optimizations. Full structural
information for the optimized phases can be found in the SI.
Figure 3 Unit cells of the different studied MnO
2polymorphs.
All calculations concerning MnO2 polymorphs
were performed
with the ground state antiferromagnetic (AFM) order, known from experiment.^13,76,84−86^ The energetic
importance of magnetic coupling in MnO2 has been explicitly
discussed in refs (13), (77), and (85−87). The energetics associated with magnetic order in
MnO2 are greater or of the same order of magnitude as those
involved with dispersion interactions and previous DFT studies have
shown that choosing a ferromagnetic (FM) order yields incorrect results
for lattice energies that the inclusion of dispersion corrections
would be unable to palliate (ΔE^FM/AFM^ ∼ 150 meV/f.u. within the same polymorph using GGA functionals
such as PBE).^78,81^ This was confirmed by our calculations
as exemplified by Table 4.
Following the discussion on TiO2, we evaluate here the effect of dispersion on the energies
of the MnO2 polymorphs. A summary of literature data and
values from our study is given in Table 4. Pyrolusite (β-MnO2) is
shown from experiment to be the ground-state polymorph^82^ and will thus be our reference for the relative
stability.
All “standard” functionals find α-MnO2 to be more stable than the experimental ground state. However,
similar to the TiO2 results, the incorrect ordering is
reversed upon inclusion of dispersion interactions. Although “standard”
PBE, HSE06, and B3LYP functionals correctly predict β-MnO2 to be stable over δ-MnO2, their D3-corrected
counterparts yield a significantly larger energy difference. Both
sets of D3-corrected results for the α- and δ-MnO2 polymorphs match experimental observations.^82^ Results are incorrect when considering R-MnO2, with the exception of B3LYP-D3. R-MnO2 is calculated
as stable relative to β-MnO2, even after inclusion
of dispersion. The energy difference, however, is substantially improved
by including dispersion.
Similar to TiO2, in MnO2 the energy contribution
of dispersion interactions is also of the same order of magnitude
as the energetic difference between polymorphic phases and contributes
to a much improved estimate of relative energies. The case of R-MnO2 may be influenced by the more complex crystal structure that
contains two nonequivalent O sites, pyramidal and planar O, both with
a CN of 3 but differentiated by their different bonding angles. The
more distorted topological connectivity of Opyr, and thus
its different relation to Mn cations compared to Oplan,
may provide an enhanced stabilization of electrostatic nature paired
with the structural complexity leading to magnetic coupling interactions
that are misrepresented by DFT, even when using hybrid exchange functionals.
Similar arguments have been proposed in refs (13), (78), and (93). The results for columbite
TiO2 and R-MnO2 thus highlight that the inclusion
of dispersion interactions needs to be paired with an already accurate
estimation of polymorphic stability to have an impact on results,
especially when energetic contributions of the same magnitude or larger,
such as magnetic coupling or change in coordination environment of
some of the ions, may affect relative stabilities.
The third system we consider is
ZnO, of interest among other applications as a photocatalyst.^93^ ZnO crystallizes as hexagonal wurtzite (HW)
and cubic zinc-blende (ZB) at standard conditions, with the former
being the ground state.^94^ We also include,
for comparison, two high pressure phases, rocksalt-type (RS) and CsCl-type
ZnO shown in Figure 4. Table 5 displays
structural information relating to the ZnO polymorphs studied here.
Full structural information on the optimized structures is provided
in the SI. Studies on the relative stability
of ZnO polymorphs include DMC^95^ and RPA^96^ calculations as well as “standard”
DFT-GGA techniques.^97−99^ All calculations find wurtzite as the ground state,
however the relative energy of rocksalt-type ZnO is largely overestimated,
and there is as yet no estimate of dispersion contribution on polymorph
stability. This is an interesting topic to examine as, unlike MnO2 polymorphs where Mn is always in an octahedral environment,
in ZnO the CN of Zn changes from 4 to 6 and 8 in the phases studied.
Figure 4 Unit cells of the different ZnO polymorphs studied.
A summary of relative energies from
literature and the present study is presented in Table 6. Similar to TiO2, there are several calorimetry experiments evaluating the relative
stability of the HW and RS polymorphs. The reference value used for
comparison to experiment here is given by Sharikov et al.,^104^ as previous studies^105^ erroneously neglect the kinetic features of the HW to RS phase transition
below 1000 K.
All functionals find HW stable over ZB and the effect of dispersion is negligible in the energy difference between these two phases, which are in practice just different stacking motifs of the same local Zn environment (also highlighted by the polymorphs’ close to identical atomic densities). The effect of dispersion is instead pronounced (over −100 meV/f.u. for B3LYP, PBE, and PBE0) when considering the relative stability of the RS polymorph, that is denser than wurtzite thanks to its higher CN. Uncorrected DFT/HF-DFT results largely overestimate the calculated calorimetric value, while after the inclusion of dispersion, the HW-RS energy difference is significantly closer to experiment, despite showing a large dependence on the functional. Once more, dispersion has similar magnitude as the energy difference between polymorphs (both of the order of 100 meV/f.u.) and must thus be accounted for in calculated thermochemistry.
A solid structure’s local
atomic environment will have the strongest impact on the absolute
value of dispersion forces, notably due to D3 dispersion coefficients
being determined by an ion’s coordination number and the largest
contributions occurring at small interatomic distances r with dispersion forces scaling as 1/r^6^. The distribution of bond distances in the crystal lattice can be
monitored through the pair radial distribution function g(r) (RDF). The RDFs for two pairs of phases whose
relative energy is affected by dispersion, i.e., rutile/anatase TiO2 and HW/RS ZnO, are compared in Figure 5. As we can see from Figure 5, there are noticeable changes in the local
environment between the polymorph pairs, including in the short-range
environments (between 2 and 4 Å). Rutile and anatase phases of
TiO2 both have Ti in 6-fold coordination and O~h~ environment. The Ti–O bond distances show
negligible changes in the RDF; however the phases are differentiated
by next neighbor relations. The O–O distances within the same
octahedron over the 2–3 Å range differ between rutile
and anatase, due to the more distorted nature of the latter. This
is exemplified by the three peaks at 2.4, 2.8, and 3.0 Å for
anatase compared to the single, high amplitude peak at 2.75 Å
for rutile. The Ti–Ti next nearest neighbor distances visible
in the peaks between 3 and 4 Å suggest shorter next nearest neighbor
bond distances overall in rutile.
Figure 5 RDFs of rutile and anatase TiO
2(top) and hexagonal wurtzite and rocksalt-type ZnO (bottom) computed with the HSE06 functional. The red lines represent rutile TiO2and HW ZnO, respectively. The blue lines represent the anatase TiO2and RS ZnO, respectively.
As the stabilizing contribution to the energy given by dispersion scales as 1/r^6^, rutile’s more contracted local environment explains why this phase is stabilized over anatase when such dispersion forces are accounted for. In the two considered ZnO phases, instead, the CN changes from 4 to 6 between the HW and RS phases. Differences in the RDF are already obvious for the nearest Zn–O distances, but also have appreciable contributions from further neighbor shells. The higher overall density of the RS phase explains its stabilization by dispersion forces relative to wurtzite.
While the RDF contains the atomic-level information necessary to
rationalize the effect of dispersion on polymorph stability, discriminating
the contribution from individual atomic pairs is challenging. It would
be more useful to be able to employ a global property of the materials
to rationalize results. Atomic density, expressed in number of atoms
per unit volume, is a simple but appropriate measure of interatomic
separations in crystal lattices. We have therefore investigated whether
correlations exist between dispersion energy and atomic density. We
first define the dispersion contribution to polymorph stability, ΔED3, as
where ΔEstd represents the relative stability between two polymorphs calculated
with “standard” functionals and ΔEcorrec represents the relative stability between two polymorphs
calculated with D3-corrected functionals. ΔED3 represents the sum of two contributions here i.e. the
influence of the D3 correction on the single-point energy of the crystal
lattice and the change in geometry stemming from the inclusion of
dispersion forces through the D3 correction.
Table 7 provides
a comparison of the calculated ΔED3 for anatase relative to rutile TiO2 and RS relative to
HW ZnO using the HSE06 functional as an example. Despite differing
contributions from dispersion based on the unit cell composition,
we observe a correlation between atomic density and contribution from
dispersion interactions. Indeed, with a less dense cell than rutile,
ΔED3 destabilizes anatase relative
to rutile, while the inverse is observed for the HW and RS ZnO polymorphs.
Figure 6 shows plots
of ΔED3 against the atomic densities
for the TiO2 polymorphs investigated in this work, using
the PBE and B3LYP functionals. Plots for other compounds (MnO2/ZnO) along with a complete plot for TiO2 can be
found in the SI. In both cases there is
a strong near-linear correlation, showing that the dispersion contribution
increases as a function of density. The difference can be of several
hundreds meV/f.u. due to a large change in density originating from
a CN change as we can see for CsCl-type ZnO, for example. Smaller
differences are observed between phases with the same CN and similar
density, but these remain non-negligible and should therefore be included
to correctly differentiate the small energy differences between these
phases. We have performed a linear fit of the calculated ΔED3 values as a function of the atomic density
for each combination of composition and functional. The gradients
ΔED3^avg^ are reported in Table 8. The slope of the best fit curve is evidently
compound and functional dependent, with the strongest dispersion contributions
occurring in TiO2 and B3LYP results. PBE, PBE0, and HSE06
show comparable values.
Figure 6 Energetic contribution from the DFT-D3 correction, ΔE
D3in meV/f.u. TiO2is plotted against the atomic density in atoms per cubic Ångström (a.Å^–3^) of each TiO2polymorph for the representative functionals, PBE and B3LYP. The data coloring is as red for PBE and blue for B3LYP.
The near-linear dependence on atomic density allows us to estimate the D3 correction energy as
where ΔED3^avg^ can be obtained
from the quasi-linear relationship between change in atomic density
and calculated lattice energies by D3-corrected functionals (gradients
of the linear fit in Figure 6 and Table 8). From its definition and the definition of its constituting parts,
the calculation of ΔED3 is evidently
system-dependent and each calculated ΔED3 value for a specific oxide will only be applicable to polymorphs
of that system, as evidenced by the results in Table 8. Despite being within the same order of
magnitude, we can observe some variations in calculated values of
ΔED3^avg^ between TiO2, MnO2 and ZnO.
Our results show that dispersion interactions are relevant to the relative stability of crystalline oxides, due to the energy contribution provided being of the same order of magnitude as the lattice energy difference between polymorphs. We have evidenced a near-linear dependence between atomic density and the energy correction provided by dispersion interactions (Figure 6 and eq 4). This is not the denser a material, the shorter its interatomic distances, r, resulting in a larger contribution from attractive dispersion interactions. The differences are far from inconsequential, often amounting to tens of meV/f.u..
A full list of structural parameters for the geometry optimized structures of all compositions, phases and functionals employed in the current work is provided in the SI. The number of results provided is too great to examine individually; we therefore employ a global analysis to examine the effect of dispersion on structural information. For each composition, phase and functional examined, we consider the error with respect to experiment of the three lattice parameters, equilibrium volume and shortest M–O bond distances (Δa, Δb, Δc, δV, δr), respectively. The errors in each lattice parameter are then used to calculate an average error, δd = . Using the equilibrium volume as an example, we define the mean and absolute errors as δV = and |δV| = , respectively.
Table 9 presents the calculated errors for each “standard” and D3-corrected functional. Results presented in Table 9 are limited to phases for which an experimental structural determination at ambient pressure conditions is available and correspond to static-lattice calculations. The inclusion of the D3 contribution always causes a contraction of the equilibrium volume; for the functionals overestimating V (PBE), adding D3 improves the equilibrium geometry. If instead V is close to the experimental value or underestimated at the DFT level (HF-DFT), D3 causes excessive contraction and the equilibrium volume deviates slightly from experiment. The effect on nearest neighbor distances is negligible (δr in Table 9) and only nonbonded distances are affected by the inclusion of D3. Overall, the effect of D3 on geometries is small, and the improvement yielded in relative stabilities largely outweighs any loss of accuracy on equilibrium geometry.
To conclude, this work highlights the importance of the inclusion of dispersion interactions for the correct computational estimate of energetic phase orderings. This is due to both the energy differences between polymorphic phases and the dispersion interactions’ energy contribution being of tens of meV/f.u.. Due to the sparse availability of measured calorimetric data, future investigations must not only involve computational studies on an even wider range of materials to observe how distinct the influence of dispersion is in different compositions, but also new experimental work to provide a more complete and reliable experimental data set. It would also be of interest to evaluate the results obtained by DFT-D4,^109,110^ the successor to DFT-D3 used in this paper or functionals self-consistently accounting for dispersion.
Via the authors’ membership of the UK’s HEC Materials Chemistry Consortium, which is funded by EPSRC (EP/R029431), this work used the Archer2 UK National Supercomputing Service (http://www.archer2.ac.uk). We are grateful to the UK Materials and Molecular Modelling Hub for computational resources, which is partially funded by EPSRC (EP/T022213/1). The authors acknowledge the use of the UCL Katheen High Performance Computing Facility (Kathleen @UCL), and associated support services, in the completion of this work.
The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jpcc.3c01013.
The authors declare no competing financial interest.