Tin Monoxide: Structural Prediction from First Principles Calculations

ARTICLE pubs.acs.org/JPCC

Tin Monoxide: Structural Prediction from First Principles Calculations with van der Waals Corrections Jeremy P. Allen,*,† David O. Scanlon,† Stephen C. Parker,*,‡ and Graeme W. Watson*,† † ‡

School of Chemistry and CRANN, Trinity College Dublin, Dublin 2, Ireland Department of Chemistry, University of Bath, Claverton Down, Bath, BA2 7AY United Kingdom ABSTRACT: Tin monoxide is a technologically important ptype material which has a layered structure dictated by nonbonded dispersion forces. As standard density functional theory (DFT) approaches are unable to account for dispersion forces properly, they routinely give rise to a poor description of the unit cell structure. This study therefore applies two forms of empirical dispersion corrections, using either atomic- or ionic-based parameters for the dispersion coefficients, to assess their ability to correctly model the atomic structure and the formation energies of the important p-type defects. Although both approaches show an improvement in the predicted unit cell structure over that with no dispersion corrections, the ionic-based parameter set shows significantly better results, with lattice vectors reproduced within 0.2% of experiment. The atomic-based parameters still predict a distorted cell though, which is carried through to the defective system. On the introduction of defects, a similar degree of structural relaxation is observed regardless of the approach. The defect formation energies, however, are seen to differ more substantially, with the atomic-based set giving an overestimation of the energies due to excessive Sn
Sn interactions. Overall, this study shows that empirical van der Waals corrections utilizing an ionic-based parameter set can be used to model SnO.

’ INTRODUCTION The tin
oxygen system has two important oxide phases, tin(II) monoxide, SnO, and tin(IV) dioxide, SnO2. The dioxide is thermodynamically more stable at standard conditions and it has therefore been more widely studied,1
3 having important applications in gas sensors4
6 and transparent conducting layers.7,8 In contrast to the n-type properties of SnO2, SnO shows good p-type performance9
11 and has received attention due to its potential use as a p-type thinfilm transistor.9,10,12
14 SnO has also recently been reported to be a bipolar material, due to the formation of n-type SnO through Sbdoping, indicating that its conduction band minimum is relatively low on an absolute scale.14 Nanoparticulate SnO has also been investigated for use as an anode material in rechargeable batteries.15 SnO can be formed via the decomposition of SnO2 above 1773 K. However, it also exhibits poor stability, in comparison to the dioxide form, with SnO oxidizing to SnO2 at around 573 K in air.16 SnO is known to crystallize in two polymorphic forms, a litharge (α) and a massicot (β) form, with the former being the most stable. The tetragonal litharge structure (space group P4/nmm) possesses 4-coordinate Sn and O atoms, with the latter showing a regular tetrahedral environment. The bonding around the Sn atoms is distorted, with the four oxygen atoms forming the base of a square pyramid and the top vertex being occupied by a lone pair. This asymmetric bonding environment gives rise to a layered structure, with ab sheets which stack in the c direction, as shown in Figure 1. The atomic structure and bonding in SnO has been the subject of a number of computational studies using density functional theory (DFT), to rationalize both the nature of the lone pairs present on the Sn(II) atoms and the stability of its layered nature over a symmetric, regular structure.18
22 The asymmetric structure r 2011 American Chemical Society

was shown to result from the lone pair density on the Sn(II) species. However, this lone pair was not inert, as classical theory predicted, but was a direct result of the bonding. The 2p states of the oxygen anion were seen to mix with the 5s states of the Sn, creating a filled antibonding orbital, rather than a nonbonding orbital. This could mix with the Sn 5p states to stabilize the antibonding states but only in asymmetric environments causing the distorted structure.23 The importance of the bonding on the formation of the lone pair is highlighted through consideration of the other Sn chalcogenides. Changing the anion from S through to Te results in reduced mixing with the Sn 5s states, due to the increasing energy of the anion p states. For Te, the 5p states are too high in energy to stabilize the antibonding orbital, and therefore the lone pair, so a symmetric structure results.21,24 Other materials where the cation oxidation state is 2 less than the group valence have also been seen to have similar methods of lone pair formation, including PbO25
29 and Bi2O3.30 Experimentally, a number of optical band gaps have been reported for SnO, ranging from 2.5 to 3.3 eV,31,32 although more recent studies tend to bring this closer to 2.7 eV.9,12 In addition, an experimental indirect band gap of 0.7 eV has been reported from diffuse reflectance spectra.9 Unsurprisingly, DFT simulations tend to predict an indirect fundamental gap of a little less than this, ranging from 0.333 to 0.6 eV,18 depending on the functional and basis set employed in the calculation. The smallest direct fundamental band gap is typically underestimated to a greater magnitude, with values reported of ∼2.0 eV.34,22 Received: June 1, 2011 Revised: August 24, 2011 Published: September 02, 2011 19916

dx.doi.org/10.1021/jp205148y | J. Phys. Chem. C 2011, 115, 19916–19924

The Journal of Physical Chemistry C

ARTICLE

where Rij is the distance between the atoms i and j, Cij6 is the dispersion coefficient between the same atoms and s6 is a scaling factor reflecting the fact that DFT functionals may include a portion of long-range attraction. fdmp is a damping function employed to screen interactions at very small distances and to avoid near-singularities DFT
D2 ðRij Þ ¼ fdmp

1 1 þ

ij

e
dðRij =Rr
1Þ

ð3Þ

where d is the damping parameter, which is related to the steepness of the damping function, and Rijr is the sum of the van der Waals radii of atoms i and j, Rir and Rjr, respectively. The dispersion coefficient between two atoms i and j, Cij6, is calculated from ij

C6 ¼

Figure 1. Atomic structure of the unit cell of litharge-structured SnO.17 Tin and oxygen atoms are colored blue and red, respectively.

Despite previous studies providing a reasonable description of the bonding in SnO, discrepancies are seen in the structural parameters.18,20,21,34 This typically occurs as an expansion of the c lattice vector relative to the a, which is also the direction of the nonbonded van der Waals, or dispersion, interactions between the SnO layers. The failure of electronic structure methods commonly applied to solids (for example, DFT, Hartree
Fock and hybridDFT) to account for van der Waals interactions is well-known35
37 and has been reported for a number of other layered materials, such as V2O538,39 and MoO3.40,41 In addition, both Christensen et al.42 and Errico33 have shown that the calculated band gap is dependent on the c/a ratio. Although standard DFT methods would not be expected to correctly predict the fundamental band gap, due to their discontinuity problem,43 this implies that the failure to model the structure will also introduce an error into the calculated band gap. A fundamental band gap, as used in other studies,44
46 is defined as the smallest possible gap of an indirect/ direct nature between the valence band and the conduction band. A number of different methods exist which work to apply dispersion corrections to DFT methodologies, as discussed by Dobson et al.36 These include, but are not limited to, the nonlocal van der Waals functional (vdW-DF) of Dion et al.,47 which modifies the correlation functional, and the DFT-D methods of Grimme and co-workers,48
50 which apply asymptotic dispersion corrections in a pairwise manner. Of the latter, the DFT-D2 method49 in particular has attracted a great deal of attention for a wide range of systems.51
54 This method acts by correcting the self-consistent Kohn
Sham energy, EKS
DFT, with an empirical dispersion correction, Edisp. The total corrected energy, EDFT
D2, is therefore given as EDFT
D2 ¼ EKS
DFT þ Edisp

ð1Þ

One of the advantages of this approach is that it can be used with a range of different density functionals due to it being a simple additive term. Furthermore, its simplicity means that is computationally inexpensive for large or complex systems. Its parametrization also makes it flexible and applicable to a range of different systems via parameters based on empirical or semiempirical physical quantities. The dispersion correction is determined via ij

Edisp ¼
s6

C6 fdmp ðRij Þ Rij 6

ð2Þ

qffiffiffiffiffiffiffiffiffiffi j Ci6 C6

ð4Þ

where the Ci6 terms are the dispersion coefficients between two atoms of the same type. The atomic Ci6 values were determined for each species using an approach derived from the London formula for dispersion, based on DFT/PBE0-calculated ionization potentials, Ip, and static dipole polarizabilities, α.55 It is equal to Ci6 ¼ 0:05NIpi αi

ð5Þ

where N has values of 2, 10, 18, 36, and 54 for atoms in rows 1
5 of the periodic table. The value of the s6 term is dependent on the functional used and the Ci6 terms. Grimme49 fitted this for a range of functionals, based on his calculated Ci6 values, using a least-squares optimization of the deviations from experiment for the calculated binding energies of 40 noncovalently bound complexes. Although the DFT-D methods of Grimme have not been directly tested on SnO before, Duan52 modeled the bulk structure and low-index surfaces of SnO using the approach of Wu and Yang.37 This method is similar to the DFT-D methods, in that the dispersion corrections are applied in an identical manner via eqs 2 and 3, but the Cij6 values are determined using ij C6

2 3 i j i j 3 4 α α Ip Ip 5 ¼ 2 Ipi þ Ipj

ð6Þ

Duan calculated these parameters using tabulated dipole polarizabilities and ionization potentials,56 finding similar Cij6 values to those of Grimme. The van der Waals radii used by Duan for Sn and O were also larger than those given by Grimme. However, the reported lattice vectors for SnO only showed limited improvement in comparison to those determined without any dispersion corrections, with the difference to experiment reducing from 1.28% to 0.19% for the a vector and from 2.70% to 1.59% for the c vector. Therefore, Duan's approach failed to fully account for the dispersion interactions existing between the layers. The purpose of this study is to study the effects of van der Waals corrections on the bulk structure of SnO, using both an alternative form of the damping function and a different set of parameters to define the dispersion corrections based on calculated ionic polarizabilities. Furthermore, the effect of the different methods on the structure and formation energies of the important p-type defects, namely a tin vacancy and an oxygen interstitial, is considered. 19917

dx.doi.org/10.1021/jp205148y |J. Phys. Chem. C 2011, 115, 19916–19924

The Journal of Physical Chemistry C

ARTICLE

’ COMPUTATIONAL METHODS All calculations contained in this study were performed using the periodic DFT code VASP.57,58 This uses a plane wave basis set to describe the valence electronic states, with interactions between the cores (Sn:[Kr] and O:[He]) and the valence electrons described using the projector-augmented wave59,60 (PAW) method. The exchange and correlation is treated with the generalized gradient approximation (GGA) via the Perdew
Burke
Ernzerhof 61 (PBE) functional. Equilibrium lattice parameters for the bulk systems were determined through a series of structural minimizations at different volumes. In each calculation, the atomic coordinates, lattice vectors and angles were allowed to fully relax within a fixed total unit cell volume. The resultant energy-volume curves were then fitted to the Murnaghan equation of state to obtain the bulk equilibrium cell volume.62 This approach avoids potential problems of Pulay stress and changes in the basis set which can accompany volume changes in plane wave calculations. The Pulay stress can be anisotropic and, to confirm the accuracy of this approach, a single-point calculation of a minimized structure at a cutoff of 1000 eV was performed, finding the pressure along the different vector directions to vary by less than 0.2 kbar. A plane wave cutoff of 500 eV was used for all calculations and Monkhorst
Pack63 k-point grids of 6  6  3 and 2  2  1 were used for the bulk and 3  3  3 supercell calculations, respectively. For each calculation, structural minimization was deemed complete when the forces on all atoms were less than 0.001 eV Å
1. For modeling the bulk structure of α-Sn metal, a k-point grid of 11  11  11 was used, whereas for gas phase O2, a large 16  16  16 Å cell was used with a 1  1  1 k-point grid. The O2 molecule was calculated using spin polarization, with the species being in the triplet state. For the defect calculations, structural minimization was achieved under a constant volume regime, with cell vectors and angles constrained to those determined for the pure bulk cell. As the cell dimensions were constrained for defect calculations, the convergence criteria for the forces on the atoms was set to 0.01 eV Å
1. In addition to the standard DFT-D2 method with the PBE functional (PBE-D2), we also used an alternative method, which will hereafter be referred to as PBE-vdW. This acts in the same way as the PBE-D2 approach but with a different damping function, similar to that used by Gonzalez and Lim,64 given by PBE
vdW fdmp ðRij Þ ¼

1 ij

1 þ e
4BðRij
Rr Þ

Table 1. Parameters Used in the PBE-D2 and PBE-vdW Methods for Modeling SnO PBE-D2 Rir (Å)

6

ion

C6 (eV Å )

Sn

401.20

1.804

20

30

O

7.26

1.342

20

30

d

cutoff (Å)

PBE-vdW ion pair

C6 (eV Å )

Rijr (Å)

B (Å
1)

cutoff (Å)

Sn
Sn

73.07

3.693

0.5

20

Sn
O

74.45

2.222

0.5

20

O
O

77.34

2.688

0.5

20

ij

6

ionic polarizability was determined using Hartree
Fock calculations. The interatomic distances (Rijr ) for the Sn
O and O
O interactions were taken from the smallest in-layer distances. However, for the Sn
Sn interaction, the across-layer distance was used. This is because the interlayer distance is controlled by the Sn
Sn van der Waals interactions and therefore it is more appropriate to use this distance to define the value of Edisp than the in-layer distance. The value of the damping parameter, B, was determined by plotting the damping function and selecting the maximum value which did not result in the force becoming negative, as will be discussed later, hence ensuring that the damping function does not affect interactions beyond nearestneighbor. The s6 scaling factor used by Grimme is 0.75, which they recommend for all calculations using the PBE functional with the DFT-D2 method. For the PBE-vdW approach, we reduced this to 0.65 through fitting to the experimental structural data. This fitting was achieved by adjusting the s6 value until the errors in the lattice vectors, c/a ratio and bond lengths compared to experiment were minimized. As these quantities are all intimately connected no complex fitting procedure was required. The defect formation energies, Ef(D), are calculated using the following equation: 1 Ef ðDÞ ¼ Etot ðDÞ
EðSnOÞ þ nO EðO2 Þ 2 þ nSn Eðα
SnÞ þ nO μO þ nSn μSn

ð7Þ

This variation uses the difference between Rij and Rijr , rather than the ratio, as used in the PBE-D2 approach. In this way, the distance over which the damping function operates is independent of the reference distance, which is not the case for the fractional distance relationship of the PBE-D2 method. The damping functions used in both the PBE-vdW and PBE-D2 method will reduce the dispersion interactions to 50% of their undamped value when the interatomic separation is equal to Rijr . However, rather than determining Rijr directly from the sum of van der Waals radii, as Grimme does, we use interatomic distances from the experimental structure.17 This allows the value of Edisp to be more specific to the system in question, where the Sn atoms have an asymmetric charge density. All the parameters for the PBE-D2 and PBE-vdW methods used in this study are detailed in Table 1, with the parameter set for the PBE-D2 method taken directly from the work of Grimme.49 For the PBE-vdW approach, the Cij6 parameters were taken from the work of Madden and co-workers,65,66 where the

i

ð8Þ

where Etot(D) is the total energy of the supercell containing the defect D and E(SnO), E(O2), and E(α-Sn) are the total energies of the pure SnO supercell, an O2 molecule and a Sn atom in bulk α-Sn, respectively, where dispersion corrections are applied to all systems. The number of oxygen and tin atoms either added to or taken away from an external reservoir are given by nO and nSn, respectively, whereas μO and μSn are the chemical potentials of oxygen and tin, respectively. The defect formation energies are calculated only at the Sn metal/SnO2 boundary due to SnO having no thermodynamic stability field, as has been suggested in other studies.67
69 This will discussed later with the results of the defect calculations.

’ RESULTS Bulk Structure. The calculated lattice vectors and bond lengths using the PBE, PBE-vdW and PBE-D2 methods are given in Table 2, along with experimental values.17 Unsurprisingly, PBE fails to model the SnO structure, as seen in previous studies.18,20,34,52 19918

dx.doi.org/10.1021/jp205148y |J. Phys. Chem. C 2011, 115, 19916–19924

The Journal of Physical Chemistry C

ARTICLE

The error in the c lattice vector, relative to experiment, is seen to be 2.7% larger than that seen for the a lattice vector. Examination of the in-layer and across-layer bond distances show that this is primarily due to an expansion between the layers, resulting from the failure to account for the nonbonded dispersion forces acting in this direction. Some studies have reported GGA-calculated c/a ratios in close agreement with experiment. However, this is likely to be a result of the convergence criteria used in these studies. Although the majority of studies do not detail the convergence criteria used, Liu et al.22 considered their cell converged when the forces on the atoms were less than 0.01 eV Å
1 and Togo et al.34 used a value of 0.05 eV Å
1. The reported c/a ratios from the work of Liu et al. and Togo et al. were 1.278 and 1.283, respectively. These ratios are both smaller than the current study (1.306) and that of Walsh and Watson20 (1.302) which were both determined using convergence criteria of 0.001 eV Å
1. This suggests that the energy curve for this system may be shallow enough that a smaller c/a can be found if the converge criteria are not rigorous enough. Similar discrepancies in the c/a ratio of the isostructural α-PbO are seen using GGA.28,70 Table 2. Experimental Lattice Vectors and Bond Lengths for SnO with Those Calculated Using PBE, PBE-vdW, and PBED2 Methodsa property

PBE-vdW

PBE-D2

experiment17

a

3.860 (1.6) 3.801 (0.0) 3.835 (0.9)

3.801

c

5.042 (4.3) 4.843 (0.2) 4.820 (
0.3)

4.835

c/a

1.306 (2.7) 1.274 (0.2) 1.257 (
1.2)

volume u Sn
Sn (in-layer)

a

PBE

75.14 (7.6) 69.97 (0.2) 70.90 (1.5)

1.272 69.85

0.2300 0.2384 0.2422 3.582 (1.2) 3.543 (0.1) 3.579 (1.1)

0.2381 3.539

Sn
O (in-layer)

2.252 (1.4) 2.224 (0.1) 2.245 (1.0)

2.222

O
O (in-layer)

2.730 (1.6) 2.687 (0.0) 2.712 (0.9)

2.688

Sn
Sn (across-layer) 3.855 (4.4) 3.694 (0.0) 3.678 (
0.4)

3.693

Sn
O (across-layer) 4.336 (4.6) 4.149 (0.1) 4.125 (
0.5)

4.145

O
O (across-layer)

4.835

5.042 (4.3) 4.843 (0.2) 4.820 (
0.3)

Percentage changes to experiment are given in parentheses and all values are in Å, except cell volume which has units of Å3 and the c/a ratio which is dimensionless. The u value relates to the fractional coordinate of the Sn atom in the unit cell, which has coordinates of (0, 0.5, u).

While the lack of dispersion forces are not the only source of error in GGA-DFT, the introduction of both types of dispersion corrections is seen to have a dramatic improvement over the PBE results. For PBE-vdW, the dispersion corrections are seen to reduce the c vector and across-layer distances more than the a vector and the inlayer distances. These reductions result in the error compared to experiment being equalized, with lattice vectors/bond lengths either being the same as experiment or at least within 0.2% of the value. This also gives a c/a ratio which is in close agreement with experiment. With the PBE-D2 approach, the c vector is also seen to be reduced from the PBE results more than the a vector, as expected. The error in the a vector reduces from 2.7% to 0.9%, whereas the c vector shows a reduction in the error from +4.3% to
0.3%. The in-layer and across-layer bond distances are also seen to reduce by equivalent amounts. Despite this improvement, the resulting unit cell is still distorted in comparison to experiment, as shown by the c/a ratio being 1.2% smaller than experiment. However, the distortion seen is the opposite to that found with PBE, with the a vector being too large and the c vector being too small. This implies that the dispersion corrections used in the PBE-D2 approach are not correcting for the dispersion forces in the correct manner, by not sufficiently accounting for the in-layer interactions but overestimating the across-layer van der Waals forces. It is also instructive to consider the effect of the van der Waals corrections on the atomic positions of Sn atoms in the unit cell, as detailed in Table 2. For PBE, the u value is too small, resulting from the overexpansion of the c vector causing Sn atoms to be further away from each other than seen experimentally. On the addition of dispersion corrections, the PBE-vdW method is again seen to perform significantly better than the PBE-D2 approach, giving to rise to a Sn position in excellent agreement with experiment. For the PBE-D2 approach, the Sn u value is too large, which, when coupled to the contracted c vector, indicates that the Sn
Sn acrosslayer distance is too small, indicating that the Sn
Sn dispersion is excessive. Overall, the structural data suggests PBE-vdW method is considerably more effective than the PBE-D2 approach. Although prediction of the atomic structure is of key importance to this study, the additional effects that it has on the calculated electronic structure also need to be considered. The calculated band structures along the high symmetry points of Bradley and Cracknell71 using the three methods are shown in Figure 2. The fundamental band gaps, both indirect and direct, are seen to be severely underestimated for all three methods in comparison to

Figure 2. Band structure of SnO as calculated using (a) PBE, (b) PBE-vdW, and (c) PBE-D2. The top of the valence band is set to 0.00 eV and indicated by the red (dashed) line. 19919

dx.doi.org/10.1021/jp205148y |J. Phys. Chem. C 2011, 115, 19916–19924

The Journal of Physical Chemistry C

ARTICLE

experiment. In accordance with previous studies, the PBE method gives a direct gap at the Γ point of 1.91 eV and an indirect gap between Γ and M of 0.45 eV. Experimentally, the optical band gap of SnO is reported as being ∼2.7 eV9,12 and the indirect band gap is 0.7 eV.9 The inclusion of dispersion corrections is seen to worsen this though, despite the improved structural reproduction. For the PBE-vdW and PBE-D2 approaches, indirect fundamental gaps between the valence band maximum (VBM) at the Γ point and the conduction band minimum (CBM) at the M point are calculated as 0.03 and 0.02 eV, respectively. Direct fundamental gaps of 1.94 and 1.90 eV are found for the PBE-vdW and PBE-D2 methods, respectively, which are similar to that predicted by PBE. However, despite being quantitatively similar, the location of the smallest direct gap is different, with both dispersion-corrected methods finding the M point to possess the smallest direct gap. The cause of this is seen to be a change in the energy of the CBM at M. Both the PBE-vdW and PBE-D2 methods show that the CBM lowers in energy relative to the rest of the band, increasing the curvature of the band. This gives rise to better agreement in the difference between the direct and indirect fundamental band gaps in the dispersioncorrected methods. As the curvature of the CBM is directly related to n-type properties, with a greater curvature indicative of a smaller band edge electron effective mass, this suggests that the n-type properties may be better than PBE calculations suggest. p-Type Defects. As the p-type defects, namely oxygen interstitials and tin vacancies, are more dominant in SnO and thus more widely studied, the effect of dispersion corrections on these defects has been considered. To determine the defect formation energies, within the stability field of SnO, the chemical potentials of Sn and O, μSn and μO, need to be calculated. These chemical potentials reflect the specific growth conditions within the global constraint of the enthalpy of formation of the parent material, ΔHSnO f SnO μSn þ μO ¼ ΔHf ð9Þ At the boundary with the metal, typically referred to as Snrich/O-poor, the limit is defined by μSn ¼ 0:00eV

ð10Þ

μO ¼ ΔHfSnO

ð11Þ

The O-rich limit is defined by the formation of SnO2 and therefore has the constraint μSn þ 2μO e ΔHfSnO2

ð12Þ

By using eqs 9 and 12, the chemical potentials at the Sn-poor/ O-rich limit can be obtained by μSn ¼ 2ΔHfSnO
ΔHfSnO2 μO ¼

ΔHfSnO2


ΔHfSnO

ð13Þ ð14Þ

The experimental enthalpies of formation for SnO and SnO2 are
2.96 and
6.02 eV, respectively.72 By substituting these into the eqs 13 and 14, chemical potentials of +0.09 and
3.06 eV are found for μSn and μO, respectively, at the Sn-poor/O-rich limit. The positive value of μSn indicates that there is no thermodynamic stability field for SnO, since the O-poor limit has μSn of 0.00 eV, and is consistent with its poor stability at standard operating temperatures.2,16 This means that it is not possible to define rich and poor limits for Sn and O. In this study, we have calculated the defect formation energies at the Sn metal/ . For the SnO2 boundary, where μSn = 0.00 eV and μO = ΔHSnO f PBE, PBE-vdW and PBE-D2 methods, μO is calculated to be
2.65,
3.03, and
3.02 eV, respectively, per formula unit.

Contrary to experiment, it is possible to define a very narrow thermodynamic stability field for SnO with GGA methods. Togo et al.34 and Oba et al.73 defined Sn-rich/O-poor and Sn-poor/Orich chemical limits in their study using the same approach as described by eqs 9 to 14. At the Sn-poor/O-rich limit, they calculated μSn and μO values of
0.23 and
2.49 eV, respectively. Using the PBE results in the current study, similar chemical potentials can be calculated to those of Togo et al. and Oba et al. Other models also give a small stability field for SnO, however, previous studies67
69 have suggested that it is not appropriate to use this erroneous stability field when it is in disagreement with the experimental results. In either case, the variation in μ is small and will therefore have little impact on the results. Togo et al.34 and Oba et al.73 successfully modeled the native defects in both neutral and charged states to help understand the origin of the reported p-type conductivity in SnO. They found that the tin vacancy is the dominant defect and thus is the main contributor to the p-type properties, through the formation of shallower acceptor states. At elevated temperatures there will be similar concentrations of the oxygen interstitial. However, these oxygen interstitials will mainly remain in their neutral state and therefore will be unlikely to give rise to p-type conductivity. As the tin vacancy and oxygen interstitial defects are the most dominant, we elected to use these defects, in their neutral state, to assess the effect of dispersion corrections on defect calculations. For the defect structures, all three methods gave the same basic structural relaxation on the introduction of the defects and are shown in Figure 3. Table 3 details bond lengths and distances, along with atomic displacements, following minimization of the defects. The introduction of the Sn vacancy, indicated by the green atom in Figure 3, panels a and b, results in four oxygen atoms becoming under-coordinated. These four oxygen atoms are then seen to relax from their original positions, colored yellow in Figure 3b, moving away from the vacancy. The value of this displacement is 0.226, 0.232, and 0.253 Å, for PBE-D2, PBEvdW, and PBE, respectively. The Oi, colored green in Figure 3, panels c and d, is positioned in the interlayer gap 2.760, 2.786, and 2.851 Å away from the Sn atoms in the layer above, for PBE-vdW, PBE-D2, and PBE, respectively. The only significant atomic displacement seen in this case is for the Sn atom to which the interstitial oxygen is directly bonded. This is shown in Figure 3d by the yellow atom which represents the position of the Sn in the pure supercell. In all three cases, the Sn atom is seen to move toward the oxygen layer by 0.3 Å, giving a Sn
Oi bond length of 1.9 Å. This in turn affects the in-layer Sn
O bonds from this Sn atom, reducing them by 0.153 Å for PBE/PBE-vdW and by 0.157 Å for PBE-D2. For the Sn vacancy, the calculated defect formation energies are 1.53, 1.42, and 2.28 eV for the PBE, PBE-vdW, and PBE-D2 methods, respectively. The corresponding formation energies for the oxygen interstitial are 1.25, 1.20, and 1.50 eV. The results of the PBE-vdW method therefore give similar, albeit slightly smaller, formation energies to those determined from the PBE approach. However, larger differences are observed for the PBED2 method. Although the O interstitial formation energy is still similar to the PBE value, the Sn vacancy energy is 0.7 eV more positive despite giving a similar defect structure to the other methods. The PBE results are also similar to those given in the GGA-PW91 study of Togo et al., which were ∼1.17 and ∼1.52 eV for the oxygen interstitial and tin vacancy, respectively. 19920

dx.doi.org/10.1021/jp205148y |J. Phys. Chem. C 2011, 115, 19916–19924

The Journal of Physical Chemistry C

ARTICLE

Figure 3. Relaxed defect structure of (a) a Sn vacancy and (c) an O interstitial in the modeled 3  3  3 supercell, where green indicates the vacancy/ interstitial. The local environments are also shown for the (b) Sn vacancy and (d) O interstitial, where yellow indicates the original position of displaced ions in the pure supercell. Sn and O ions are colored blue and red, respectively.

Table 3. Bond Distances and Atomic Displacements Following the Formation of a Sn Vacancy and O Interstitial in SnO Using PBE, PBE-vdW, and PBE-D2 Methodsa vacancy type

a

distance

PBE

PBE-vdW

PBE-D2

Sn vacancy

O displacements

0.253

0.232

0.226

O interstitial

Sn displacements

0.300

0.318

0.336

Oi
Sn

1.932

1.920

1.920

Oi
Snacross
layer

2.851

2.760

2.786

Sn
Oin
layer

2.099

2.071

2.088

All values are given in Å.

’ DISCUSSION For the calculation of the bulk structure of SnO, the PBE-vdW and PBE-D2 approaches are both seen to give significant improvements over the standard PBE methodology. However, the PBE-vdW results are seen to be in much closer agreement with experiment than those of PBE-D2, which still predicts a distortion in the unit cell. Changes in the minimum energy structure on the introduction of defects are similar for all three methods, although the distortion introduced into the unit cell on minimization of course remains. These results question the validity of the PBE-D2 method for modeling bulk SnO. It is instructive to consider both the form of the damping functions used to apply these corrections and the parameter set used to assess the cause of differences between the methods. The PBE-vdW and PBE-D2 methods use different damping functions in their calculation of the energy correction to account for dispersion interactions, given previously in eqs 7 and 3, respectively. To assess their effects, Figure 4 shows plots of the values of Edisp and the associated force determined using the two approaches for the Sn
Sn interaction as a function of distance. The PBE-vdW energy and force is seen to be less damped than those using the PBE-D2 damping function. In addition, the damping function used by the PBE-D2 method causes an unphysical increase in the damped energy, which gives rise to a sign reversal in the damped force, resulting in it converging to 0.00 eV at low values of Rij. This sign reversal results from the steepness of the damping function, as previously noted by Wu and Yang.37 The damping employed in the PBE-vdW approach of this study is smaller and therefore no sign reversal is seen using this method. Furthermore, for the PBE-D2 approach49

Figure 4. Plots of the Sn
Sn interaction using the (a) PBE-vdW and (b) PBE-D2 methods. Undamped (solid) and damped (dashed) energies (black) and forces (red) are shown, with 0 eV represented by the dotted blue line.

(Figure 4b), the force is seen to change its sign within the range of Sn
Sn distances found in this system (3.539
3.693 Å), leading to this interaction being modeled incorrectly. A van der Waals interaction is an attractive force in a system and therefore the purpose of the dispersion correction is to modify the energy to account for this. By making the force change sign, with the corresponding rise in energy, the damping function is applying an erroneous repulsive force. Grimme states that the need for the damping function is to avoid near-singularities at very small distances.49 In fact, the requirement of the damping function is to account for, and reproduce, the real effect of charge densities overlapping as the value of Rij reduces. As the overlap increases, the polarizability is reduced and hence so is the van der Waals interaction, changing the R
6 dependence.37,74,75 The Cij6 and Rijr parameters used in the two models, as given previously in Table 1, are seen to differ in a number of ways. The Cij6 parameter used in the PBE-D2 approach for the Sn
Sn interaction has a considerably larger dispersion coefficient, 401.20 eV Å6, than those for the Sn
O and O
O interactions, 53.97 and 7.26 eV Å6, respectively. This implies that the Sn ion is considerably more polarizable than the O ion. For the PBE-vdW approach the dispersion coefficients are similar for all three interactions, ranging from 73.07 to 77.34 eV Å6. The cause of the difference is due to the 19921

dx.doi.org/10.1021/jp205148y |J. Phys. Chem. C 2011, 115, 19916–19924

The Journal of Physical Chemistry C origin of the parameters. The PBE-D2 coefficients were determined from atomic polarizabilities suitable for the neutral species. This therefore leads to a large difference between the Sn
Sn and , and O
O terms. The result of this is that the value of the EPBE
D2 disp therefore the van der Waals interactions, will be dominated by the Sn
Sn interactions, with the Sn
O and, particularly, the O
O interactions having a negligible influence. Those used in the PBEvdW method, however, are based on Hartree
Fock-calculated polarizabilities specific to SnO. This results in all three interactions having an importance, as would be expected in a semi-ionic material with ions of similar polarizabilities. Differences are also seen in the values of Rijr . Grimme’s approach involves determining Rijr from the sum of van der Waals radii49 whereas those used in the PBE-vdW method are taken from experimental distances. Although the values for Sn
Sn and O
O interactions are roughly similar for the two approaches, the Sn
O distance is not, which is a consequence of using van der Waals radii. The experimental in-layer Sn
O distance is 2.222 Å, as used in the PBE-vdW method, but that resulting from the van der Waals radii is 3.146 Å. This means that the in-layer Sn
O interactions are strongly damped with the PBE-D2 method. The underperformance of the PBE-D2 method is a direct consequence of the parameter set used, in terms of the large Sn
Sn interaction, and the sign reversal of the damped force. These two factors affect the in-layer and across-layer distances, and therefore the lattice vectors, differently due to the variation in the bond length. For the smaller in-layer distances, the large dispersion correction reduces the bond distances from that predicted by PBE. However, the sign reversal in the force introduces a repulsive component into the damped energy for nearest neighbor Sn
Sn interactions. This results in the Sn
Sn distance not being reduced sufficiently and remaining larger than experiment. This is additionally affected by the value of Rijr for the Sn
O interaction being 0.924 Å larger than the in-layer distance, meaning that any interactions between nearest neighbors in the layer will effectively be fully damped, and therefore the van der Waals interaction will have a significantly reduced influence. For the longer across-layer Sn
Sn distance, although the interaction is still within the region with the sign reversal, it is less repulsive. This means that the excessive Sn
Sn interaction for non-nearest neighbors dominates, causing the distances to reduce from the PBE value by too much, resulting in a Sn
Sn distance shorter than experiment. Overall, these effects result in the PBE-D2 predicted unit cell being larger than experiment in the ab plane, but shorter in the c vector. This also causes a discrepancy in the Sn position, with a u value of 0.2422 rather than the experimental value of 0.2381. For the PBE-vdW method, the magnitude of Cij6 reflects the ionic polarizabilities and the use of experimental bond distances result in excellent reproduction of the experimental structure. The minimized vectors are much closer to experiment, with the percentage difference in the a and c lattice vectors being 0.0 and 0.2%, respectively, giving rise to a c/a ratio 0.2% greater than experiment. The bond distances also show similar errors, with the in- and across-layer distances being within 0.1% and 0.2%, respectively, of the experimental value. For the defect systems, the different methods of applying the dispersion corrections are seen to have a greater effect on the energetics of defect formation than on the structure. Although the structure surrounding a Sn vacancy is similar, the formation energy is considerably larger for the PBE-D2 method than those predicted by either PBE or PBE-vdW. This increased formation energy results from the large Sn
Sn dispersion corrections in the bulk

ARTICLE

SnO material, leading to excessive van der Waals interactions between the Sn atoms. Therefore, the energy required to remove a Sn atom to form the vacancy becomes much larger. Concerns with the atomic-based parameters used in the DFTD2 method have also been raised by other studies for ionic/semiionic systems. For example, Conesa76 recently used it to model the relative stabilities of different TiO2 and Al2O3 polymorphs. Although they used the PBE-D2 method, finding it to perform well with consistent trends to experiment, they rejected the parameter set suggested by Grimme49 for use with oxide materials. This was due to a larger Cij6 value for Ti
Ti than for O
O, suggesting that a titanium cation is more polarizable than an oxygen anion. Instead, they determined ionic-based parameters from the refractive index of the material and used ionic radii instead of van der Waals radii. The DFT-D2 approach has also been used to model a He atom on the MgO (100) surface, with a range of DFT and hybrid-DFT functionals, by MartinezCasado et al.77 Their results showed that the application of this method, while obtaining a qualitatively correct binding of He, fails to determine the binding energy correctly by overbinding it to the surface.

’ CONCLUSIONS In conclusion, this study has assessed the effects of dispersion corrections to the PBE-calculated structures and p-type defects of SnO. Two approaches were used, PBE-D2 and PBE-vdW, allowing for a comparison of atomic- and ionic-based dispersion coefficients, respectively. Although the results show that both approaches give rise to an improved structure over that predicted with PBE, the ionic-based parameters give results in much better agreement with experiment. The atomic polarizabilities and PBE-D2 method lead to a poorly described, distorted unit cell, with a 0.3% contracted c vector and 0.9% expanded a vector relative to experiment. This leads to an associated c/a ratio which is 1.2% smaller than experiment. In comparison, the PBE-vdW method and ionic-based parameter set predict differences to experiment in the a and c vectors of 0.0% and 0.2%. For the defect systems, the PBE-D2 approach is seen to give a formation energy for a neutral Sn vacancy significantly larger than those predicted by PBE and PBE-vdW. The inherent failure of PBE-D2 to model the SnO system is due to the atomic-based parameters used in this approach, which is unsuitable for ionic/semi-ionic systems. The unit cell distortions are caused by an excessive Sn
Sn interaction and very small O
O interaction, coupled to a Sn
O Rij6 term which strongly damps in-layer Sn
O interactions, giving an incorrect representation of dispersion forces. Furthermore, the parameter set used to define the damping function causes an unphysical sign reversal in the force, which, in the case of Sn
Sn, occurs within the range of Sn
Sn distances found in the structure. The excessive Sn
Sn dispersion correction also gives rise to an overestimation of the Sn vacancy formation energy. Overall, this study shows that the use of an ionic-based parameter set, with the PBE-vdW method, offers a significant improvement to the modeling of the layered structure of SnO, which is dictated by van der Waals forces. This may also extend to the simulation of other layered ionic/semi-ionic materials. ’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected]; [email protected]; [email protected]. 19922

dx.doi.org/10.1021/jp205148y |J. Phys. Chem. C 2011, 115, 19916–19924

The Journal of Physical Chemistry C

’ ACKNOWLEDGMENT This work was supported by Science Foundation Ireland through the Principal Investigators program (PI Grant Nos. 06/IN.1/I92 and 06/IN.1/I92/EC07). Calculations were performed on the IITAC and Lonsdale supercomputers as maintained by TCHPC, and the Stoney supercomputer as maintained by ICHEC. ’ REFERENCES (1) Rumyantseva, M. N.; Safonova, O. V.; Boulova, M. N.; Ryabova, L. I.; Gas’kov, A. M. Russ. Chem. Bull. 2003, 52, 1217–1238. (2) Batzill, M.; Diebold, U. Prog. Surf. Sci. 2005, 79, 47–154. (3) Godinho, K. G.; Walsh, A.; Watson, G. W. J. Phys. Chem. C 2008, 113, 439–448. (4) Li, G.-J.; Kawi, S. Mater. Lett. 1998, 34, 99–102. (5) Batzill, M. Sensors 2006, 6, 1345–1366. (6) Choi, Y.-J.; Hwang, I.-S.; Park, J.-G.; Choi, K. J.; Park, J.-H.; Lee, J.-H. Nanotechnology 2008, 19, 095508. (7) Jarzebski, Z. M.; Marton, J. P. J. Electrochem. Soc. 1976, 123, 299C–310C. (8) Klc-, C-.; Zunger, A. Phys. Rev. Lett. 2002, 88, 095501. (9) Ogo, Y.; Hiramatsu, H.; Nomura, K.; Yanagi, H.; Kamiya, T.; Hirano, M.; Hosono, H. Appl. Phys. Lett. 2008, 93, 032113. (10) Ogo, Y.; Hiramatsu, H.; Nomura, K.; Yanagi, H.; Kamiya, T.; Kimura, M.; Hirano, M.; Hosono, H. Phys. Status Solidi A 2009, 206, 2187–2191. (11) Guo, W.; Fu, L.; Zhang, Y.; Zhang, K.; Liang, L. Y.; Liu, Z. M.; Cao, H. T.; Pan, X. Q. Appl. Phys. Lett. 2010, 96, 042113. (12) Liang, L. Y.; Liu, Z. M.; Cao, H. T.; Yu, Z.; Shi, Y. Y.; Chen, A. H.; Zhang, H. Z.; Fang, Y. Q.; Sun, X. L. J. Electrochem. Soc. 2010, 157, H598–H602. (13) Liang, L. Y.; Liu, Z. M.; Cao, H. T.; Pan, X. Q. ACS Appl. Mater. Interfaces 2010, 2, 1060–1065. (14) Hosono, H.; Ogo, Y.; Yanagi, H.; Kamiya, T. Electrochem. Sol. Stat. Lett. 2011, 14, H13–H16. (15) Odani, A.; Nimberger, A.; Markovsky, B.; Sominski, E.; Levi, E.; Kumar, G. G.; Motiei, M.; Gedanken, A.; Dan, P.; Aurbach, D. J. Power Sources 2003, 119, 517–521. (16) Sabnis, A. G. J. Vac. Sci. Technol. 1978, 15, 1565–1567. (17) Moreno, M. S.; Mercader, R. C. Phys. Rev. B 1994, 50, 9875– 9881. (18) Watson, G. W. J. Chem. Phys. 2001, 114, 758–763. (19) Raulot, J.-M.; Baldinozzi, G.; Seshadri, R.; Cortona, P. Solid State Sci. 2002, 4, 467–474. (20) Walsh, A.; Watson, G. W. Phys. Rev. B 2004, 70, 235114–7. (21) Walsh, A.; Watson, G. W. J. Phys. Chem. B 2005, 109, 18868– 18875. (22) Liu, Q.-J.; Liu, Z.-T.; Feng, L.-P. Comput. Mater. Sci. 2010, 47, 1016–1022. (23) Walsh, A.; Payne, D. J.; Egdell, R. G.; Watson, G. W. Chem. Soc. Rev. 2011, DOI: 10.1039/c1cs15098g (24) Waghmere, U. V.; Spaldin, S. A.; Kandpal, H. C.; Seshadri, R. Phys. Rev. B 2003, 67, 125111. (25) Terpstra, H. J.; de Groot, R. A.; Haas, C. Phys. Rev. B 1995, 52, 11690–11697. (26) Watson, G. W.; Parker, S. C. J. Phys. Chem. B 1999, 103, 1258–1262. (27) Watson, G. W.; Parker, S. C.; Kresse, G. Phys. Rev. B 1999, 59, 8481–8486. (28) Walsh, A.; Watson, G. W. J. Solid State Chem. 2005, 178, 1422 1428. (29) Payne, D. J.; Egdell, R. G.; Law, D. S. L.; Glans, P. A.; Learmonth, T.; Smith, K. E.; Guo, J. H.; Walsh, A.; Watson, G. W. J. Mater. Chem. 2007, 178, 1422–1428. (30) Walsh, A.; Watson, G. W.; Payne, D. J.; Edgell, R. G.; Guo, J.; Glans, P.-A.; Learmonth, T.; Smith, K. E. Phys. Rev. B 2006, 73, 235104.

ARTICLE

(31) Sivaramasubramaniam, R.; Muhamad, M. R.; Radhakrishna Phys. Status Solidi A 1993, 136, 215–222. (32) Geurts, J.; Rau, S.; Richter, W.; Schmitte, F. J. Thin Solid Films 1984, 121, 217–225. (33) Errico, L. A. Phys. B Cond. Matt. 2007, 389, 140–144. (34) Togo, A.; Oba, F.; Tanaka, I.; Tatsumi, K. Phys. Rev. B 2006, 74, 195128. (35) Koch, W.; Jan Baerends, E.; Holthausen, M. C. A Chemist’s Guide to Density Functional Theory; Wiley-VCH Verlag GmbH: Weinheim, Germany, 2001. (36) Dobson, J. F.; McLennan, K.; Rubio, A.; Wang, J.; Gould, T.; Le, H. M.; Dinte, B. P. Aust. J. Chem. 2001, 54, 513–527. (37) Wu, Q.; Yang, W. J. Chem. Phys. 2002, 116, 515–524. (38) Scanlon, D. O.; Walsh, A.; Morgan, B. J.; Watson, G. W. J. Phys. Chem. C 2008, 112, 9903–9911. (39) Londero, E.; Schr€oder, E. Phys. Rev. B 2010, 82, 054116. (40) Cora, F.; Patel, A.; Harrison, N. M.; Roetti, C.; Catlow, C. R. A. J. Mater. Chem. 1997, 7, 959–967. (41) Scanlon, D. O.; Watson, G. W.; Payne, D. J.; Atkinson, G. R.; Egdell, R. G.; Law, D. S. L. J. Phys. Chem. C 2010, 114, 4636–4645. (42) Christensen, N. E.; Svane, A.; Peltzer y Blanca, E. L. Phys. Rev. B 2005, 72, 014109. (43) Godby, R. W.; Schl€uter, M.; Sham, L. J. Phys. Rev. Lett. 1986, 56, 2415. (44) Walsh, A.; Da Silva, J. L. F.; Wei, S. H.; Korber, C.; Klein, A.; Piper, L. F. J.; DeMasi, A.; Smith, K. E.; Panaccione, G.; Torelli, P.; Payne, D. J.; Bourlange, A.; Egdell, R. G. Phys. Rev. Lett. 2008, 100, 167402. (45) Kehoe, A. B.; Scanlon, D. O.; Watson, G. W. Phys. Rev. B 2011, 83, 233202. (46) Scanlon, D. O.; Watson, G. W. Phys. Chem. Chem. Phys. 2011, 13, 9667–9675. (47) Dion, M.; Ryberg, H.; Schr€oder, E.; Langreth, D. C.; Lundqvist, B. I. Phys. Rev. Lett. 2004, 92, 246401. (48) Grimme, S. J. Comput. Chem. 2004, 25, 1463–1473. (49) Grimme, S. J. Comput. Chem. 2006, 27, 1787–1799. (50) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. J. Chem. Phys. 2010, 132, 154104. (51) Kerber, T.; Sierka, M.; Sauer, J. J. Comput. Chem. 2008, 29, 2088–2097. (52) Duan, Y. Phys. Rev. B 2008, 77, 045332. (53) Foster, M. E.; Sohlberg, K. Phys. Chem. Chem. Phys. 2010, 12, 307–322.  ngyan, J. G. J. Phys. Chem. A (54) Bucko, T.; Hafner, J.; Lebegue, S.; A 2010, 114, 11814–11824. (55) Adamo, C.; Barone, V. J. Chem. Phys. 1999, 110, 6158–6170. (56) Lide, D. R., Ed. CRC Handbook of Chemistry of Physics; CRC Press: Boca Raton, FL, 2002. (57) Kresse, G.; Hafner, J. Phys. Rev. B 1994, 49, 14251–14271. (58) Kresse, G.; Furthm€uller, J. Phys. Rev. B 1996, 54, 11169–11186. (59) Bl€ochl, P. E. Phys. Rev. B 1994, 50, 17953. (60) Kresse, G.; Joubert, D. Phys. Rev. B 1999, 59, 1758–1775. (61) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865. (62) Murnaghan, F. D. Proc. Natl. Acad. Sci. U.S.A. 1944, 30, 244–247. (63) Monkhorst, H. J.; Pack, J. D. Phys. Rev. B 1976, 13, 5188–5192. (64) Gonzalez, C.; Lim, E. C. J. Phys. Chem. A 2003, 107, 10105– 10110. (65) Wilson, M.; Madden, P.; Peebles, S. A.; Fowler, P. W. Mol. Phys. 1996, 88, 1143–1153. (66) Bernasconi, L.; Wilson, M.; Madden, P. A. Comput. Mater. Sci. 2001, 22, 94–98. (67) Singh, A. K.; Janotti, A.; Scheffler, M.; Van de Walle, C. G. Phys. Rev. Lett. 2008, 101, 055502. (68) Varley, J. B.; Janotti, A.; Singh, A. K.; Van de Walle, C. G. Phys. Rev. B 2009, 79, 245206. (69) Varley, J. B.; Janotti, A.; Van de Walle, C. G. Phys. Rev. B 2010, 81, 245216. 19923

dx.doi.org/10.1021/jp205148y |J. Phys. Chem. C 2011, 115, 19916–19924

The Journal of Physical Chemistry C

ARTICLE

(70) Noguchi, Y.; Takahashi, M.; Miyayama, M. J. Ceram. Soc. Jpn. 2004, 112, 50–56. (71) Bradley, C. J.; Cracknell, A. P. Mathematical Theory of Symmetry in Solids; Oxford Univeristy Press, Oxford, U.K., 1972. (72) Stark, J. G.; Wallace, H. G. Chemistry Data Book; John Murray Ltd.: London, U.K., 1982. (73) Oba, F.; Choi, M.; Togo, A.; Seko, A.; Tanaka, I. J. Phys.: Condens. Matter 2010, 22, 384211. (74) Zimmerli, U.; Parrinello, M.; Koumoutsakos, P. J. Chem. Phys. 2004, 120, 2693–2699. (75) Housden, M. P.; Pyper, N. C. J. Phys.: Condens. Matter 2008, 20, 085222. (76) Conesa, J. C. J. Phys. Chem. C 2010, 114, 22718–22726. (77) Martinez-Casado, R.; Mallia, G.; Usvyat, D.; Maschio, L.; Casassa, S.; Sch€utz, M.; Harrison, N. M. J. Chem. Phys. 2011, 134, 014706.

19924

dx.doi.org/10.1021/jp205148y |J. Phys. Chem. C 2011, 115, 19916–19924