7120
J . Phys. Chem. 1989, 93, 7 120-7 130
Ab I nftio Calculations of Polarizablltties and Second Hyperpdarizabilities of Organic Molecules with Extended ?r-Electron Conjugation Pratibha Chopra, Louis Carlacci, Harry F. King, and Paras N. Prasad* Department of Chemistry, State University of New York at Buffalo, Buffalo, New York 14214 (Received: April 5, 1989)
Static polarizability and second hyperpolarizability tensors are computed for a series of polyenes, polyynes, and cumulenes by ab initio SCF theory. Numerically stable finite field (FF) calculations can be achieved by using polynomial fits of either energy or induced dipole moment as a function of field strength. The nonlinear expansion coefficientsfrom these fits correspond to the microscopic nonlinear optical property. Our results from fully coupled (FF) ab initio calculations for polarizability are in good agreement with those derived from uncoupled (sum-over-states) ab initio methods. The hyperpolarizabilities do not compare as well. A qualitative description of the chain length dependence of polarizability and hyperpolarizability for moderately long chains is discussed in terms of an empirical function. Diffuse orbital basis functions are required for qualitatively correct hyperpolarizabilitiesof small conjugated n systems or for that matter any small molecule. For example, the average second hyperpolarizability, 7,of ethylene is computed to be -13, 1.7, and 726 au with STO-3G, 6-31G, and augmented 3-21G basis sets, respectively. The value computed with the augmented basis set agrees, within a factor of 2, with the experimental value of 1500 au. The valence set for the carbon atoms is augmented with diffuse s, p, and two Cartesian d sets, or subsets of these. The inclusion of diffuse polarization functions drastically alters the computed second hyperpolarizabilities. The results are nearly insensitive to the choice of valence set but are highly dependent on the basis set quality. We also describe the use of a corresponding orbital analysis to aid in the interpretation of ab initio results obtained by either FF or analytic derivative methods. The computed polarizability and hyperpolarizabilityand the n and u components indicate that the contribution due to the x orbitals is much more significant than the u orbitals. Both contributions change sign in going from the ethylenic and acetylenic chains to the cumulenic systems. The polarization of n-electron density by the field is illustrated by contour surfaces of derivative n-electron-density functions. Contour maps of the first derivative of charge density with respect to the field (a)for an acetylenic chain are nearly periodic, corresponding to localized polarization, whereas the third derivative density (y) corresponds to longer range charge shifts.
I. Introduction The prospect of using nonlinear optical effects for optical switching and logic operations for optical computing and signal processing has attracted considerable attention.' Organic materials show promise for these applications because of their relatively large nonresonant nonlinear optical properties and subpicosecond response time^.^,^ As a result, interest in organic nonlinear optics has grown tremendously in recent years. Proposed schemes of optical switching and logic operations utilize third-order nonlinear optical effects which are characterized by microscopic second hyperpolarizabilities y and corresponding bulk susceptibilities x ( ~ ) . The *-conjugated polymeric structures have shown the largest values of nonresonant third-order nonlinearity with the largest tensor component being along the chain direction. This large y has been shown theoretically and experimentally to reflect the effective x-electron conjugation. For weak intermolecular coupling, one can use an oriented gas model and relate the microscopic nonlinearity y to the bulk nonlinearity x ( ~using ) the Lorentz approximation for the local field c ~ r r e c t i o n . ~In such a case X(3)(-W4;W~,WzW3)
= ~ ( w I F(wz) ) F ( w ~F(w4)C(yn(6,4)) ) (1) n
where F(w,)represents the local field at frequency w,. The summation n runs over all molecular sites to give an orientationally averaged ( y ). The two limits are (i) all polymer chains aligned in the same direction in which case x n ( y n ( 0 , 4 ) )= Ny,,,, and = (ii) a completely amorphous polymer for which xn(-yn(6,4)) '/5Nyr,rrr assuming that y,,,, is the dominant tensor element. In either case, a large value of y produces a large value for xt3). Therefore, third-order effects are determined primarily by the microscopic nonlinearity provided local field effects are described (1) Clymer, 8.; Coolins, S. A,, Jr. Opt. Eng. 1985, 24, 74. (2) Williams, D. J., Ed. Nonlinear Optical Properties of Organic and Polymer Materials; American Chemical Society: Washington, DC, 1983. (3) Prasad, P. N., Ulrich, D. R., Eds. Nonlinear Optical and Electroactiue Polymer; Plenum Press: New York, 1988. (4) Shen, Y. R. The Principles of Nonlinear Optics; Wiley: New York, 1984.
0022-36S4/89/2093-7120$01 S O / O
by the oriented gas model. This appears to be a reasonable model for conjugated polymeric systems for which bulk polarization is expected to be weak. Interest centers, therefore, on enhancing the molecular hyperpolarizability y,and it becomes of paramount importance to understand the relationship between molecular structure and hyperpolarizability in order to predict structures with large third-order response. The sign of the nonlinearity is also important because it determines whether self-focusing or self-defocusing will occur in a self-action e ~ p e r i m e n t . ~ Although several theoretical studies of linear and nonlinear susceptibilities of conjugated organic molecules have been reported in recent years,5-" results of rigorous ab initio quantum mechanical calculations for such molecules are only now beginning to appear?*'0 It is unlikely, however, that they will completely supplant the various semiempirical methods in current use. The two approaches have rather different limitations, and so are complementary. Ab initio computations are subject to systematic errors due to the use of finite orbital basis sets and finite many-body basis sets, Le., truncation of the configuration interaction expansion. Their strength lies in the fact that they build on a sound theoretical foundation using well-defined mathematical approximations that can be improved by straightforward, albeit computationally expensive, extensions of the method. Thus they provide a reliable testing ground for empirical and semiempirical methods. (5) McIntyre, E. F.; Hameka, H. F. J. Chem. Phys. 1978, 68, 3481. Zamani-Khamiri, 0.; Hameka, H. F. J. Chem. Phys. 1980, 72, 5903; Ibid. 1980, 73, 5693. (6) Zyss, J. J. Chem. Phys. 1979, 70, 3333. (7) Lalama, S. J.; Garito, A. F. Phys. Reu. A . 1979, 20, 1179. Garito, A. F.; Teng, C. C.; Wong, K. Y.; Zamani-Khamiri, 0. Mol. Cryst. Liq. Cryst. 1984, 106, 219. (8) Papadopoulos, M. G.; Waite, J.; Nicolaides, C. A. J. Chem. Phys. 1982, 77, 2527. Waite, J.; Papadopoulos, M. G. J. Chem. Phys. 1985, 82, 1427. (9) Andre, J. M.; Barbier, C.; Bcdart, V. P.; Delhalle, J. Nonlinear Optical Properties of Organic Molecules and Crystals; Chemla, D. S., Zyss, J., Eds.; Academic Press: London, 1987; Vol. 2, p 137. (10) Hurst, G. J. B.; Dupuis, M.; Clementi, E. J. Chem. Phys. 1988, 89, 385. (11) deMelo, C. P.; Silbey, R. Chem. Phys. Left. 1987, 140, 537.
0 1989 American Chemical Society
Hyperpolarizabilities of Organic Molecules The theory of an isolated molecule interacting with a general electric field has been formulated by Buckingham and others.I2 Here we confine our attention to the special case of a uniform static electric field. The quantum mechanical calculation of static polarizability, a,and hyperpolarizability, y, has been carried out exactly for the hydrogen atomI3 and to high accuracy, using variational wave functions with explicit dependence upon interelectron distances, for the helium atomI4 and hydrogen m01ecule.I~ Nelin16 et al. report accurate calculations of a for F, F, and Ne using complete active space multiconfiguration self-consistent field (CASSCF) theory. In a recent comparative study, Bauschlicher and TaylorI7 have computed polarizabilities for F, HF, CH2, and SiH2 using several different ab initio formalisms in current use including CASSCF, single configuration SCF, configuration interaction with all single and double excitations out of a single configuration reference function (SDCI) or out of a multireference function (MRCI), coupled pair functional (CPF), and full configuration interaction (FCI) methods. This systematic study provides an essential benchmark for testing electron correlation effects. Few of these formalisms can presently be applied to the large organic molecules of interest, but that situation can change rapidly as a result of current developments in computer technology and computational methods, notably, in the development of analytic derivative methods.'* Analytic methods, pioneered by Gerratt and Mills,Ig Pulay,Ig Schlegel and Pople,zo and others,21 avoid the infinite sums over excited states, which arise in the spectral expansion of traditional perturbation theory, by using modern computational techniques for solving large linear systemsm For example, static polarizabilities, as well as infrared frequencies and intensities, have been routinely calculated in recent years for molecules of 10-30 atoms by analytically solving the fully coupled ab initio S C F equations using GRADSCF codes on a CRAY comp ~ t e rand , ~very ~ ~recently ~ Hurst, Dupuis, and Clementi'O have developed closely related methods for hyperpolarizabilities and applied them to polyenes as large as C22H24. The development of a b initio theories of frequency-dependent properties is in p r o g r e s ~ . ~Dykstra ~ * ~ ~ et al. have developed a version of fully coupled S C F known as derivative Hartree-Fock (DHF) which has been used to compute high-order hyperpolarizabilities for a number of first-row hydrides.26 Other reported rigorous ab initio computations of second-order hyperpolarizabilities include be-
(12) Buckingham, A. D. Adu. Chem. Phys. 1967, 12, 107. McLean, A. D.; Yoshimine, M. J. Chem. Phys. 1967, 47, 1927. (13) Sewell, G. W. Proc. Cambridge Philos. SOC.1949, 45, 678. (14) Grasso, M. N.; Chung, K. T.; Hurst, R. P. Phys. Reu. 1968,167, 1. Sitz, R.; Yaris, R. J. Chem. Phys. 1968, 49, 3546. (15) Bishop, D. M.; Cheung, L. M. Phys. Reu. A 1979,20, 1310. Hunt, J. L.; Poll, J. D.; Wolniewicz, L. Can. J . Phys. 1984, 62, 1719. (16) Nelin, C.; Roos, B. 0.;Sadlej, A. J.; Siegbahn, P. E. M. J. Chem. Phys. 1982, 77, 3607, and references therein. (17) Bauschlicher, C. W., Jr.; Taylor, P. R. Theor. Chim. Acta 1987, 71, 263. (1 8) Jsrgensen, P., Simons, J., Eds., Geometrical Deriuatiues of Energy Surfaces and Molecular Properties; Reidel: Dordrecht, 1986. (19) Gerratt, J.; Mills, I. M. J . Chem. Phys. 1968,49, 1719, 1730. Pulay, P. Mol. Phys. 1969, 17, 197. (20) Pople, J. A.; Krishnan, R.; Schlegel, H. B.; Binkley, J. S. Int. J . Quanlum Chem. Symp. 1879, 13, 225. (21) Diercksen, G.; McWeeny, R. J. Chem. Phys. 1966, 44, 3554. Komornicki, A.; Ishida, K.; Morokuma, K.; Ditchfield, R.; Conrad, M. Chem. Phys. Lett. 1977,45,595. Dodds, J. L.; McWeeny, R.;Raynes, W. T.; Riley, J. P. Mol. Phys. 1977,33,611. Brooks, B. R.; Laidig, W. D.; Saxe, P.; Handy, N. C.; Schaefer, H. F. Phys. Scr. 1980, 21, 312. (22) GRADSCF is ab initio gradient program system designed and written by A. K. Komornicki at Polyatomics Research Institute and supported on
grants through NASA. (23) King, H. F.; Komornicki, A. K. Geometrical Deriuatiues ofEnergy Surfaces and Molecular Properties; Jsrgensen, P., Simons, J., Eds.; Reidel: Dordrecht, 1986; p 207. (24) Sekino, H.; Bartlett, R. J. J. Chem. Phys. 1986,84, 2726; Ibid. 1986, 85, 976. (25) Olsen, J. Geometrical Deriuatiues of Energy Surfaces and Molecular Properties; Jsrgensen, P., Simons, J., Eds.; Reidel: Dordrecht, 1986; p 157. Reinsch, E. A. J . Chem. Phys. 1985, 83, 5784. (26) (a) Dykstra, C. E.; Jasien. P. G. Chem. Phys. Lert. 1984, 109, 388. (b) Dykstra, C. E. J. Chem. Phys. 1985, 82, 4120. (c) Liu, S. Y.;Dykstra, C. E.; Malik, D. J. Chem. Phys. Lett. 1986, 130, 403.
The Journal of Physical Chemistry, Vol. 93, No. 20, 1989 7121 ryllium, neon and argon atoms,27water,28m e t h a r ~ e ,fluoro~~,~~ me thane^:^ diatomic nitrogen,30hydrogen cyanide,30acetylene,30 and diacetylet~e.~~ The study of hydrogen fluoride by Sekino and BartlettZ4appears to be the currently most careful comparison of theoretical and experimental molecular hyperpolarizabilities and concludes that an unresolved discrepancy of about 30% remains. Early theoretical studies of nonlinear susceptibilities of polyenes by Hameka and co-workers5employed Pariser-Parr-Pople (PPP)" and extended Hiickel Hamiltonians. Zyss used the intermediate neglect of differential overlap (INDO) approximation to study substituted benzenes.6 Garito et al. computed excited states by SDCI in the complete neglect of differential overlap (CNDO) approximation and substituted the results into sum-over-states (SOS) formulas of perturbation theory to compute hyperpolarizabilities for substituted nitroanilines and other m o l e c ~ l e s . ~ Papadopoulos et a1.8 modified the C N D O parameters to fit polarizabilities and applied them using McWeeny's version of the coupled Hartree-Fock equations2' to compute y for polyenes, aromatic amines, and related molecules. Andre et al. carried out ab initio S C F calculations of a for the polyene series ethylene through octatetraeneg and computed y tensors for these molecules by an SOS formula equivalent to the ab initio uncoupled Hartree-Fock formalism.32 Melo and Silbeyl' have recently computed first and second hyperpolarizabilities for neutral and electrically charged polyene chains using an analytical derivative method based on McWeeny's method21 and the PPP a p p r ~ x i m a t i o n . ~ ~ In this paper, three computational aspects are discussed in detail. The first concerns the numerical stability of finite field techn i q u e applied ~ ~ ~ ~to~the ~ computation of linear and nonlinear susceptibilities using power series expansion of either the energy or the dipole moment. This is explored at the S C F level of theory, but the motivation is the prospect of using finite field methods for correlated wave functions not presently amenable to analytical techniques. Next we inquire into orbital basis set requirements with emphasis on the need to include diffuse polarization functions. Finally, we ask what additional information can be extracted from a wave function that might shed light on the relationship between molecular structure and electrical susceptibility. In this spirit we look at the responses of individual orbitals to an applied field and at the spatial regions of a molecule that most contribute to the induced electrical dipole moment. Of particular interest is the partition of the induced moment of a conjugated hydrocarbon into its CT and T components. 11. General Theoretical Considerations
The quantum mechanical electronic Hamiltonian operator for an isolated molecule in a uniform static electric field F is
H(F) =
+ CF-r, - CZ,F.R, I c
where CT and I label electrons and nuclei, respectively, 2,is the atomic number of the Ith nucleus, and Ho is the field-free Hamiltonian. The molecular energy and dipole moment are the appropriate expectation values evaluated for the many-electron wave functions, \k(F), determined in the presence of the field. (3) (27) Sitter, R. E., Jr.; Hurst, R. P. Phys. Rev.A 1972, 5, 5. (28) Lazzeretti, P.; Zansi, R. Chem. Phys. Lett. 1976,39,323. Purvis, G. D.;Bartlett, R. J. Phys. Reu. A 1981, 23, 1594. (29) Harrison, R. J.; Fitzgerald, G. B.; Laidig, W. D.; Bartlett, R. J. Chem. Phys. Lett. 1986, 124, 291. (30) Jameson, C. J.; Fowler, P. W. Chem. Phys. 1986.85, 3423. (31) Pariser, R.; Parr, R. G. J. Chem. Phys. 1953, 21, 767. (32) Caves, T. C.; Karplus, M. J. Chem. Phys. 1%9,52,3649. Liebmann, S. P.; Moskowitz, J. W. J . Chem. Phys. 1971, 54, 3622. Nakatsuji, H.; Musher, J. I. J. Chem. Phys. 1974, 61, 3737. (33) Cohen, H. D.; Roothaan, C. C. J. J . Chem. Phys. 1965, 43, S34. Hush, N. S.; Williams, M. L. Chem. Phys. Lett. 1970, 5 , 507. (34) Gready, J. E.; Bacskay, G. B.; Hush, N. S. Chem. Phys. 1977, 22, 141.
7122 The Journal of Physical Chemistry, Vol. 93, No. 20, 1989 r(F) = (Q(F)I(CZ& I
- CrJIQ(F))
(4)
0
The essence of a fully coupled method is that an approximate wave function, Q(F), is employed which is uniformly valid over the range of field strengths of interest, in contrast to methods in which the zero-field wave function plays a unique role. The advantages of fully coupled methods for susceptibility calculations have long been r e ~ o g n i z e d . ~In~ the fully coupled S C F formalism, Q is approximated by a single Slater determinant whose orbitals are reoptimized for each F. In fuly coupled SDCI the orbitals of the reference determinant and the coefficients in the configuration interaction expansion are all reoptimized for each F. Aside from matters of computational efficiency and numerical stability, it does not matter whether the response of the orbitals to the changing electric field is determined by explicit recomputation of Q for each field strength, the finite field method (FF),33,34or by analytical solution of the response equations, Le., the coupled perturbed Hartree-Fock (CPHF) e q ~ a t i o n s . ' ~ - ~ ' ~ ~ ~ Differentiation of (3) gives
Chopra et al. turbation theory in the spectral expansion formalism. This leads to the familiar SOS formulas which offer certain advantages for frequency-dependent polarizabilities and provide a more direct relationship between polarizability and spectroscopic properties such as HOMO-LUMO band gaps. On the other hand, an adequate description of these spectroscopic states can be a much more difficult theoretical problem than that of calculating ground-state properties. It is known, for example, that extensive configuration interaction is essential for a proper description of low-lying excited singlet states of a polyene chain.37 In contrast, the wave function in our study need only provide an adequate description of the ground electronic state or, to be more precise, an adequate description of how a molecule in its ground state responds to an applied static external field. Some SOS and fully coupled S C F results are compared in section IV below. 111. Method of Approach For calculations reported in this paper we use the closed-shell, single configuration, S C F wave function, namely, an antisymmetrized product of n doubly occupied molecular orbitals @
= A[ha(l) hP(2) h 4 3 ) hP(4)
... ~ 3 ( 2 n ) l
(9)
where 2n is the number of electrons. The spatial orbitals, f$k, are expanded as linear combinations of m contracted Gaussian basis functions, According to the Hellmann-Feynman t h e ~ r e m , ~the ' first and third terms on the right-hand side of ( 5 ) are zero if \k is the true wave function. The theorem is satisfied by various approximate wave functions including SCF, CASSCF, and FCI wave functions computed by using a finite orbital basis. Fully coupled methods tend to satisfy the theorem. The dipole moment is the negative gradient of energy with respect to field.
The Hellmann-Feynman theorem asserts that (3), (4),and ( 6 ) are compatible and that the following two series expansions have common coefficients.
m
4k =
a=
I
xucuk(F)
(10)
with coefficients, Cuk(F),chosen to minimze the energy (3) in the presence of a field F. In matrix notation (10) becomes
4 = XC(F)
(1 1)
where 4 and x and are row vectors of molecular orbitals and atomic orbital basis functions, respectively. In the finite field method, either energy or dipole moment is computed for several different field strengths and fitted by a polynomial according to (7) and (8). In the analytic method, the field dependence of the MO coefficient matrix, C(F), is introduced through an antisymmetric matrix, X(F), according to3*
where Co represents zerdield orbitals and X is block off-diagonal. X = I -M' 0
The Hellmann-Feynman theorem is not satisfied, for example, in Mdler-Plesset theory, or by an SDCI calculation using fixed, zero-field SCF orbitals. The popular opinion16.Mis that (7) is more valid than (8) when using an approximate \k that does not satisfy the theorem, but there is little computational evidence in support of this opinion. Tensor elements in (7) and (8) are symmetric with respect to permutation of indices. The terms in (7) with odd powers of F are zero for centrosymmetric molecules and also for various other molecular symmetries.I2 We confine our attention in this paper to such symmetric molecules. In comparing results reported in the literature, it is well to keep in mind that in a widely used alternative second hyperpolarizability tensor elements Yjjk( are defined to be fourth partial derivatives of the energy with respect to field. Such yVklhave values 6 times greater than those in (7) and (8). Many computations of polarizability and hyperpolarizability reported in the literature, particularly those for molecules with extended T c o n j u g a t i ~ n , ~employ -~ Rayleigh-Schrdinger per-
We refer to M as the mixing matrix. Element Mij indicates the extent to which the ith occupied zero-field orbital is mixed with thejth virtual by the given field. Derivatives of M with respect to field are obtained by solution of the CPHF equation^.'^.^^ The block off-diagonal form of X provides a complete set of linearly independent variational parameters for an S C F ~ a l c u l a t i o nbut ~~ does not preserve canonical S C F orbital form. In other words, (12) properly describes the field dependence of Q(F) and thus of the energy and charge density, but the orbitals described by (12) are not eigenfunctions of the one-electron Hartree-Fock operator. The concept of corresponding orbitals39provides a useful description of how S C F orbitals respond to the application of an external electric field.40 The S C F wave function, (9), is invariant (37) Hosteny, R. P.; Dunning, T. H., Jr.; Gilman, R. R.;Pipano, A. J . Chem. Phys. 1975,62,4764. Schulten, K.; Okmine, I.; Karplus, M. J. Chem. Phys. 1976,64,4422. Ohmine, I.; Karplus, M.; Schulten, K. J . Chem. Phys. 1978,68, 2298. 51. Hirschfelder, J. 0.;Byers Brown, W.; Epstein, S. T. Adu. Quantum Chem. 1964, 1 , 256. (38) King, H. F.; Champ, R. N.; McIver, J. W., Jr. Chem. Phys. 1983, Rn 1 1~. - 1 1_ . __, (39) Amos, A. T.; Hall, G. G. Proc. R. Soc. London, Ser. A 1961, 263, 483. King, H. F.; Stanton, R.E.; Kim, H.; Wyatt, R.E.; Parr, R.G. J. Chem. Phys. 1967, 47, 1936. (40) Martin, R. L.; Davidson, E. R.; Eggers, D . F. Chem. Phys. 1979, 38, ~
(35) Feynman, R.Phys. Rev. 1939,56,340. Epstein, S. T. Am. J . Phys. 1954, 22, 613. Stanton, R. E. J . Chem. Phys. 1962, 36, 1298. (36) Nerbrant, P. 0.;Rws, B. 0.; Sadlej, A. J. Int. J . Quantum Chem. 1979, 15, 135. Diercksen, G . H. F.; Roos, B. 0.;Sadlej, A. J. Chem. Phys. 1981, 59, 29.
MI 0
341.
The Journal of Physical Chemistry, Vol. 93, No. 20, 1989 7123
Hyperpolarizabilities of Organic Molecules under a unitary transformation that does not mix occupied orbitals with virtuals. Such a transformation can be applied to the finite field orbitals, 4(F), to produce an equivalent set, &F). Similarly, the zero-field orbitals can be transformed into an equivalent set, $ O , such that occupied orbital, $"k, mixes only with its corresponding virtual, $'k+n. The effect of the applied field is then described by n two-by-two rotation^.^^.^' k 5n $k(F) = COS 6k$ko Sin &&+KO (14)
+
When analytic methods are employed, the same result is achieved by diagonalizing the mixing matrix
d = U+MV
(16)
where d is a real, nonnegative, diagonal matrix39with elements dkk = e k , and U and V are the appropriate unitary transformations that bring the occupied and virtual zero-field orbitals, respectively, into correspondence.
C" =C'I
u o o v
1
(17)
In the weak field limit the rotation angle, e k , is proportional to field strength and the corresponding occupied-virtual pairs, $ko and &+ko, are stationary with respect to small variations of field strength. Let $ko in (1 4) and (1 5) be an occupied corresponding orbital for a field applied along, say, the z axis. Then substituting (14) and (1 5) into (4) and (6) and extracting the term linear in F, gives a polarizability tensor element as a sum of contributions from individual corresponding orbital pairs, e.g. k= I
The dipole moment matrix elements, ( $k+,,o(il$kO ), are negative by construction. The most polarizable orbitals are those with large rotation angles and large dipole moment matrix elements. As discussed briefly in section IVD below, nonlinear susceptibility arises from the field dependence of both factors on the right-hand side of (1 8). Linear and nonlinear susceptibilities can also be interpreted in terms of the spatial regions which contribute to the field dependence of the electron charge density. The charge density function can be expanded in powers of field components in much the same way as energy and dipole moment have been expressed in (7) and (8). To simplify the notation, consider a fixed field direction and express the electron charge density function, p(r,F), as a Taylor series in field strength, F = 1FI.
The dipole moment expansion (8) can be recovered by integration of (19), e.g.
Polarizability and hyperpolarizability tensor elements can be expressed as integrals over the appropriate density derivatives. For example, if the field is directed along the z axis, then cy,,
= j z p ( ' ) ( r ) dr3
and
We have computed the charge density derivative functions, p(m), by finite field method. The one-particle density matrix, D, is computed at each of several field strengths, and the individual matrix elements are fitted by polynomials. (41) Camp, R. N.; King, H.F. J . Chem. Phys. 1981, 75,268. See eq 22 and 23 and appendices.
Density derivatives are then computed by using the coefficients of the polynomial fits. d m ) ( r )=
C@:' X , ( d
Xp(4
(24)
U,P
It is also possible to partition the density derivatives into contributions from individual orbitals or groups of orbitals, e.g., u and x contributions. Contour maps of derivative densities34are displayed in section IV.
IV. Results and Discussion A . Finite Field Method. The finite field method provides a straightforward computational technique applicable to almost any quantum chemical formalism, but its applicability to the computation of hyperpolarizabilities for moderately large conjugated x systems has been que~tioned.~ Indeed, the method does require care to avoid numerical differencing errors and is less efficient than analytical methods where they are available, but we have found it to be a viable method for all of the several conjugated molecules to which it has been applied. Success is attributed to three factors: computation of highly precise energies and density matrices, careful selection of field strengths, and polynomial fitting to several data points. The accuracy to which a susceptibility tensor is computed is determined by the numerical stability of the finite field method, as well as by the quality of the orbital and many-body bases. We are presently concerned only with the issue of numerical stability. Table I gives detailed results of a finite field calculation of a, and -y,z using the induced moment expansion and a 3-21G orbital basis4* for the S C F wave function computed with a suitably modified HONDO code.43 The molecule is octatetrayne, CsH2. The second column of Table I gives the strength, in atomic units, of an electric field applied along the molecular z axis. The induced dipole moment is reported to eight decimal places in column three. Since dipole moment integrals are of about unit magnitude and since the convergence tolerance on the density matrix has been set to lo-' in the S C F calculations, we expect the dipole moment values to be precise to at least six decimal places. A least-squares linear fit of the 12 data points gives a standard deviation of 3 X more than 3 orders of magnitude greater than the anticipated computational error. Deviations from the linear fit for each of the 12 points are reported in column 4 of Table I. It is clear that these are systematic deviations from a linear function and that points 5-12 contain at least four significant figures of information about the nonlinearity. A two-parameter fit with linear and cubic terms reduces the standard deviation to but even these small deviations from a cubic polynomial, reported in column 5 of the Table, do not appear to be random fluctuations. Adding an F5 term to the fitting polynomial further reduces the standard deviation by 2 orders of magnitude and provides our best estimates of aZzand ~ r r r , . Deviations from three- and four-parameter fits, reported in the last two columns of Table I, represent the noise in the computed dipole moment values. To achieve this precision, one, must avoid nearly linearly seven decimal places in ~ ( 9 dependent orbital basis sets. The smallest eigenvalue of the overlap in this example. The nearly redundant function matrix is 2 X was retained, but very nearly the same a and 7 are obtained without it, provided that one either accepts it or rejects it consistently for all fields. Addition of diffuse functions leads to even more severe basis set redundancies, in which case the nearly redundant linear combinations must be rejected. The aZzand Y~~~~values obtained from the three-parameter fit of the induced moment, a,, = 206.5176 and yrzZz= 19388, are nearly identical with those from the four-parameter fit, and are in excellent agreement with the results cy,, = 206.5176 and Y~~~~ = 19 387 given by a four-term polynomial fit of the SCF energies. (42) Binkley, J. S.,Jr.; Pople, J. A,; Hehre, W. J. J . Am. Chem. Soc. 1980, 102, 939.
(43) Dupuis, M.; Rys, J.; King, H. F. J . Chem. Phys. 1976, 65, 111.
7124
The Journal of Physical Chemistry, Vol. 93, No. 20, 1989
Chopra et al.
TABLE I: Polynomial Fits of the Dipole Moment Induced in %Ha by an Electric Field Applied along the Molecular z Axis'
deviation: r(F) polynomiald point
P
1 2 3 4 5 6 7 8 9
0.0001 0.0003 0.0005 0.0007 0.0010 0.0020 0.0025 0.0040 0.0050 0.0070 0.0090 0.0100
10 11 12
0.020 651 52 0.061 95581 0.103261 22 0.144 568 97 0.206 537 09 0.413 19045 0.5 16 597 3 1 0.827 31443 1.035021 14 1.452 324 5 1 1.872971 13 2.084 867 94
1
1+3
1+3+5
1+3+5+7
-0.000 138 32 -0.000 41 3 70 -0.000 687 97 -0.000 959 90 -0.001 361 29 -0.002 606 3 1 -0,003 148 64 -0.004 279 08 -0,004 470 75 -0.002 964 14 0.001 886 32 0.005 884 15
o.Oo0 000 47 0.000002 18 0.000003 60 0.000 004 98 0.000 007 03 0.000 01 1 86 0.000 012 99 0.000 009 40 0.000 001 25 -0.000 020 80 -0.00001646 0.000018 07
-0.000 000 26
-0.000 000 26 0.000 OOO 000 -0.000 00001 -0.000000 02 0.00000008 0.000 000 00
o.oooOOOo00 -0.000 00002 -0.Ooo 000 02 0.000 OOO 07 -0.00000002 0.000 OOO 00 -0.000 00003 0.00000002 0.00000003 -0.000 000 04 0.000 000 02
0.000 000 01 -0.000 00003 0.000 000 00 0.00000001 -0.00000001 0.000 000 00
coefficient of
no. of fitting param
F
1 2 3 4
207.898 38c 206.510 3 l C 206.5 17 63' 206.5 17 62c
P
Fs
F3
19 747f 19 388/ 19 3 8 q
3.04 X lo6 3.02 X lo6
std dev 3.4 x 1.2 x 9.3 x 9.7 x
1.4 X lo8
10-3 10-5
10-8 10-8
'Computed by SCF theory using a 3-21G basis. Nuclear coordinates: f1.3214, *3.5567,16.2085, *8.4308, f10.447 bohr. 1 bohr = 0.529 167 1 au = 1.7152 X IO' esu (statvolt/cm). *Ab initio dipole moment in atomic units. 1 au = 2.5416 X esu (statcoulombmn). dThe fitting polynomial contains one to four terms in odd powers of F. cPolarizability.CY,, in atomic units. 1 au = 1.48176 X esu (cm3). fHyperpolarizability, yzzrz,in atomic units. Note that many authors adopt a convention in which all elements of y are 6 times larger than those in this paper. 1 au = 5.0372 X lo* esu ("/erg).
A. *Field strength in atomic units.
TABLE 11: Comparison of Polarizability and Second Hyperpolarizability Tensors Computed for Ethylene by Finite Field and Sum-Over-States SCF Methods STO-3G basis 6-3 1G basis this work'
tensor elem
e=
00
FF
ffxx
10.83
ffX I
0.00
ffYY
ffzz
a YXJXX
Yxxxr YXJYY
Yxxzz YXYYZ YXZZZ
YYYYY YYYZZ YZZZZ
Y
2.51 19.42 10.92 -4.10 0.00 3.45 -4.31 0.00 0.00 -1.29 12.34 -82.50 -12.99
this work'
ref 9
e = 300 FF
FF
12.97 3.72 2.51 17.27 10.92 -12.31 -10.40 5.67 -15.70 3.85 -23.55 -1.29 10.12 -51.51 -12.99
13.20
sos 11.11
e=
00
FF
19.50
0.00 2.53 17.48 1 1.07
1.93 15.12 9.39 4.03
7.18 32.84 19.84 61.25
0.00 22.01 -83.85 -2.72 43.92 -21 1.86 -49.28
1.52 66.27 0.00 0.00 2.68 14.67 -220.20 1.73
e = 300
ref 9
FF
sos
22.84 5.78 7.18 29.51 19.84 95.25 -0.69 4.80 -38.09 5.69 -121.19 2.68 11.38 -45.48 1.73
18.43 6.22 24.31 16.34 459.23 46.99 -77.31 27.13 94.53 -23.58 118.24
'SCF 6-31G optimized geometry with bond lengths CH 1.0807 A, CC 1.3200 A, and HCH angle 121.065O. Molecule lies in the xz plane with angle 0 between the CC bond and the z axis. The principal axis frame is 0 = Oo. The molecular orientation in ref 9 is equal to or very close to 0 = 300. Tensor elements are given in atomic units The same value, azz= 206.5 176, is computed analytically using the GRADSCF programzz with the same basis set and molecular geometry. The coefficient of is determined by this calculation to be 3.0 X lo6. Twelve well-chosen field strengths were used for the polynomial fits reported in Table I. Not all finite field computations reported in this paper have been carried out so carefully, but in no other sense is this example atypical. Obviously, high numerical stability can be achieved for moderately large conjugated systems using the finite field method. In general, to determine a and y one must sample the linear range and the beginning of the nonlinear range in p(F). We suggest that one estimate the field strength F'for which the fourth term in (8) is 3 orders of magnitude greater than the noise level, and pick a few points in the range F' to 3F' plus a few more scattered lower field points. Obviously, this requires some experience. Essentially the same numerical precision is achieved for noncentrosymmetric molecules, but with twice the computational effort, by computing p ( F ) + ~(-9. Either the energy expansion or the moment expansion can be used, but the moment expansion offers various advantages. One obtains all three com-
ponents of the induced moment from one SCF calculation, so once yzzzzhas been computed, -yxzzzz and yyzzzcan be obtained (if nonzero) with negligible additional computational effort. For a molecule such as CsHz, the yxxrrtensor element is obtained by applying a field of varying magnitude but fixed direction in the xz plane. Even a small ymz tensor element is computed accurately from the polynomial fit of r,(F)by picking a field direction fairly close to the z axis. Selection of an optimal set of field points for the computation of hyperpolarizabilities requires care and experience, but it has a decisive effect on the quality of the results. B. Comparison of Finite Field and SOS Computations. We computed polarizability and hyperpolarizability tensors for ethylene by the SCF finite field (FF) method for the purpose of comparing them with the sum-over-states (SOS) results of Andre, Bodart, and Delhalle? We use the same basis sets and compute all integrals as did Andre et ai., so that any observed differences in the results can be attributed to differences between the FF and SOS methodologies. Table I1 reports a and 7 tensors computed with minimal (STO-3G) and split-valence (6-3 1G) Gaussian orbital basis sets" for the ethylene molecule in its 6-31G S C F
The Journal of Physical Chemistry, Vol. 93, No. 20, 1989 7125
Hyperpolarizabilities of Organic Molecules
TABLE 111: Polarizability and Second HyperpolarizabilityTensors Computed for trans-Polyenes Using Coupled SCF Theory with Minimal and Split-Valence Orbital Basis Sets tensor elem ax,
a,, aYY
a,, Yxxxx Yxxxr Yxxrz Yxzrr YYYYY Yzzzz
C2H4"vb STO-3G 12.97 3.72 2.51 17.27 -12.3 -10.4 -1 5.7 -23.6 -1.3 -5 1.5
3-21G 21.95 5.43 5.81 28.23 -3.2 -37.8 -28.4
-60.0 -1.6 -116.1
C4H6' STO-3G 23.98 9.86 4.64 48.1 1 -23.8 -34.9 -76.1 -158.1 -2.2 534.8
CaH8' 3-2 1G 38.71 13.27 10.69 72.10 94.4 -61.1 -182.3 -478.2 -1.8 906.0
STO-3G 35.23 17.8 1 6.75 94.72 -40.1 -80.2 -230.0 -523.0 -2.1 4976.0
C8H,0a 3-21G 55.60 23.40 15.47 137.60 112.6 -165.0 -638.0 -1507.0 -2.0 8805.0
STO-3G 46.73 26.54 8.91 150.72 -61.1 -146.6 -550.0 -945.0 -3.5 18152.0
3-21G 72.81 34.81 20.28 217.94 116.5 -324.0 -1 312.0 -3193.0 -2.5 32200.0
'SCF 6-31G optimized geometry. Molecule lies in xz plane with z axis passing through midpoints of CC bonds. The defined z axis does not coincide with a principal polarizability axis. In the case of CsHlothe angle between them is 13'. axy= ayr= yxxxy= yxxyr= yxyyy= Y~~~~= yyyyz = 30° as in Table 11. Y~~~~ = 0.
--
optimized geometry. Here we focus on how the results vary with respect to methodology for a given basis set and return in the next section to discuss the basis set effects themselves. STO-3G results are given in columns 2 and 3, and 6-31G results are given in columns 6 and 7 of Table 11. Tensor elements reported in columns 2 and 6 (0 = Oo) are expressed in the symmetry-adapted frame. Those in columns 3 and 7 (e = 30") were obtained by a tensor transformation corresponding to a 30° rotation of the molecule about its y axis, as in ref 9. The precise molecular geometry is not reported there, but Figure 2 in that paper portrays the molecule lying in the x z plane with its carbon-carbon bond titled approximately 30" away from the z axis as in a trans-polyene chain with 120° bond angles. The average polarizability and average second hyperpolarizability Y =
'/s(Yxxxx
+ Yyyyy + Yrrrr + 2Yxxyy + 2Yxxzr + 2Yyyzr) (26)
are invariants of the tensor transformation.12 The finite field polarizabilities computed by Andre et al.? reproduced in column four of Table 11, are about 1% larger than ours in column three. This can be. attributed to small differences in choice of bond lengths and angles and is hardly significant. Andre et al. noted that the SOS method underestimates STO-3G polarizabilities by 16-19% percent. We observe similar behavior for the 6-31G basis. Far more remarkable are the observed discrepancies between F F and SOS hyperpolarizabilities. The SOS method overestimates y, (26), by a factor of 3.5 in the STO-3G calculation and by nearly 2 orders of magnitude in the 6-31G example. Trends in the values of individual tensor elements are poorly reproduced and one instance of a disagreement in sign is observed, namely, yxxrr= -1 2.3 1 (FF) versus yxxxx= 4.03 (SOS). The disagreement between F F and SOS is much worse for y than for a. The discrepancies displayed in Table I1 reflect well-known differences between fully coupled and uncoupled Hartree-Fock theory.9~27~32-34 That the uncoupled equations were solved by spectral expansion (SOS) is not of fundamental importance; they could have been solved with the same result by an alternative technique, e.g., by minimizing the Hylleraas function.4s The essential distinction is whether the Coulomb and exchange operators are fixed (uncoupled) of field dependent (coupled). In the uncoupled formalism, application of the field perturbs the occupied orbitals and so destroys self-consistency between the eigenfunctions of the Hartree-Fock operator and the orbitals defining that operator. It follows that the uncoupled formalism is less valid for the molecule in a field than it is for the field-free molecule. In coupled Hartree-Fock theory, as the occupied orbitals respond to the applied field, so do the Coulomb and exchange operators. The results in Table I1 suggest that this is a minor (44) Hehre, W. J.; Ditchfield, R.; Pople, J. A. J. Chem. Phys. 1972, 56, 2257. (45) Hylleraas, E. A. Z . Phys. 1930, 65, 209. Hylleraas, E. A. Ado. Quantum Chem. 1964, 1 , 1.
correction for the calculation of second derivatives of energy with respect to field but is of major importance for the calculation of fourth derivatives. C. Orbital Basis Sets. In addition to differences between F F and SOS results discussed above, Table I1 also exhibits remarkably large orbital basis set effects. Note, for example, that the minimal and split-valence bases give appreciably different average polarizabilities, namely, a = 10.9 and 19.8 au, respectively, neither of which agrees closely with the experimental value$6 a = 28.5 au. Hyperpolarizability results are seen to be even more sensitive to choice of orbital basis set with y values -13, 1.7, and 1500 au given by minimal basis, split-valence basis, and respectively. Sitter and discuss basis set requirements for atomic calculations, and show that the F-r operator acting on an atomic orbital with angular momentum 1 generates components with angular momentum I f 1 and I f 2 in the first- and second-order perturbation functions, respectively. In accord with the Wigner 2n 1 rule,s1the first-order wave function is required for a (or p) and the second-order wave function for y (or 6). These basis set requirements are the same for F F or analytic derivative calculations. One concludes that for an atom with s- and pvalence orbitals (e.g., neon) diffuse s, p, and d orbitals are required for the computation of a and additional f orbitals are required for y , in agreement with computational e ~ p e r i e n c e .Similar ~~ considerations apply to a molecule, but y for a typical conjugated molecule is far greater than the sum of contributions from its elements, so purely atomic polarization constitutes only a small fraction of the total effect. Dykstra26refers to the primary effect as interatomic valence polarization which is qualitatively described by the usual valence basis sets. On the other hand, polarization is generally thought to take place mainly in the tails of an orbital where the nuclear electric field is relatively weak. The literature contains a number of reports of systematic attempts to select diffuse Gaussian orbital exponents for computing molecular pol a r i z a b i l i t i e ~and , ~ ~ all arrive at very similar conclusions. We tentatively adopt Chong'ss3 suggested polarization functions,
+
(46) Bridge, N. J.; Buckingham, A. D. Proc. R.Soc. London, Ser. A 1966, 295, 334. (47) Ward, J. F.; Elliott, D. S.J. Chem. Phys. 1978, 69, 5438. (48) Buckingham, A. D.; Bogaard, M. P.; Dunmur, D. A,; Hobbs, C. P.; Orr, B. J. Trans. Faraday SOC.1970, 66, 1548. (49) Bogaard, M. P.; Orr, B. J. Int. Rev.Sci.: Phys. Chem. Ser. Two 1975, 2, 149. (50) The average second hyperpolarizabilityof ethylene in the gas phase, y(694 nm) = 1505 34 au, measured by Ward and Elliott using the technique of static electric field induced second harmonic generation (EFISH) agrees, within estimated experimental error, with third harmonic generation (THG) results (ref 47), but not with Kerr cell results, y(633 nm) = -100 au, reported by Buckingham et al. (ref 48). The experiments measure different combinations of electronic, vibrational, and rotational contributions as is discussed in ref 53. For discussion of these experiments see ref 14 and 49. (51) Hirschfelder, J. 0.; Byers Brown, W.; Epstein, S.T. Ado. Quantum Chem. 1964, 1 , 256. (52) Chong, D. P. Int. J. Quantum Chem. 1983,23,447; J . Chem. Phys. 1983, 78, 5287. (53) Shelton, D. P. Mol. Phys. 1987, 60, 65.
*
7126 The Journal of Physical Chemistry, Vol. 93, No. 20, 1989
to those, 1.66 and 1.54, obtained from similar minimal basis and split-valence ab initio calculationsg by Andre et al. Hurst et al.1° suggest an alternative empirical formula that is valid in the limit of inifitly long chains. Their calculated polarizabilities correspond to powers of n varying from 1.4 to 1.6 for various basis sets and chain lengths. Our values for longitudinal hyperpolarizability, X = 3.27 and 3.14, are somewhat lower than semiempirical results but in reasonable agreement with other a b initio values. Hurst et a1.l0 find that the power of n decreases from 4 to 3 over the range n = 6-1 1 and presumably approaches unity in the infinite chain limit. Let us accept, as a working hypothesis, that (27) provides a valid extrapolation from short to moderately long chains. We have already seen that a minimal basis set calculation hopelessly underestimates y for ethylene. If the main effect of improving the basis set were to change only the A parameter in (27), then it would not be of great importance for long chains, but that is not the case. The B parameter for yzrrrmore than doubles in going from a minimal to a split-valence basis set. Inspection of Table 111 reveals this and various other features such as the rather erratic behavior of yxxxx.The improved radial flexibility provided by the split valence basis is obviously important. Table V shows a marked difference for longitudinal tensor components computed with minimal and split-valence sets but relatively little difference among the various split-valence results themselves, Le., among 3-21G, 6-31G, and 6-31G*. The later basis contains a d set on carbon with Gaussian exponential parameter a = 0.8, appropriate for polarization associated with bond formation. It is not sufficiently diffuse to describe atomic polarization by an external field. The main effect of adding these d functions is probably to slightly contract the valence orbitals and so reduce polarizability. Transverse hyperpolarizability components, not shown in Table V, are somewhat more sensitive to choice of split valence basis than are longitudinal components. Tables VI and VI1 report the results of augmenting 3-21G and 6-31G valence sets with one set of diffuse s and p functions and two sets of diffuse Cartesian d's on carbon atoms. Parameters are given in Table VIII. Addition of these diffuse functions has a profound effect on both CY and 7 , particularly for the latter. As indicated in Table VI, the augmented 3-21G set gives very nearly the same results as does the augmented 6-31G set. We have tested these same diffuse functions with extended valence basis sets for polarizability calculations and again find the results to be nearly insensitive to choice of valence set. Even the augmented STO-3G set is reasonably good for this purpose, although it would be
TABLE I V Parameters in (27) for Selected Polarizability and Hyperpolarizability Tensor Elements in Table 111 tensor elem CY,
aXx
CY,, CY,^
Y~~~~ ~rrrr
parameters"
orbital basis
A
B
8
x
STO-3G 3-21G STO-3G 3-21G STO-3G 3-21G
10.60 10.36 16.82 24.1 1 -51.6 -116.1
10.80 16.24 29.43 35.95 461.5 1018.1
0.77 0.28 0.954 0.778 0.924 1.000
1.03 1.025 1.36 1.44 3.27 3.14
"The A and 6 parameters become redundant in the limit X
-
I.
keeping in mind the results of Sitter and Hurst, namely, that functions suitable for polarizabilities are not necessarily adequate for hyperpolarizabilities. The hyperpolarizability of a polyene chain increases much faster than linearly with chain length. Even a minimal basis set calculation reproduces this qualitative feature.1° Table 111 reports CY and 7 tensors for ethylene, butadiene, hexatriene, and octatetraene computed with STO-3G and 3-21G basis sets. We used the empirical formula
f, = A
Chopra et al.
+ B(n - 6)X
to express chain-length dependencies in this series of C2nH2n+2 molecules. Separate A, B, 6, and X parameters were used to fit each individual tensor element. This functional form is not preserved under a tensor transformation and is physically unrealistic in the limit n a. It is intended to be used for moderately large oligomers in a coordinate frame with z as the longitudinal chain axis (not a uniquely defined direction for short polyene chains). Equation 27 achieves good fits of longitudinal tensor components, but less good for transverse components. Some parameter values are reported in Table IV. We arrived at (27) starting with the simpler relationshipf, = BnX,thought to be valid for long chains, and adding A and 6 parameters to take account of end effects for short chains. A free-electron-gas predicts that CY,, = Bn3 and y,,,, = Bn5 in the limit of large n. Early Huckel = B&28 and and the PPP calculations by Hameka et aL5 gave CY,, y,,,, = Bn5.3. More recent PPP calculations by de Melo and for trans-polyenes. Like (27), these Silbey" give y,,,, = Bn4,25 are physically unrealistic in the limit of an infinite chain; they also correspond to X values higher than our values reported in Table IV. Our X values for a,,,1.36 and 1.44, are much closer
-
TABLE V: Orbital Basis Set Dependence of Selected Polarizability and Hywrmlarizability Tensor Components C2H4a3b basis set STO-3G 3-21G 6-31G 6-31G*f
C~H~"B~
ffzz
Ymr
19.42 3 1.36 32.84 32.30
-82.5 -191.7 -220.2 -194.6
a, 48.1 1 72. IO 74.80 74.20
C2HzLd
CqH2b*'
Yzm
ffzz
Yzm
534.8 906.0 1048.0 992.0
17.59 26.49 27.55 27.29
-20.0 1.9 14.6 16.7
%I
47.92 69.17 71.14 71.25
Ymr
38 1.7 841.1 962.6 882.29
"6-31G geometry as in Tables I1 and 111. bCarbon atoms lie on the z axis. CButadienegeometry as in Table 111. "STO-3G geometry. Nuclear coordinates: f l . 1 0 4 0 , f3.1173 bohr. eSTO-3G geometry. Nuclear coordinates: f1.3297, f3.5493, f5.5646 bohr. fCarbon d orbital exponent is 0.8 bohrm2.
TABLE VI: Polarizability and Hyperpolarizability Tensor Elements for Ethylene and Acetylene Computed Using Orbital Bases with Diffuse Polarization Functions C2H2b orbital basisC
C2H4a
orbital basis' CYXX
ffxx
ffzz %z YXXXX
Yxxxx YZZIZ
Yzm
none 18.82 ( 19.50) 3 1.36 (32.84) 33.96 (61.25) -191.7 (-220.2)
SP 19.30 (19.72) 34.83 (35.66) 141 ( 1 63) 59.8 (40.1)
Pdd 23.49 (23.66) 36.36 (36.53) 243 (255) 525 (516)
spdd 23.55 (23.69) 36.38 (36.53) 272 (271) 528 (516)
none 5.41 (6.55) 26.49 (27.55) -0.72 (3.34) 1.9 (14.6)
SP 10.14 (10.48) 29.10 (29.70) 680 (673) 191.7 (1 92.2)
pdd 17.94 (1 7.99) 30.21 (30.08) 559.8 (605) 515.5 (504.8)
spdd 18.26 (18.22) 30.20 (30.08) 853 (857) 515.1 (502.4)
"Geometry is given in Table 11, 0 = Oo. bGeometry is given in Table V. COrbital basis consists of a 3-21G (or 6-31G) valence basis augmented with the indicated diffuse functions on carbon atoms. See Table VI11 for orbital parameters.
The Journal of Physical Chemistry, Vol. 93, No. 20, 1989 7127
Hyperpolarizabilities of Organic Molecules TABLE VII: Polarizability and Hyperpolarizability Tensor Elements for Butadiene and Butadiyne Computed Using Orbital Bases with Diffuse Polarization Functions C4H6'
C&,6
orbital basis' sp 38.71 41.72 10.69 26.73 arr 72.10 79.82 ( a ) 40.50 49.42 yxxxx 94.4 576 yyyYy -1.8 2488 yzzzr 906 3450 (y) 148.3
none
a,, ayy
pdd 46.21 36.13 83.37 55.24 828 1865 4111 23 12d
9.22 9.22 69.17 8.44 8.44 841.1
16.78 16.78 74.86 772 772 3301
29.20 29.20 76.27 929 929 3105
OGeometry is given in Table 111. bGeometry is given in Table V. 'orbital basis consists of a 3-21G valence basis augmented with the indicated diffuse functions on carbon atoms. See Table VI11 for orbital parameters. dExperimentalvalue of ( 7 ) is 4722 au (ref 47). TABLE VIII: Exponential Gaussian Orbital Parameters for Diffuse Carbon Atom Basis Functions this work HDC' ref 10 JF ref 31 a, aP ad
ar
0.07284b O.0338gb 0.02789 0.03610 0.17020 0.07907
U
orbital basis' none sp pdd
0.0214b
0.0343b
0.0500 0.0500
0.0556e 0.0800 0.1000
OThe 6-31G+PD basis of Hurst, Dupuis, and Clementi.Io Those authors also report other basis sets. b ACartesian d shell contains an s component with effectiveexponent u, = 3ad/7, (95.2% overlap between s and d,, + dyy + dzz). ' A Cartesian f shell contains a p component with effective exponent a p = 5ar/9, (97.7% overlap between p, and f,, + fxyy + fxzz). atrocious for many other purposes such as for geometry optimization. The 3-21G valence set augmented with a diffuse pdd basis gives average polarizabilities a = 27.30 and 55.24 au for ethylene and trans-butadiene, respectively. The calculated value for ethylene is about 5% below the experimental value4 a = 28.48 1 esu). The discrepancy between theory and au (4.220 X experiment is reduced to the point where a more precise comparison requires consideration of vibrational and rotational corr e c t i o n ~ .Our ~ ~ calculated values of y for ethylene and transbutadiene are 726 and 23 12, about half the experimental values esu) and 4722 au (2.38 X of 1505 au (0.758 X esu), respectively, reported by Ward and Elli~tt.~'Our tram-butadiene results, 55.24 and 2752, are in good agreement with those, 53.27 and 2474, obtained by Hurst et a1.I0 with a 6-31G valence set augmented with the diffuse p and d functions given in Table VIII. Our calculations show that the transverse hyperpolarizability component yxxxxis larger than yzzrzin acetylene. This appears to be a common feature of small linear molecules, but certainly not of long chains. As shown in Tables VI and VII, and observed previously by Jameson and Fowler,30 the situation is reversed already in going from acetylene to diacetylene. Basis set requirements for the computation of transverse components tend to be more stringent than those for longitudinal components, so calculation of a and y for a short chain requires a better basis than for a long one. A study of polyalkynes using diffuse polarization functions, reported in part in Tables VI and VII, shows that transverse polarization is nearly a linear function of chain length, axx= 8 + 1 In, as one would expect for purely atomic polarization. Results in Tables VI and VI1 show that diffuse d orbitals contribute significantly to a, in much the same way as they do for first-row atoms. These diffuse d functions are less important for azr A diffuse d set has u, x , and 6 components under cylindrical symmetry. Only the u and x components contribute to azz,and it appears that that role is fulfilled by diffuse s and p orbitals on adjacent atoms. On the other hand, Table VI shows an increase from about 192 to 5 15 au in yzzzl when diffuse d
/.
14.0~
10.0
W
t
a
-
G
6.0 40
LT 20
0
0.02
0.04
0.06
0.08
FIELD (0.u)
Figure 1. Rotation angles for corresponding orbital pairs in ethylene as
a function of field strength. Electric field is applied along the CC bond. Linear response is indicated by solid lines; calculated points are indicated by small circles. Each corresponding orbital is given the labels of the major contributing canonical orbitals. Orbital basis is 3-21G. orbitals are added to the diffuse sp basis for acetylene. The diffuse s and p orbitals have the correct symmetry but, apparently, not sufficient flexibility to describe the second-order orbital response. The diffuse s function can be eliminated with practically no effect on yzzzz, since the s space is well spanned by the Cartesian d sets. It comes as a considerable surprise, therefore, to see the yxxxx tensor component increase from about 560 to 853 when the diffuse s function is added to the diffuse pdd basis. Our results for acetylene, 5 15 and 853 for yzzrzand y-,, respectively, agree well with the values, 486 and 850, obtained by Jameson and Fowler30 using a [5s,4p,ld] valence basis augmented with the diffuse Cartesian d and f sets given in Table VIII. As expected, they observe that the diffuse f set is considerably more important for yxxxxthan for yzzzz. Our basis is deficient in lacking the diffuse f set, but somewhat more complete than theirs in lower angular momentum orbital components. Our 3-21G valence set augmented with diffuse pdd sets for diacetylene yields 3105 and 929 for yzzzz and yxxXx, respectively, as reported in Table VII. Jameson and Fowler obtain 3248 and 923, respectively, with their valence basis augmented with only a diffuse d set. Further study of basis set requirements is needed, but it is already clear that diffuse basis functions are essential for a qualitatively correct description of hyperpolarizabilities of short-chain conjugated molecules. It remains to be seen whether or not they are required to correctly predict qualitative trends in long molecules. D. Corresponding Orbital Analysis. Figure 1 is a plot of rotation angle versus field strength for each of the eight corresponding orbital pairs in ethylene. A 3-21G basis is used and the field is applied along the CC bond axis. The plot exhibits the expected close relationship between orbital energies and extent of orbital mixing. Each plot of B,(F) is labeled to show the principal canonical orbital components. The largest rotation angle x * ) . The next is observed for the HOMO-LUMO pair ( x lower canonical orbital (number 7, a linear combination of CH u bonds) undergoes the next largest rotation, and so on down to k-shell orbitals which are essentially unperturbed by the applied field. The initial slope, (dB/dF),, for the HOMO-LUMO pair is 170' per unit field. This pair contributes 54% of the total molecular polarizability, azz,and makes a large negative contribution to yrzzzwhich is partially canceled by positive contributions from the u orbitals, resulting in the value ~ z z r z= -191.7 au reported in Tables V and VI for this basis set. This is shown graphically in Figure 2. The leading nonlinear term in the power series expansion of B k ( F ) is cubic, since quadratic terms are symmetry forbidden for ethylene. The HOMO-LUMO plot in Figure 1 shows a perceptible negative cubic term, characteristic of basis sets with only a few virtual orbitals. With a minimal basis set, the x orbital and its corresponding x* orbital for ethylene are completely determined by symmetry. In that case the rotation angle, e k , is the only field-dependent variable in (14) and (1 5 ) ,
-
7128 The Journal of Physical Chemistry, Vol. 93, No. 20, 1989
Chopra et al.
TABLE I X
u- and *-Orbital Contributions to Longitudinal Polarizabilities and Hyperpolarizabilities of Linear Acetylenic Chain Molecules Computed Using Minimal and Split-Valence Orbital Basis Sets
C2H2 STO-3G 3-21G 9.56 16.93 26.49 52.78 -50.88 1.90
tensor elem
CsH2
C6H2
C4H2
STO-3G 12.40 35.51 47.92
3-21G 14.59 54.58 69.17
-52.7
17.4
434.4
823.7 841.1
381.7
STO-3G 17.74 73.82 91.55 -206 306 1 2855
3-21G 20.17 110.46 130.62 -185 5785 5600
STO-3G 23.57 121.83 145.41 -450 10 530
3-21G 26.19 180.33 206.52 -630 20018 19 388
10080
2.0
0.0
-2.0
4.0t
"
"
-
"
FIELD (a.u.1
Figure 2. Total induced dipole moment and its u and A components. See caption to Figure 1.
'
'
O
'
;
I
"
"
"
'
.
.
,
k
,e-.
F
- 4 0 4 , . ,bo
Po
.
, ,bo
.
,
/%o
.
! 00
.
, V0
,
,P
,
Figure 4. Contour plot of the third derivative of *-electron charge density with respect to field weighted by distance from internuclear axis, rpu). Molecule is butadiyne, C4H2,with field applied along CC bonds from left -40
/*
0
4.0
.
, /b
.
O
"
, /
.
,
3
bo
'
' 0.
"
! . O0
'
'
I
, I .
0 %
'
"
,
,
Po
"
, bo
I
c *O
"
I. . ! . . I . , . 0. , po O0 VO .? b' eo Figure 3. Contour plot of the first derivative of ?r-electroncharge density with respect to field weighted by distance from internuclear axis, rp('). Molecule is butadiyne, C4H2,with field applied along CC bonds from left to right. Tick marks indicate nuclear positions. Orbital basis is STO-3G (top), 3-21 G (bottom). Broken lines indicate negative density. Contour interval is 0.1 au. -4.0
0 t'
.
,
P'O
.
,
/b
to right. Tick marks indicate nuclear positions. Orbital basis is STO-3G (top), 3-21G (bottom). Contour interval is 1.0 au.
.
o
and polarization saturation occurs; Le., the plot of B k ( F ) approaches a limiting value less than 90'. Similar behavior is exhibited by the 3-21G basis set calculation for ethylene as well as by a number of other small basis set calculations reported above. To obtain a positive hyperpolarizability value requires a low-energy second-order orbital response function lying mainly outside the space spanned by the first-order and unperturbed orbitals. A more
complete analysis of nonlinear behavior requires other quantities in addition to B k ( F ) , and it is not yet obvious to us which of these are conceptually most useful. E . Derivatives of the Charge Density with Respect to Field. Shown in Figure 3 are contour plots of the first derivative of the .rr-electrondensity with respect to field strength, weighted by the distance from the internuclear axis, rp('). The split-valence and minimal basis contour maps are similar in that both show local polarization of the CC triple bonds with only a small amount of charge transfer from one end of the chain to the other. The split-valence map extends out further from the internuclear axis giving an integrated polarizability contribution, (21), 54% larger than that for the minimal basis. The values, given in Table VIII, are cu,,(n) = 35.51 and 54.58 au for STO-3G and 3-21G bases, respectively. Similarly, the weighted third-derivative density is plotted in Figure 4. The split-valence basis gives a charge redistribution which is more diffuse and more intense than that given by the STO-3G calculation. The third-derivative function corresponds to electron transfer from the right half of the molecule to the left half, with large charge transfer between the interior carbon atoms. Figure 5 shows both rp(I)and rp(3)for octatetrayne, CBHZ,obtained with a 3-21G basis. Comparing the first-derivative function for this molecule with that for C4Hz(lower plot in Figure 3) shows that polarization is largely a local effect. The localized CC bond polarization is simply repeated for all four triple bonds in the longer chain molecule. A completely additive effect would correspond to (27) with unit exponent, close to the value X = 1.55 actually obtained when (27) is fitted to the 3-21G basis a,,(.rr) results in Table IX. The other fitting parameters are A = 8.97,
Th e Journal of Physical Chemistry, Vol. 93, No. 20, 1989 7129
Hyperpolarizabilities of Organic Molecules
TABLE X: Orbital Analysis of urrand yZznin the Cumulene System
tensor elem 4r)"
%A*') *
4CHdC %(u)d - 4 0 1 ,
,v
0
,
,
0 ,9'
. I
/B
,
0
I .
. I
,9'
0
.
!
0 0'
.
1 .
0 9'
.I
.
0
0'
I .
.
0
9'
.
. V .
0
/
cu,,(total) Yzzzr(")a YZZZA*'
Ib
Yzzrz(CH2)b Yzzz*(dd
Yzrzr(total)
C,HAC 16.51 8.35 5.92 30.77 -229.8 26.1 32.8 -170.9
CAHd' 62.32 12.89 18.40 7.66 101.27 -1482 -1127 914 145 -1551
GHAC 141.55 44.69 32.08 9.15 227.47 -6462 -8100 5240 499 -8822
CnHd' 259.64 113.45 38.19 10.50 422.37 -23340 -28500 13730 1384 -36726
Contributions from n *-type orbitals perpendicular to molecular plane. bContributions from (n - 1) *'-type orbitals lying in the molecular plane. cContributions from the pair of antisymmetric CH, u orbitals with the same symmetry as the in-plane r' orbitals. dContributions from the 4n + 1 remaining u orbitals. C A 3-21G orbital basis set is used at the STO-3G geometry. Carbon atoms lie on (I
."' 0
0
/q
'
0 /@
'
0
0
3'
0'
0
3'
0
b'
0
9'
0
Figure 5. Contour plots of the first (top) and third (bottom) derivatives of r-electron charge density with respect to field weighted by distance from internuclear axis, r # ) and r ~ ( ~Molecule ). is octatetrayne, CsH2, with field applied along CC bonds from left to right. Tick marks indicate nuclear positions. Orbital basis is 3-21G. Contour intervals are 0.1 au (top), 10.0 au (bottom).
B = 24.84,and 6 = 0.52. Comparing third-derivative functions for CsH2and C4H2shows that hyperpolarizability is a cooperative effect with substantial charge transfer along the entire length of ~ ~given ( in a Table ) the chain. A tit of the 3-21G basis ~ ~ ~results, IX,shows a 3.5 power dependence on chain length for a-bond hyperpolarizability. The fitting parameters are A = -63.5, B = 240.5,6 = 0.56,and A = 3.58 in this case. The a-orbital contribution to ~ z 2 2 2becomes negligible in comparison to the a contribution for the longer chains reported in Table IX, but its chain length dependence is remarkable. An acceptable fit of the 3-21G basis ~ ~ ~results ~ ~is obtained ( u ) with the parameters A = 53,B = -38, 6 = 1, and X = 2.65. This value of X is higher than we expected for nonconjugated u orbitals. Both the A and B parameters differ in sign from those for the a orbitals. In other words, the effect of the u system is to partially cancel hyperpolarizability contributions from the a system. This feature could be an artifact of a small basis set calculation and should be reinvestigated using diffuse orbital basis functions. F. Orbital Analysis Applied to Cumulenes. Bodart et aLs4 used minimal basis sets in their theoretical studies of the polarizability of cumulenes. All carbon-carbon bonds in a cumulene chain are double bonds with typical C=C double bond lengths of about 1.32 A for the early members of the series ethylene, allene, butatriene, etc. Longer cumulene chains are thought to spontaneously revert to an alkyne structure, Le., the =C=C=C= bonding gradually changes into =C-C=C-. It is of some interest to compare the electrical properties of these two different conjugated linear a systems. Table X reports longitudinal polarizability and hyperpolarizability for the even cumulenes, CbH4, n = 1-4. These molecules have point group symmetry Dzhwhich reduces to C2, in the presence of an electric field applied along the CC bond axis. The 6n 2 occupied canonical molecular orbitals are divided into four groups: (a) n a orbitals (perpendicular to the molecular plane), (b) (n - 1) a' orbitals (lying in the molecular plane), (c) two antisymmetric CH2 u orbitals (with the same point group symmetry as the a' orbitals), and (d) 4n 1 remaining u orbitals (completely symmetric under Czv). Results reported in Table IX show that the average in-plane a' orbital contributes slightly more to yzzZz than does an out-of-plane a orbital. The antisymmetric CH2bonds contribute much more than all other u orbitals combined, and their contribution increases
+
+
(54) Bodart, V. P.;Delhalle, J.; Dory, M.; Fripiat, J. G.; Andre, J. M. Opt. SOC.Am. 1987, 84, 1047.
the z axis.
4"'
rapidly with increasing chain length. In marked contrast to alkenes and alkynes, the a orbitais in cumulenes make negative contributions and the u orbitals make positive contributions to the longitudinal hyperpolarizability.
V. Conclusions Numerically stable computations of static polarizabilities and hyperpolarizabilities can be carried out for conjugated a-electron systems by the finite field method using polynomial fits of either energy or induced dipole moment. The F F method is, however, more difficult and less efficient than analytic derivative methods, where they are available, e.g., for SCF.'0-26 To the best of our knowledge FF is the only available technique for computing hyprpolarizabilities in low-order Mraller-Plesset (MP) theory. That quantum chemical formalism provides a practical, sizeconsistent approach to electron correlation effects in moderately large molecules. In comparing the results of fully coupled ab initio S C F calculations with those from ab initio SOS calculations using the same orbital basis, we find the two methods to be in fair agreement for a! but in poor agreement for 7.The small basis SOS results reported in Table I1 are in slightly better agreement with experiment than are the fully coupled S C F results, due to a fortuitous partial cancellation of errors due to methodology and basis set truncation. There is no theoretical justification for believing this to be a general result. A similar failure may or may not be inherent in other applications of the SOS method. In semiempirical calculations, for example, Coulomb and exchange integrals are replaced by parameters. Are these parameters, like the Coulomb and exchange integrals in the coupled HartreeFock operator, field dependent? When computing an energy derivative, should one include terms involving derivatives of these parameters with respect to field? It would be extremely useful to establish formal relations between empirical and ab initio parameters. In this context it is well to recognize the changing role of diffuse basis functions and the changing relative importance of u and a orbitals for the early members of a series of polyenes. When calibrating semiempirical parameters to be used in long molecules against available experimental data for small conjugated molecules, one should be aware that the dominant physical effects change significantly with chain length. Fully coupled S C F calculations with augmented basis sets reported in Tables VI and VI1 give polarizabilities and hyperpolarizabilities that are about 95% and 50%, respectively, of the accepted experimental values for ethylene and trans-butadiene. The work of Sekino and Bartlett suggests that a large part of the remaining correction will come from low order perturbation theory,24but further study of orbital basis set effects is probably of equal importance. In future studies of orbital basis set requirements we suggest that one try to establish whether a basis function is contributing mainly to the first- or second-order orbital response function. This information would be of use with modern
7130
J . Phys. Chem. 1989, 93, 7130-7134
analytic derivative methods developed according to the early ideas of Gready et This would significantly extend the size of molecules amenable to a b initio computation. In any case we anticipate that a b initio theory will play a significant role in the study of nonlinear properties of large conjugated molecules.
Acknowledgment. This work is supported by the Air Force Office of Scientific Research through Contract No. F4962087C0042 and by a grant of computer time on the CrayXMP at the National Science Foundation Pittsburgh Supercomputer Center.
Ground-State and Excited-State Electron-Transfer Reactions of Zinc Cytochrome c Edmond Magner and George MeLendon* Department of Chemistry, University of Rochester, Rochester, New York 14627 (Received: October 24, 1988; In Final Form: March 17, 1989)
Zinc-substituted cytochrome c, Zn cyt c, in which the heme iron atom is substituted by Zn(II), has been used in a number of basic studies of electron-transfer reactions of proteins. In this paper we report detailed studies of the ground-state and excited-state redox behavior of Zn cyt c which provide experimental estimates of key parameters including the oxidation potential (Zn cyt c)/Zn c, Eo z 0.8 V, and the (relative) reaction rates of Zn cyt c/Fe cyt c with a common reactant (corrected for reaction AG),kll(Zn)/k,,(Fe) = 5 i 1. The implications of these results for electron-transfer studies of proteins are briefly discussed.
Introduction Replacement of the metal atom in the active site of a protein has been used extensively to study the structure and function of proteins.' To this end a number of derivatives of cytochrome c,Z hemoglobin: and cytochrome c peroxidase4have been prepared and characterized. These derivatives serve as useful probes for the study of electron transfer (ET). For example, long-distance electron transfer has been studied in zinc/iron hybrid hemoglob i n ~and ~ ~in, a~variety of protein complexes such as cytochrome c/cytochrome c peroxida~e,~ cytochrome clcytochrome bS,5and zinc hemoglobin/cytochrome bS.3d A number of metal-substituted derivatives of cytochrome c (zinc, tin, manganese, cobalt) have been prepared and characterized* and used to study long-distance electron transfer between cytochrome c and its partner^.^.^ For example, zinc and free base porphyrin cytochrome c have been used to study the rate of electron transfer in the cytochrome c/cytochrome bS complex as a function of exothermicity.S A key requirement in such studies is the maintenance of the structure of the native protein; i.e., removal and substitution of the iron atom should not perturb the protein. A number of studies have been carried out to ensure that this requirement is fulfilled.2 For example, the circular dichroism spectrum of porphyrin cytochrome c (por c) is similar to that of (1) Hoffman, B. M. In The Porphyrins; Dolphin, D., Ed.; Academic Press: New York, 1979; Vol. VII, Chapter 9. (2) (a) Fisher, W. R.; Taniuchi, H.; Anfinsen, C. B. J . Eiol. Chem. 1973, 248,3188-3195. (b) Flatmark, R.; Robinson, A. B. In Structure and Function of Cytochromes; Okunuki, K., Kamen, M. D., Sekuzu, J., Eds.; University Park Press: Baltimore, 1967; p 318. (c) Moore, G. R.; Williams, R. J. P.; Chien, J. C. W.; Dickson, L. C. J . Inorg. Biochem. 1980, 12, 1-15. (d) Vanderkooi, J. M.; Erecinska, M. Eur. J . Eiochem. 1975,60, 199-207. ( e ) Eltis, L.; Mauk, A. G. J . Eiol. Chem. 1988, 27, 3103, Abstract No. 132. (3) (a) Peterson-Kennedy, S. E.; McGourty, J. L.; Kalweit, J. A.; Hoffman, B. M.J. Am. Chem.Soc. 1986,108,1739-1746. (b) McGourty, J. L.; Blough, N. V.; Hoffman, B. M.J. Am. Chem. SOC.1983, 105,4470-4472. (c) Simolo, K. P.; McLendon, G. L.; Mauk, M. R.; Mauk, A. G. J. Am. Chem. Soc. 1984,106, 5012-5013. (d) Simolo, K. P.; Stucky, G.; Chen, S.;Bailey, M.; Scholes, C. P.; McLendon, G. J . Am. Chem. SOC.1985, 107, 2865. (4) (a) Liang, N.; Mauk, A. G.; Pielak, G. J.; Johnson, J. A.; Smith, M.; Hoffman, B. M. Science 1988, 240, 311-313. (b) Cheung, E.; Taylor, K.; Kornblatt, J. A.; English, A. M.; McLendon, G.; Miller, J. R. Proc. Natl. Acad. Sci. U S A 1986,83, 1330-1333. (c) Taylor Conklin, K.; McLendon, G. L. J. Am. Chem. Soc. 1988,110,3345-3350. (d) McLendon, G. L.; Miller, J. R. J. Am. Chem. Soc. 1985, 107, 7811-7816. ( 5 ) Elias, H.; Chou, M. H.; Winkler, Jr. J. Am. Chem. Soc. 1988,110,429. ( 6 ) Armstrong, F. A.; Hill, H.A. 0.; Walton, N. J. Q.Reo. Eiophys. 1986, 18(3). 261-322.
0022-3654/89/2093-7 130$01.50/0
TABLE I: Quenchers k,/M-l
wencher
no. 1
H3C
s-l
E01V7 (uH7, 0.1 M P;) -0.78
3.0 X lo4
-0.65
3.4 x 105
-0.55
2.4 X lo6
-0.49
1.0 x 107
-0.18
2.6 x 107
WCH3 (CH:),
2cq N+\
(CH/),
ferricytochrome c , ~indicating ~ * ~ that the helical contents of the two proteins are nearly identical. By use of high-field N M R , it has been shown that the structure of the zinc cytochrome c is similar to that of the native protein.2c Vanderkooi" showed that addition of porphyrin cytochrome c does not stimulate the oxygen uptake of cytochrome c depleted mitochrondria, whereas cytochrome c does. However, por c does prevent the enhancement 0 1989 American Chemical Society
