J. Phys. Chem. 1996, 100, 16105-16108
16105
Local Density-Functional Polarizabilities and Hyperpolarizabilities at the Basis-Set Limit Ross M. Dickson† and Axel D. Becke* Department of Chemistry, Queen’s UniVersity, Kingston, Ontario K7L 3N6, Canada ReceiVed: February 27, 1996; In Final Form: July 28, 1996X
We have performed finite-field calculations of the dipole polarizabilities and hyperpolarizabilities of some first-row compounds in the local spin-density approximation (LDA) at the basis-set limit using the numerical density-functional code NUMOL. We show that the calculations of Guan et al. [J. Chem. Phys. 1993, 98, 4753] on dipole moments and polarizabilities are well converged with their high-quality double-zeta basis set which includes field-induced polarization functions. However, hyperpolarizabilities are not similarly well converged. We provide reference values for first hyperpolarizabilities and estimated LDA values of some second hyperpolarizabilities.
1. Introduction Electrical properties such as multipole moments and polarizabilities of molecules are of interest in a number of fields, most notably in nonlinear optics and the theory of molecular interactions. The ab initio computation of such properties is a challenging enterprise in which electron correlation, basis set demands, and vibrational and solvent effects all play significant roles.1-3 Density-functional theory is viewed these days as a practical alternative to traditional treatments of electron correlation. It is especially useful in medium and large systems where postHartree-Fock methods become prohibitively expensive. But most density-functional determinations of electrical properties are carried out with basis-set expansions whose convergence must be laboriously verified for each new property calculated.4-13 Therefore, in the present work we given multipole moments, static polarizabilities, and hyperpolarizabilities for a few benchmark systems using our unique basis-set-free density-functional code NUMOL.14,15 While the full time-dependent density-functional theory necessary to calculate frequency-dependent polarizabilities is available,16 calculations of static polarizabilities are possible with the conventional time-independent theory. Coupled-perturbed schemes have been implemented,8 but finite-field schemes require only the most trivial of modifications to existing computer codes and in the present work we employ the latter. We calculate quantities accessible with a constant electric field F in the local spin-density approximation (LDA). We use the “traditional” symbols for the various terms in the electric field expansion of the energy17
E ) E(0) - ∑ µiFi i
1
1
βijkFiFjFk ∑ RijFiFj - 6 ∑ 2 ij ijk 1 24
γijklFiFjFkFl - ... ∑ ijkl
or, equivalently, the dipole moment:
µi ) µi(0) + ∑ RijFj + j
1
1
γijklFjFkFl + ... ∑ βijkFjFk + 6 ∑ 2 jk jkl
† Present address: Hypercube, Inc., 419 Phillip St., Waterloo, Ontario N2L 3X2, Canada. X Abstract published in AdVance ACS Abstracts, September 1, 1996.
S0022-3654(96)00596-5 CCC: $12.00
We also calculate first and second moments of the charge distribution. The first or dipole moment is calculated:
µi ) ∫ F(r)ri d3r where summation over nuclear charges is implicitly included in the integration. With regard to the second moment there are different conventions in use; we prefer the Cartesian convention of Dykstra et al.,1 Viz.:
Qij ) ∫ F(r)rirj d3r These can be easily converted to the more familiar traceless quadrupole moments Θ17 to facilitate comparison with experiment:
Θij ) 1/2 ∫ F(r)(3rirj - r2δij) d3r Atomic units are used throughout; see refs 1, 2, or 7 for conversion factors to SI or CGS units. 2. Procedures All calculations have been carried out with the fully numerical density-functional code NUMOL14,15 using the local correlation functional of Perdew and Wang.18 Experimental equilibrium geometries are used in all calculations and are as follows: H2, r ) 1.401a0; N2, r ) 2.075a0; O2, r ) 2.282a0; HF, r ) 1.733a0; CO, r ) 2.132a0; H2O, r ) 1.809a0, and ∠HOH ) 104.51°; NH3, r ) 1.912a0 and ∠HNH ) 106.7°; CH4, r ) 2.052a0. The molecules have been oriented with respect to the Cartesian axes so that the highest symmetry axis coincides with the z direction for all but CH4, which is inscribed in a cube in the usual manner. H2O has the hydrogens lying in the xz plane. Signs are fixed by the convention that the permanent dipole moment (if any) is directed along the positive z-axis. O2 is the only open-shell molecule studied: The 3∑ ground state is represented by occupying each of the two frontier π* orbitals with one spin-up electron. The spatial orbitals are unrestricted. We follow loosely the polynomial fitting method of Guan et al.7 We choose, as they do, to work from the dipole moment calculated as the expectation value of r rather than taking explicit derivatives of the energy, since we thereby reduce the order of numerical differentiation required. This introduces no error since the Hellmann-Feynman theorem is satisfied by self© 1996 American Chemical Society
16106 J. Phys. Chem., Vol. 100, No. 40, 1996
Dickson and Becke
TABLE 1: Polynomial Fits of µ(Fx) for Hydrogen Fluoridea Standard Mesh, 0 e Fx e 0.010 (11 Data) c0
c2
0.706427 0.706427 0.706427
-1.160 -1.152 -1.152
c4
c6
R
-78.4 -76.4
-54.5
c1 5.9351 5.9225 5.9226 5.9226
c3
c5
c7
191.1 187.6 187.7
2.9 × 104 2.8 × 104
8.8 × 106
7.8 × 10-8 3.8 × 10-10 3.9 × 10-12 3.4 × 10-5 1.4 × 10-7 1.2 × 10-9 1.2 × 10-11
Standard Mesh, 0 e Fx e 0.024 (13 Data) c0
c2
c4
0.706430 0.706427 0.706427
-1.197 -1.149 -1.153
-90.6 -73.5
c1 5.9981 5.9193 5.9229 5.9226
c3
c5
210.9 184.2 188.2
3.9 × 104 2.5 × 104
c6
R
2.0 × 104
2.9 × 10-6 1.0 × 10-7 6.8 × 10-9
1.5 × 107
5.0 × 10-4 1.5 × 10-5 8.6 × 10-7 5.9 × 10-8
rank of the least-squares design matrix for a given range of field strengths.19 (2) We plot the residuals and estimate visually whether an increase in polynomial order would simply amount to fitting “noise” rather than “signal”. We use these two diagnostics to decide when to terminate the sequence of increasing-order fitting polynomials, or to determine whether modification of the range of field strengths is appropriate. Hydrogen fluoride was used to test the sensitivity of the results to (i) the range of field strengths, (ii) the order of fitting polynomial, and (iii) the number of mesh points used in the NUMOL calculations. To test the range of field strengths, we took evenly-spaced fields ranging from 0 to either 0.010 or 0.024 au, comprising 11 and 13 data points, respectively. To test the effect of the molecular mesh, we performed the calculations first with the standard NUMOL mesh of 40 × 194 points on first-row atoms and 20 × 110 points on hydrogens. [The notation X × Y points indicates X radial points or spheres about the nucleus, and Y angular points or rays.] Then we used a fine mesh of 100 × 302 points on first-row atoms and 60 × 194 points on hydrogen. The resulting polynomials are given in Tables 1 and 2. We also tabulate the standard error R as a measure of goodness of fit:
Fine Mesh, 0 e Fx e 0.007 (9 Data) c0
c2
c4
0.706437 0.706437 0.706437
-1.155 -1.145 -1.140
-221.2 -567.8
c1 5.9284 5.9221 5.9221
c3
c5
189.5 189.4
1697
c6
R
4.9 × 106
5.9 × 10-8 1.9 × 10-8 1.0 × 10-8 1.1 × 10-5 4.5 × 10-8 4.8 × 10-8
a Even-order fits are to µ and so β z zxx ) 2c2. Odd-order fits are to µx and so Rxx ) c1. “Standard” and “fine” meshes are defined in the text. R is fit standard error.
consistent solutions to the Kohn-Sham equations; see, for instance, the appendix to ref 7. Consequently, the polynomials employed are N
µ ) ∑ cnFn n)0
and we have µ ) c0, R ) c1, β ) 2c2, and γ ) 6c3. We observe that in most cases, two-sided finite differences with field strengths of 0.001 au give values of R essentially identical to those obtained by fitting. Three-point finite differences with larger fields, (0.01 au, give values of β that vary by at most 1% from values obtained by fitting. We vary from the procedure described by Guan et al., though, in two ways: (1) We employ a singular-value decomposition least-squares scheme which allows us to diagnose the effective
n
R2 )
(yifit - yidata)2 ∑ i)1 n-m
where n - m is the number of degrees of freedom, i.e., the number of data minus the number of fit parameters. Convergence with respect to polynomial order is fairly clear in Tables 1 and 2, as is the small dependence on data range and the larger dependence on the molecular mesh. µ and R are stable to within parts per thousand. βzxx varies about 1% between the standard and fine meshes. Variation with data range depends on the magnitude of the higher-order coefficients, which may be different in other molecules, in which case we adjust the data range to ensure that higher-order terms do not contaminate the coefficients of interest. In the fit of µx(Fx) in Table 1 it would appear that γxxxx ) 6c3 is reasonably stable with respect to variations in fitting range and molecular mesh, so we are able to determine LDA values for this property in certain molecules (such as HF). µz(Fz) was not, however, smooth enough to allow us to estimate γzzzz for N2, nor were we able to obtain any γ for NH3. Uncertainties in the last digit are given in parentheses in the tables when the uncertainty is greater than (1 in the last digit shown. It should be emphasized that these are estimates of the error arising from the fitting procedure; they are not expected to embrace experimental values. Neither do we explicitly estimate deviations from the LDA limit due to mesh inadequacy,
TABLE 2: Polynomial Fits of µz(Fz) for Hydrogen Fluoride c0
c1
c2
c3
c4
c5
c6
Standard Mesh, -0.01 e Fz e 0.01 (21 Data) 0.706229 0.706427 0.706427 0.706427 0.706427 0.706427
6.8587 6.8587 6.8520 6.8520 6.8520 6.8520
0.706381 0.706437 0.706437 0.706437 0.706437 0.706437
6.8543 6.8543 6.8515 6.8515 6.8515 6.8515
-5.3959 -5.3959 -5.3764 -5.3764 -5.3765
103.2 103.2 102.1 102.1
-207.9 -207.9 -207.1
8531 8531
-5723
1.9 × 10-4 1.9 × 10-5 2.1 × 10-7 4.5 × 10-8 1.8 × 10-10 9.0 × 10-11
-4 × 105
7.9 × 10-5 4.5 × 10-6 3.5 × 10-8 3.1 × 10-9 3.2 × 10-10 2.6 × 10-11
Fine Mesh, -0.006 e Fz e 0.006 (11 Data) -5.3821 -5.3821 -5.3740 -5.3740 -5.3742
102.6 102.6 102.2 102.2
-231.8 -231.8 -211.4
7961 7961
R
Local Density-Functional Polarizabilities and Hyperpolarizabilities
au. The experimental value is larger still, 0.66 au22 but includes vibrational and dispersion effects of unknown size. First dipole hyperpolarizabilities (β) are displayed in Table 6. Here the NUMOL results are significantly different from the basis-set results. It is evident from the work of Lee and Colwell that basis-set incompleteness usually has a much larger effect on β than the use of gradient-corrected functionals,11 so we fell that the use of the LDA here is probably a much smaller source of error than the neglect of vibration and dispersion in comparing with experiment. The second-harmonic-generation experiments of Ward and coworkers provide orientational averages of β,24,25
TABLE 3: Dipole Moment Magnitudes (in au)
a
molecule
numericala
DMolb
exptc
CO HF H2O NH3
0.0907 0.7064 0.7300 0.601
0.099 0.699 0.727 0.604
0.048 0.707 0.727 0.579
J. Phys. Chem., Vol. 100, No. 40, 1996 16107
Present work. b Reference 7, basis DNP+. c Collected in ref 20.
although we feel that the present work more nearly represents that limit than any existing basis-set calculations. We draw the reader’s attention to the hydrogen fluoride test case, where we feel the mesh dependence is not atypical.
β h ) 3/5(βzxx + βzyy + βzzz)
3. Results Permanent dipole moments are displayed in Table 3 and compared with the LDA calculations of Guan et al.7 and with experiment. It is clear that the DMol calculations are well converged for this basic property. Agreement with experiment is also quite good, as has been widely reported already. Quadrupole moments (Table 4) are also in excellent agreement with the large-basis LDA calculations of Duffy, Chong, and Dupuis.13 Agreement with experiment is also very good, with the exception of O2 where we obtain a negative sign for Θzz while the experimental value is positive. Multireference CI calculations give a value of -0.185 for this quantity, and all other calculations cited in ref 20 agree that this value is small but negative, so it would seem the experimental value of Θzz for O2 is in need of reevaluation. Dipole polarizabilities (R) are displayed in Table 5. The agreement between the basis-set-free and the DMol results is also quite good for this property. Experimental results are taken from a variety of different experiments: molecular beam studies, Rayleigh scattering, and dielectric constants.7,21,22 Individual tensor components are not directly accessible by experiment, but rather the mean dipole polarizability R j,
which show notable qualitative agreement with our calculations. But Adamowicz and Bartlett, in a numerical Hartree-Fock study of hydrogen fluoride, calculated that vibrational averaging was sufficient to reverse the sign of βzzz and increase γzzzz by almost an order of magnitude,26 and so the present agreement with experiment should not be viewed as anything more than accidental. The β values for NH3 are noteworthy in light of the factor of 2 discrepancy with the DMol calculations of βzzz.7 NUMOL calculations with standard and fine meshes yield similar results, and the error estimates in Table 6 embrace both these. We therefore feel that these numbers give a true representation of the LDA β for ammonia. We note that Guan et al. obtain LDA estimates for βzzz ranging from 5.76 to 82.6 au with various basis sets.7 More recent calculations on CH4 with a much larger basis have βxyz ) 7.52 au,10 in contrast to -13.9 au shown in Table 6 calculated with a double-zeta polarized basis, and in excellent agreement with the present results (-7.5 au). The sign is effectively arbitrary. In Table 7 we give second dipole hyperpolarizabilities (γ) for those molecules for which we could reasonably deduce them from our data. Although the data in Tables 1 and 2 shows remarkable mesh stability for γ in hydrogen fluoride, we find it is generally more sensitive to both mesh and fitting procedure than R or β. The data in Table 7 should not be considered definitive LDA values to better than (10%. The numerical SCF solutions for NH3 were not sufficiently smooth with respect to field variation to yield values of γ, and for N2 only γxxxx could be determined, with error bounds of at least (20%. These are not complete specifications of γ; there are other independent components in the complete tensors necessary for a comparison with experiment. We have only calculated those which appear in the polynomial fit analysis for R. Others would require additional NUMOL runs for additional data points. For H2, however, Bishop, Pipin, and Cybulski27 have carried out definitive theoretical calculations of individual components of γ which are of sufficient accuracy to challenge experimental results, and more directly comparable to our vibrationless, zerofrequency results. They obtain γxxxx ) 576 au and γzzzz ) 683
R j ) 1/3(Rxx + Ryy + Rzz) and the polarizability anisotropy |∆R|
|∆R|2 )
3 tr R2 - (tr R)2 2
where tr R is the trace of R. In this admittedly restricted set of molecules, the LDA consistently overestimates R j and underestimates |∆R|. The polarizability anisotropy of H2O is surprising. Guan et al. find |∆R| ) 0.43 with their best basis,7 Sim, Salahub, and Chin 0.328 with theirs,6 and McDowell et al. 0.237,23 while our numerical calculations indicate that the polarizability is very nearly isotropic (|∆R| ) 0.02). We investigated the possibility of numerical error by reproducing the calculations with a fine mesh and several different field strengths, but our most exacting calculations increase the anisotropy a negligible amount, to 0.03
TABLE 4: Quadrupole Moments (in au) with Respect to the Center of Mass numericala
a
molecule
Qxx, Qyy
Qzz
Θzz
H2 N2 O2 CO HF NH3 H 2O
-1.614 -7.704 -7.563 -7.729 -4.443 -4.810 -3.378, -5.884
-1.177 -8.842 -7.919 -9.208 -2.736 -7.061 -4.746
0.437 -1.137 -0.356 -1.478 1.707 -2.251 -0.111, 1.937d
Present work. b Reference 13, ANO basis. c Collected in ref 20. d Θxx.
DeMonb Θzz -1.129 -1.471 1.676 -2.201 n/a, 1.913d
exptc Θzz -1.09 0.25 -1.44 1.75 -2.45 -0.10, -1.96d
16108 J. Phys. Chem., Vol. 100, No. 40, 1996
Dickson and Becke
TABLE 5: Dipole Polarizabilities (in au) molecule H2
N2
O2
CO
HF
H2O
NH3
CH4
component
numericala
Rxx Rzz R j |∆R| Rxx Rzz R j |∆R| Rxx Rzz R j |∆R| Rxx Rzz R j |∆R| Rxx Rzz R |∆R| Rxx Ryy Rzz R j |∆R| Rxx Rzz R j |∆R| Rxx
5.28 7.18 5.91 1.9 10.73 15.4(3) 12.3 4.6(3) 8.67 15.47 10.93 6.8 12.61 15.88 13.70 3.3 5.92 6.85 6.23 0.93 10.59 10.61 10.59 10.60 0.02 14.59 17.45(5) 15.54(2) 2.9 17.69
DMolb
expt
5.20 7.11 5.84 1.91 9.98 15.13 11.70 5.15
5.43c 2.0d 11.74c 4.7e 10.78f 7.43f
12.06 15.77 13.30 3.71 5.77 6.94(2) 6.16 1.17 10.58 10.12(1) 10.51(1) 10.40 0.43 14.44 16.76 15.21 2.32 17.50(2)
13.08c 3.6f 5.60c 1.3g
9.64c 0.66g 14.56c 1.94g 17.27c
Present work. b DNP+ basis, ref 7. c Static, collected in ref 21. Static, from molecular beam studies.28 e Static, by extrapolation.29f At 632.8 nm, from ref 30. g At various frequencies, collected in ref 22. a
d
TABLE 6: First Dipole Hyperpolarizabilities (in au) molecule
component
numericala
DMolb
CO
βzxx βzzz β h βzxx βzzz β h βzxx βzyy βzzz β h βxyy βzxx βzzz β h βxyz
8.6 33.7 30.5 -2.3 -10.8 -9.2 -11.9 -9.2 -20.1 -24.8 8. -14.(3) -64.(4) -55.(6) -7.5
8.84 34.90(8) 31.5 -0.93(9) -11.1(5) -7.8 -15.2 -4.26(8) -13.2 -19.6 11.98(8) -16.89(9) -28.2(3) -37.2 -13.9
HF H2O
NH3
CH4 a
exptc
30 (3) -11
-22
-48
Present work. b Reference 7, DNP+ basis. c References 24 and 25.
TABLE 7: LDA Second Dipole Polarizabilities (in au)a molecule H2 N2 O2 CO HF H2O CH4 a
component
value
γxxxx γzzzz γxxxx γxxxx γzzzz γxxxx γzzzz γxxxx γzzzz γxxxx γyyyy γzzzz γzzzz
1100 1300 1100 ( 200 1100 1700 2100 2700 1100 610 1400 5100 2600 3700
Uncertainties are approximately (10%, except as noted.
au. Our results, 1100 and 1300 au, respectively, are clearly too large by a factor of 2, which must be attributed to the local spin-density approximation.
4. Conclusions We confirm the conclusions of Guan et al.7 with regard to basis set requirements for electrical properties, namely, that a good-quality double-zeta basis set augmented with field-induced polarization functions gives well-converged dipole moments µ and polarizabilities R. The local spin-density approximation yields qualitatively correct predictions for these properties at a computational cost comparable to or smaller than the HartreeFock method. Higher polarizabilities β and γ, though, clearly require larger basis sets. We have therefore calculated reference values for these properties in the local spin-density approximation. We observe surprising agreement with experiment for β using only the simple LDA, which must be considered fortuitous until further investigations into vibrational effects are carried out. We have included O2 in our data to illustrate the ease with which density-functional theory treats this important molecule. Traditional ab initio calculations on such open-shell molecules require multiconfiguration techniques which are both complicated and expensive for such derivative properties as hyperpolarizabilities. With density-functional theory, both the cost and accuracy of the treatment are comparable to that of a closedshell molecule. Acknowledgment. This work is supported by the Natural Sciences and Engineering Research Council (NSERC) of Canada. References and Notes (1) Dykstra, C. E.; Liu, S.-Y.; Malik, D. J. AdV. Chem. Phys. 1989, 75, 37. (2) Bishop, D. M. AdV. Quantum Chem. 1994, 25, 1. (3) Luo, Y.; A° gren, H.; Jørgensen, P.; Mikkelsen, K. V. AdV. Quantum Chem. 1995, 26, 165. (4) Moullet, I.; Martins, J. L. J. Chem. Phys. 1990, 92, 527. (5) Jasien, P. G.; Fitzgerald, G. J. Chem. Phys. 1990, 93, 2554. (6) Sim, F.; Salahub, D. R.; Chin, S. Int. J. Quantum Chem. 1992, 43, 463. (7) Guan, J.; Duffy, P.; Carter, J. T.; Chong, D. P.; Casida, K. C.; Casida, M. E.; Wrinn, M. J. Chem. Phys. 1993, 98, 4753. (8) Colwell, S. M.; Murray, C. W.; Handy, N. C.; Amos, R. D. Chem. Phys. Lett. 1993, 210, 261. (9) Matsuzawa, N.; Dixon, D. A. J. Phys. Chem. 1994, 98, 2545. (10) Chong, D. P. Chem. Phys. Lett. 1994, 217, 539. (11) Lee, A. M.; Colwell, S. M. J. Chem. Phys. 1994, 101, 9704. (12) Barone, V. Chem. Phys. Lett. 1994, 226, 392. (13) Duffy, P.; Chong, D. P.; Dupuis, M. J. Chem. Phys. 1995, 102, 3312. (14) Becke, A. D.; Dickson, R. M. J. Chem. Phys. 1990, 92, 3610. (15) Becke, A. D. Int. J. Quantum Chem. Quantum Chem. Symp. 1989, 23, 599. (16) Gross, E. K. U.; Kohn, W. AdV. Quantum Chem. 1990, 21, 255. (17) Buckingham, A. D. AdV. Chem. Phys. 1967, 12, 107. (18) Perdew, J. P.; Wang, Y. Phys. ReV. B 1992, 45, 13244. (19) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in Fortran, 2nd ed.; Cambridge University Press: Cambridge, 1992; Chapter 15 and references therein. Earlier editions of this work may contain bugs in the singular-value decomposition code; we used instead the corresponding LINPACK routine. (20) Bu¨ndgen, P.; Grein, F.; Thakkar, A. J. J. Mol. Struct. (THEOCHEM) 1995, 334, 7. (21) Spackman, M. A. J. Chem. Phys. 1991, 94, 1288. (22) Spackman, M. A. J. Chem. Phys. 1989, 93, 7594. (23) McDowell, S. A. C.; Amos, R. D.; Handy, N. C. Chem. Phys. Lett. 1995, 235, 1. (24) Ward, J. F.; Miller, C. K. Phys. ReV. A 1979, 19, 826. (25) Dudley II, J. W.; Ward, J. F. J. Chem. Phys. 1985, 82, 4673. (26) Adamowicz, L.; Bartlett, R. J. J. Chem. Phys. 1986, 84, 4988. (27) Bishop, D. M.; Pipin, J.; Cybulski, S. M. Phys. ReV. A 1991, 43, 4845. (28) MacAdam, K. B.; Ramsey, N. F. Phys. ReV. A 1972, 6, 898. (29) Alms, G. R.; Burnham, A. K.; Flygare, W. H. J. Chem. Phys. 1975, 64, 3321. (30) Bridge, N. J.; Buckingham, A. D. Proc. R. Soc. London, Ser. A 1966, 295, 334.
JP9605966
