Density Functional Energetics of α-Quartz for Calibration of SiO2

4168

J. Phys. Chem. B 2005, 109, 4168-4171

Density Functional Energetics of r-Quartz for Calibration of SiO2 Interatomic Potentials N. Flocke, Wuming Zhu, and S. B. Trickey* Quantum Theory Project, Departments of Physics and of Chemistry, UniVersity of Florida, GainesVille Florida 32611-8435 ReceiVed: July 12, 2004; In Final Form: NoVember 22, 2004

All-electron, Gaussian basis density functional theory (DFT) total energies for R-quartz were generated on a substantial grid of points for eventual use in a first-principles check of calibration of a widely used effective interionic potential for SiO2. Unlike a previous all-electron, Gaussian basis calculation, we do not find a weak double minimum with respect to cell volume. Rather there is a weak double minimum with respect to internal cell parameters that is at the margins of numerical precision and, hence, may only be the signal of a very flat potential surface. The results show typical dependence on the choice of approximate DFT functional and differ from experimental values enough that both aspects may be problematic for parametrization of potentials.

Motivation and Approach Multiscale material simulations using molecular dynamics (MD) embed a so-called quantum mechanical (QM) region in a classical region. In the QM region, the forces are obtained as gradients of quantum mechanical total electronic energies. The classical region uses a potential function. The fidelity of the classical potential to QM results therefore is crucial. Two physically plausible pair potentials widely used in simulations of SiO2 were calibrated to a combination of first-principles electronic-structure data from a small cluster and experimental data on R-quartz.1-3 The detailed form for these “TTAM” and “BKS” potentials is irrelevant to the present investigation. (TTAM and BKS are commonly used identifiers based on the author names for ref 1 and refs 2 and 3, respectively.) What we study here is whether calculated R-quartz data could replace the calibration to experimental data. Implicitly the issue is whether an all-computed version of the TTAM/BKS potential could be constructed that is as successful as the partly empirical one. Elsewhere we have shown that the portion of the calibration that relies on a single small cluster treated with the restricted Hartree-Fock approximation is not reliable.4 This paper examines the opposite extreme, the variability in computed periodic (extended) system inputs alone. A third paper5 will report the effects on MD calculations from fitting parameter sets to the data of this paper and the data from ref 4. Both TTAM and BKS used R-quartz as the experimental reference for their calibration. From among the sets of parameter fits for their respective cluster calculations, they chose the set that fits experimental R-quartz equilibrium properties best. Physically, the intrinsic interest in R-quartz is motivated by the six parameters that determine the Wigner-Seitz unit cell geometry. These comprise the two hexagonal cell constants (a, c), the bond length and bond angle in one SiO2 unit, the offset of a silicon atom from a cell wall, and the dihedral angle of an SiO2 unit with respect to the cell basal plane. This diverse collection of conformational parameters provides a very demanding grid on which to fit the potential parameters. * Author to whom correspondence should be addressed. Phone: (352) 392-6978. Fax: (352) 392-8722. E-mail: [email protected].

There are at least four prior all-electron, full-potential, density functional theory (DFT) calculations of the equilibrium geometry of R-quartz6-9 as well as a very extensive nonlocal pseudopotential treatment.10 Even those five studies differ notably. There also are large variations among the more extensive set of other pseudopotential DFT calculations. To eliminate one major source of variability, we focus mostly on all-electron results. Experimental data are found in Lager et al.11 For local-density approximation (LDA) calculations, the reported lattice parameters range over 4.840 Å e a e 4.916 Å (T ) 13 K experiment is 4.9020 Å). (Here and throughout we have converted Bohr radii (au of length) to angstroms as 1 au ) 0.529177 Å.) These cell constants correspond to 1.0988 e c/a e 1.114 (T ) 13 K experiment is 1.1016), with the corresponding range of 36.22 Å3 e Ω e 37.64 Å3 for the cell volume per formula unit (T) 13 K experiment is 37.5 Å3). Catti et al.8 reported a double minimum at volumes of 36.5 and 36.9 Å3, separated in energy by only about 100 µhartree (∼2.7 meV) per formula unit. For gradient-dependent DFT approximations, less has been reported, but the Perdew-Burke-Ernzerhof (PBE) GDA12 gives Ω ) 39.42 Å3 versus Ω ) 40.19 Å3 for PW9113 and Ω ) 37.675 Å3 for BLYP.14 Catti et al. report that the hybrid B3LYP functional15 also gives a double minimum: Ω ) 37.9 Å3 and 38.5 Å3. The variation among the calculations is obvious and larger than might be expected from experience regarding the effect of simple methodological differences on computed crystalline properties. Without that striking variation, we would have simply redone the calculation to get energies on a grid of geometries around equilibrium and proceeded to fitting the potentials. With such variation, a more detailed restudy seemed advisable. Methodology Our approach also is via DFT, with both the local and the gradient-dependent approximations for the exchange-correlation energy Exc and its corresponding functional derivative Vxc. We used the Perdew-Burke-Ernzerhof (PBE) gradient-dependent approximation12 and both the Hedin-Lundqvist (HL)16 and the Perdew-Zunger (PZ)17 local models. We solved the Kohn-Sham equations in a Gaussian-type basis set, with variational charge fitting methodology to ac-

10.1021/jp0469158 CCC: $30.25 © 2005 American Chemical Society Published on Web 02/12/2005

DFT Energetics of R-Quartz

J. Phys. Chem. B, Vol. 109, No. 9, 2005 4169

TABLE 1: KS, Q, and XC Basis Exponents (inverse au2) and Contraction Coefficients (KS only)a type

KS exponent

KS contr. coefficient.

Si s

0.26737389d+05 0.40763786d+04 0.95329315d+03 0.27458441d+03 0.90681987d+02

0.0010354 0.0077158 0.0376661 0.1368397 0.2433422

Q exponent

1.00000000d+03 2.60000000d+02 7.20000000d+01 2.40000000d+01 8.10000000d+00 2.70000000d+00 0.90681987d+02 0.0924999 0.90000000d+00 0.33529057d+02 0.4327388 0.35000000d+00 0.13457629d+02 0.1933028 0.13000000d+00 0.40506916d+01 1.0000000 0.14841839d+01 1.0000000 0.27044323d+00 1.0000000 0.09932161d+00 1.0000000 Si p 0.16373270 d+03 0.0118853 6.00000000d-01 0.38352081d+02 0.0805337 1.80000000d-01 0.12021114d+02 0.2705100 0.41845733d+01 0.4640613 0.14827170d+01 1.0000000 0.33498569d+00 1.0000000 0.09699384d+00 1.0000000 Si d 0.45000000 d+00 1.0000000 a

XC exponent

type

KS exponent

KS contr. coefficient.

3.60000000d+02 9.00000000d+01 2.70000000d+01 9.00000000d+00 3.00000000d+00 1.00000000d+00 0.36000000d+00 0.17000000d+00 0.08000000d+00

Os

0.31195600d+05 0.46693800d+04 0.10626200d+04 0.30142600d+03 0.98515300d+02 0.35460900d+02

0.0002080 0.0016183 0.0083949 0.0340909 0.1094802 0.2664469

6.00000000d-01 1.80000000d-01

Op

Od

Q exponent

9.00000000d+02 2.40000000d+02 7.00000000d+01 2.20000000d+01 7.00000000d+00 2.40000000d+00 1.10000000d+00 0.13617900d+02 1.0000000 0.50000000d+00 0.53861800d+01 1.0000000 0.22000000d+00 0.15387300d+01 1.0000000 0.65000000d+00 1.0000000 0.30000000d+00 1.0000000 0.14000000d+00 1.0000000 0.11486300d+03 0.0023269 6.00000000d-01 0.26876700d+02 0.0182863 1.80000000d-01 0.83207700d+01 0.0797585 0.29723700d+01 0.2176610 0.11284800d+01 1.0000000 0.42360000d+00 1.0000000 0.18000000d+00 1.0000000 0.30000000d+00 1.0000000

XC exponent 9.00000000d+02 2.40000000d+02 7.00000000d+01 2.20000000d+01 7.00000000d+00 2.40000000d+00 1.10000000d+00 0.50000000d+00 0.22000000d+00

6.00000000d-01 1.80000000d-01 2.00000000d+00 0.60000000d+00

Functions for Q and XC denoted as p type actually are p2 except the last two for O are pz; see text.

celerate the calculation. This “linear combination of Gaussiantype orbitals with fitting functions” methodology is embodied in the code GTOFF.18 The method is all-electron (no pseudopotentials) and imposes no shape approximations on the potential. Three Gaussian-type basis set expansions are used: the “KS” basis for the orbitals, the “Q” basis for the electron density n(r), and the “XC” basis both for the exchangecorrelation kernel that depends nonlinearly on n(r) in Exc and for Vxc(r) itself. Because of both basis set sensitivity issues and experience, which shows that high-precision calculation of periodic systems requires somewhat different basis sets from those common in molecular studies, we give the orbital exponents and the contraction coefficients (used only in the KS basis) in Table 1. The two KS basis sets (12s7p1d/6s4p1d for Si and 12s7p1d/ 7s4p1d for O) are crystalline adaptations by Boettger,19,20 of published basis sets.21,22 This basis is somewhat richer and more diffuse than that used by Catti et al.8 Superficially the KS basis is unbalanced in favor of oxygen. In fact, this is appropriate because the rather ionic character of the Si-O interaction depletes electrons from the silicon. Two examples make the point. First, the molecular cluster associated with the BKS potential yields effective Si charges of 2.5 f 2.75 qe depending on basis set and symmetries.4 Second, both the BKS and the TTAM potentials use an effective Si charge of 2.4 qe.1-3 The Q and XC fitting bases also are from Boettger, modified by dropping the unneeded Si f functions and adding two p2 functions (literally, the square of p functions). For O, the XC basis was supplemented further with two pz functions. The result is a 9s2p2 Q basis for both Si and O, a 9s2p2 XC basis for Si, and a 9s2p22pz XC basis for O. The Brillouin zone (BZ) scan was 6 × 6 × 6 in the full zone corresponding to the direct lattice hexagonal unit cell. This grid, in combination with the variational charge fitting, determines the Fermi energy to a fractional error in the integrated total charge of 3.3 × 10-15. GTOFF does the numerical integrals needed in the XC fitting procedure via a spherical grid around each nuclear site and an interstitial grid. To improve the precision of the interstitial integrals, two empty sphere grids were added. These were centered on the c-axis, halfway from the cell midplane to the z-coordinate of the uppermost and

lowermost O atoms. The grid in each sphere and interstitially is somewhat complicated; a full description is in ref 18. Energy minimization was done by numerical gradient extrapolation. Starting at the T ) 296 K experimental geometry, Etot per cell for the HL local spin density approximation (LSDA) was computed on a grid of eight regularly spaced points in the internal coordinates around that geometry (i.e., at fixed c, a). In terms of initial values for the internal fractional displacement ui, Si-O bond length d, O-Si-O angle θi, and SiO2 plane tilt angle φi, the internal grid at fixed c and a was ui ( 0.01, di ( 0.1 au (di ( 0.0529 Å), θi ( 1°, and φi ( 1°. The resulting total energies ranged over a few mhartree (out of about 1300 hartree per unit cell depending on the particular XC approximation). New values of those four internal parameters were obtained by interpolation, and then new energies were computed on a four-point grid in a and c around the original ai and ci, namely, ai ( 0.1 au, ci ( 0.1 au (ai ( 0.0529 Å and ci ( 0.0529 Å). These were extrapolated to new values of a and c, and another eight-point internal coordinate grid was calculated. Three such cycles, with declining increments, yielded the minimum energy HL LDA configuration (as well as the energy vs conformation data for eventual fitting of potential parameters). The precision is about (2 µhartree in total energy per cell. We did not do further optimization of the structure in the PBE XC approximation since our main objective was to generate a reliable set of energies on a grid of geometries for potential fitting. Therefore the lowest PBE energy reported here is not guaranteed to be the PBE minimum but only the lowest-energy point on the grid searched. Results Table 2 compares our results for the equilibrium structure with previous all-electron calculations6-9 and, for reference, the most recent pseudopotential calculation.10 Our PZ and HL results are essentially identical, so we report only the latter. Note that these values are at an actual calculated data point, not interpolated. Reference 6 only reports the lattice parameters at experimental values but gives an interpolated minimum cell volume. Reference 8 gives a double minimum. For the B3LYP hybrid DFT functional, we list the lower energy structure they reported. For LDA, we give both.

4170 J. Phys. Chem. B, Vol. 109, No. 9, 2005

Flocke et al.

TABLE 2: Calculated Cell Parameters (Å), c/a, Si Displacement Parameter u (As a Fraction of a), Si-O Bond Length d (Å), O-Si-O Angle θ (deg), Angle O of SiO2 Plane Relative to Cell Basal Plane (deg), and Volume per Formula Unit Ω (Å3) for r-Quartza source expt present present second ref 6 ref 7 ref 8 ref 8 second ref 10 present ref 7 ref 10 ref 9 ref 8 a

DFT model HL HL HL PW92 SVWN SVWN PZ PBE PBE PW91 BLYP B3LYP

a

c

c/a

4.9020 4.9313 4.9313 4.9160

5.4000 5.3606 5.3606 5.4054

4.8400 4.8700 4.8992 4.9210

5.3940 5.386 5.3832 5.4215

5.0271 4.9138 4.9550

5.5089 5.4052 5.4280

1.1016 1.0871 1.0871 1.10996 1.104 1.1145 1.1060 1.0988 1.1017 1.1099 1.0958 1.1000 1.0955

θ

u

d

0.0331 0.0326

1.613 1.628 1.632

108.836 108.827

24.213 24.218

0.0321

1.626 1.630 1.602 1.638

109.43 108.836

24.213

1.615

109.28

1.636

φ

Ω 37.5 37.63 37.63 37.64 36.22 36.5 36.9 37.30 37.90 39.42 40.2 37.68 38.5

See text for comments regarding present PBE and present second HL results. Experimental data (T ) 13 K) are from ref 11 as quoted in ref

8.

Typically (though not always) LDA contracts internuclear distances relative to experiment. Our calculations give a slightly lengthened a (0.6%) and slightly shortened c (0.7%) relative to experiment. The molar volume thus is enlarged relative to experiment by 0.3%, and the calculated c/a is low by 1.3%. These results are reasonably consistent with ref 6. Those authors used full-potential, linearized augmented plane wave (FLAPW) methodology and the same LDA to find essentially the same a as ours but a slightly longer (0.8%) c and a correspondingly slightly higher molar volume. Given the complete independence of the two calculations, it is clear that the LDA prediction of the equilibrium R-quartz a and c parameters is reasonably wellsettled. For the HL LDA approximation, our calculations yield a second set of internal cell parameters at the same a and c that has an energy about 65 µhartree per cell above the calculated minimum or 45 parts in 109. (The corresponding PZ difference is 85 µhartree per cell.) That slightly unfavored configuration is listed as “present second” in Table 2. The energy difference between the two configurations is only 1 order of magnitude larger than our estimate of the limits of numerical precision. Even that may be deceptive, since precise determination of the internal coordinates is quite difficult. To illustrate, at a and c corresponding to the HL minimum, our calculations are unable to distinguish two internal configurations: u ) 0.0331, d ) 1.628 Å, θ ) 108.836°, φ ) 24.213° and u ) 0.0321, d ) 1.628 Å, θ ) 108.827°, φ ) 24.218°. These two differ by 2 µhartree per cell, i.e., less than 2 parts in 109. The double minimum we find differs qualitatively, however, from the one reported by Catti et al.8 They found significant lattice parameter shifts between their two minima; compare the listing for ref 8 and ref 8 second. As discussed just above, this finding is inconsistent with both our results and the methodologically completely independent result of ref 6. In part, we agree with the assessment of Catti et al. that the energy surface is sufficiently flat around the minimum that nearly identical calculated minima are possible. Sorting among such minima is generally agreed to be at the margins of current electronicstructure calculations (particularly in view of the differences introduced by the various Exc models). We disagree implicitly, however, on the physics of the situation. In our calculation, the cell parameters themselves are uniquely determined, and the soft modes of motion are internal to the cell. We suggest that such a picture is more consistent with the brittleness of quartz than a potential surface that is flat with respect to a and c. Our calculations further suggest that basis set size, BZ scan density,

and convergence tolerances may have led to the Catti et al. results. For example, we found it essential to use both much tighter and much more diffuse functions than apparently were used in ref 8. Again, since our focus is on input to potential parametrization, we did not pursue the analysis of causes. The calculated cohesive energy with respect to separated atoms for HL LDA is 21.74 eV compared to the experimental value (as given in ref 23) of 19.23 eV. Reference 6 gives 24.57 eV for the same LDA. We get about that value if the decontracted crystalline basis sets are used to calculate the atomic total energies. But it is well-known that such basis sets are too constricted to describe the atom properly. Because ref 6 used the FLAPW basis, it is not clear whether such constriction in their basis set causes the difference. A finding consistent with that suggestion is the Liu et al.23 ultrasoft pseudopotential calculation with a slightly different LDA; it gets 22.42 eV. Consistent with general experience on other materials, the PBE GGA gives lattice parameters a and c that are slightly too long with respect to experiment. Within the range of lattice and internal parameters tested, we found no double minimum for PBE. Recall, however, our earlier remarks on the limitations of our PBE search. In particular, we do not claim that the values shown in Table 2 are the optimized equilibrium ground-state parameters for R-quartz from PBE. Rather, they are the values corresponding to the lowest energy point (in the PBE model) of a grid of conformations generated in searching for the energy minimum in LDA. It is interesting that our PBE results and the Gaussian basis BLYP results of ref 9 compare very favorably. In contrast, the PW91 functional with pseudopotentials gives values of a and c that are far too large even though c/a comes out rather well. Because the fitting function methodology avoids all fourcenter integrals, GTOFF does not include hybrid functionals. Therefore we cannot compare with the B3LYP calculation presented in ref 8. Table 3 provides a comparison of the calculated equilibrium Si-O bond length and O-Si-O angle for these DFT calculations on R-quartz versus the geometry-optimized and “tetrahedrally constrained” H4SiO4 cluster results from highly refined molecular calculations4 and the H4SiO4 results published by BKS. Tetrahedrally constrained means that the SiO4 core was held in tetrahedral symmetry. “CCSD(T)” stands for coupled cluster single and double excitations with perturbative (noniterative) triple excitations, essentially the state of the art in quantum chemistry today for treatment of electron-correlation

DFT Energetics of R-Quartz

J. Phys. Chem. B, Vol. 109, No. 9, 2005 4171

TABLE 3: Comparison of Si-O Bond Length (Å) and O-Si-O Angle (deg) from the Present r-Quartz Calculations with Experiment and Cluster Calculationsa calculation

Si-O length

angle O-Si-O

R-quartz (PBE) R-quartz (HL) R-quartz (exptl) CCSD(T)/aug-cc-pVDZ CCSD(T)/aug-cc-pVDZ (tc) BKS

1.638 1.628 1.613 1.675 1.6806 1.63

108.836 108.836 100.3 109.47 109.47

T ) 13 K in the R-quartz experiment.11 The optimized and tetrahedrally constrained (tc) BKS clusters H4SiO411 and the H4SiO4 are as reported in the BKS papers.2,3 See text. a

contributions. “Aug-cc-pVDZ” is a Gaussian basis set description that for the present purposes can be interpreted simply as quite large. Table 3 illustrates an intrinsic difficulty in the logic of the BKS/TTAM approach. Irrespective of whether the R-quartz reference data are experimental or computed and irrespective of whether the cluster data are for the fully optimized or tetrahedrally constrained H4SiO4, the crystalline Si-O distance is much shorter than that for the cluster. But the published parametrization procedure is to use the cluster to generate multiple parameter sets and then select from among those sets according to best fit to the crystal. Clearly the procedure will not work unless the cluster and crystal Si-O bond lengths are close, a condition that obviously is not met if the calculations are done to high accuracy. Table 3 also illustrates how coincidence can influence potential fitting. If one had decided to do a low-level (restricted Hartree-Fock) cluster calculation (as BKS did) and combine it with an LDA study of quartz, the result would, at least superficially, have been consistent. Notice the bond lengths and angles in Table 3. But combining an improved cluster calculation with a DFT crystal would not work. Summary These calculations show that the DFT (LDA) values of the ground-state cell parameters of R-quartz differ slightly from the T) 13 K experimental data. Consistent with ref 6, we find no sign of the double minimum in the a and c parameters found by ref 8. There is typical dependence on the choice of DFT exchange-correlation approximation. The calculated Si-O R-quartz bond length is substantially shorter than those in small clusters that have SiO2 local bonding features. All of these features suggest that an “all-computed” potential of the BKS/ TTAM form would be difficult if not impossible to achieve by using the best available calculations for molecules together with those for crystals. Separately we will present a study of the effects of calibration to various choices of first principles calculated data.5

Acknowledgment. We thank J. C. Boettger for advice and assistance, particularly about the basis sets. We thank K. Runge, R. J. Bartlett, H.-P. Cheng, F. E. Harris, J. H. Simmons, and S. Yip for helpful discussions. This work was supported by the U.S. National Science Foundation under KDI Award DMR9980015 and ITR Award DMR-0325553. References and Notes (1) Tsuneyuki, S.; Tsukada, M.; Aoki, H.; Matsui, Y. Phys. ReV. Lett. 1988, 61, 869. (2) van Beest, B. W. H.; Kramer, G. J.; van Santen, R. A. Phys. ReV. Lett. 1990, 64, 1955. (3) Kramer, G. J.; Farragher, N. P.; van Beest, B. W. H.; van Santen, R. A. Phys. ReV. B 1991, 43, 5068. (4) Al-Derzi, A. R.; Cory, M. G.; Runge, K.; Trickey, S. B. J. Phys. Chem. A, in press. (5) Zhu, Wuming; Taylor, C.; Al-Derzi, A. R.; Runge, K.; Trickey, S. B.; Zhu, T.; Yip, S. Unpublished work. (6) Di Pomponio, A.; Continenza, A. Phys. ReV. B 2003, 48, 12558. (7) Zupan, A.; Blaha, P.; Schwarz, K.; Perdew, J. P. Phys. ReV. B 1998, 58, 11266. (8) Catti, M.; Civalleri, B.; Ugliengo, P. J. Phys. Chem. B 1999, 104, 7259. (9) Gnani, E.; Reggiani, S.; Colle, R.; Rudan, M. IEEE Trans. Electron DeVices 2000, 47, 1795. (10) Demuth, Th.; Jeanvoine, Y.; Hafner, J.; A Ä ngya´n, J. G. J. Phys.: Condens. Matter 1999, 11, 3833. (11) Lager, G. A.; Jorgensen, J. D.; Rotella, F. J. J. Appl. Phys. 1982, 53, 6751. (12) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. ReV. Lett. 1997, 77, 3865. Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. ReV. Lett. 1997, 78, 1396 (erratum). (13) (a) Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Phys. ReV. B 1992, 46, 6671. (b) Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Phys. ReV. B 1993, 48, 4978 (erratum) and references therein. (14) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (15) Stephens, P. J.; Devlin, J. F.; Chablowski, C. F.; Frisch, M. J. J. Phys. Chem. 1994, 98, 11623. (16) Hedin, L.; Lundqvist, B. I. J. Phys. C 1971, 4, 2064. (17) Perdew, J. P.; Zunger, A. Phys. ReV. B 1981, 23, 5048. (18) Trickey, S. B.; Alford, J. A.; Boettger, J. C. In Computational Materials Science, Theoretical and Computational Chemistry; Leszczynski, J., Ed.; Elsevier: Amsterdam, 2004; Vol. 15, p 171 and references therein. Boettger, J. C.; Trickey, S. B. J. Mol. Struct. (THEOCHEM) 2000, 50102, 285. (19) Boettger, J. C Int. J. Quantum Chem. 1996, 60, 1345. (Note that the particular number of the journal is mispaginated inside the issue; there this paper appears as “S 30, 133”.) (20) Boettger, J. C.; Ray, A. K. Int. J. Quantum Chem. 2000, 80, 824. (21) Huzinaga, S. Approximate Atomic Functions II, Department of Chemistry, University of Alberta. Unpublished work, 1971. (22) van Duijneveldt, F. B. IBM Research Report RJ945. Unpublished work, 1971. (23) Liu, F.; Garoflini, S. H.; King-Smith, D.; Vanderbilt, D. Phys. ReV. B 1994, 49, 12528.