9230
J . Phys. Chem. 1991, 95, 9230-9238
various intermolecular dimers in the previous paper16 and this work. With this approach, the overcorrection of BSSE by the Boys and Bernardi’s C P correction method has been evaluated for intermolecular potential energy calculations. First, we conclude that to compare with the total BSSE by the C P method, the amount of the overcompensation of BSSE at the second-order correlation potential energies generally is quite small at the minimum of the intermolecular potentials. Second, all results indicate that for accurate weak intermolecular potential energy calculations, the use of C P method is extremely important. In particular, a correction is essential at the electron correlation energy level because the BSSE is rather large while the BSSE at the Hartree-Fock level is relatively easy to reduce with practical large basis functions. In summary, an
accurate weak intermolecular interaction, such as van der Waals interaction between molecular, can be obtained with the C P method including electron correlations. Third, all results demonstrate that although the overestimation is small, the contribution of overcompensation of the BSSE from each subsystem depends on the number of occupied orbitals on the ghost atoms. Also the overcorrection from each subsystem as well as its BSSE strongly depends on the choice of basis functions, the electronegativity of atom, and the location and orientation of atoms. Acknowledgment. Work supported by Department of Energy grant from the Division of Chemical Sciences, Office of Basic Energy Studies.
Electromagnetic Properties of Molecules from a Uniform Procedure for Differentiation of Molecular Wave Functions to High Order Joseph D. Augspurger and Clifford E. Dykstra* Department of Chemistry, Indiana University-Purdue University at Indianapolis, 1125 East 38th Street, Indianapolis, Indiana 46205 (Receiaed: February 22, 1991)
Differentiation of the Schriidinger equation can be accomplished to arbitrarily high order for essentially any type of molecular wave function so as to obtain a variety of molecular properties analytically. The uniformity of derivative Schrijdinger equations makes possible computational implementation which is entirely open-ended with respect to order of differentiation. Techniques are developed here for electronic wave functions at the Hartree-Fock level and for one-dimensional vibrational wave functions in order to obtain electrical, magnetic, and electromagnetic response properties, and the derivatives of these with respect to geometrical parameters. These properties characterize the molecular response to static electric and magnetic fields, and representative values are reported for a number of small molecules.
Introduction Numerous molecular properties are derivatives of the energy of the molecular eigenstate. The dipole polarizability is the second derivative with respect to the strength of an applied electric field, for example. Magnetic susceptibilities, harmonic vibrational frequencies, and chemical shielding tensors are still other examples of properties that are formally defined as second energy derivatives. Hyperpolarizabilities, hypersusceptibilities, and so on are still higher order derivatives. The analytical evaluation of these energy derivatives follows from differentiating the Schrijdinger equation and solving the resulting equations, which are generally inhomogeneous differential equations. It is implicit that the energy and wave function, in the presence of a perturbation, are expressed as Taylor series, and so the coefficients of the series are the derivatives. Explicit, recursive expressions for these derivatives are given as well by Rayleigh-Schrijdinger perturbation theory (RSPT),’ and it is possible, in principle, to organize the solutions of the RSPT inhomogeneous equations in an open-ended manner. Alternatively, derivative quantities may be obtained by finite differences, but usually an analytical procedure is preferred because of the possible numerical errors in finite difference procedures and/or the substantial number of differences required for higher order property values. Organizing the differentiation of the Schriidinger equation can be done to achieve open-endedness in the computation of derivatives, and this is developed here. Two such methods have already been worked out, one for the differentiation of electronic selfconsistent-field wave function^^.^ and one for differentiation of ( I ) Bartlett, R. J.; Brandas, E. J. J . Chem. Phys. 1972, 56, 5467. Bartlett, R. J.; Shavitt, 1. Chem. Phys. Leu. 1977, 50. 190. (2) Dykstra. C. E.; Jasien. P. G.Chem. Phys. Lett. 1984, 109. 3 8 8 .
0022-3654/9l/2095-9230$02.50/0
numerical vibrational wave function^.^,^ The generalization to complex operator problems and the general use of the “2N 1” rule are incorporated here. These developments make it possible to calculate essentially any conceivable property, including certain important mixed derivative properties that are rarely calculated. To demonstrate the capability, applications are reported for problems of electric and magnetic field effects for both equilibrium structures of small molecules and individual vibrational states of diatomics. As experimental techniques evolve to probe molecules in greater detail and on shorter time scale^,^ the detailed molecular response to fields and other perturbations becomes more at issue and increasingly measurable. Purely theoretical studies have pointed out the possibility of tuning vibrational transition frequencies and transition strengths in small molecules, at least, by application of electric fields6 It is the vibrational response to an applied field that comes into play. As another case, the subtle variation of chemical shielding with respect to an applied electric field has recently provided a basis to model distal ligand effects in carbonmonoxyheme proteins.’ The ligands electrically perturb the heme CO, and for a series of myoglobins, hemoglobins, and peroxidases, the chemical shifts for I3C and ” 0 directly reflect the electrical environment. Overall, the simultaneous interactions of electric and magnetic fields plus the role of vibration underlie many phenomena, and complete understanding calls for a theo-
+
(3) Dykstra, C. E. Ab Initio Calculation of the Structures and Properties of Molecules; Elsevier: Holland, 1988. (4) Dykstra, C. E.; Malik, D. J . J . Chem. Phys. 1987, 87, 2806. ( 5 ) Zewail, A. H.; Bernstein, R. B. Chem. Eng. News 1988, 66, 24. ( 6 ) Malik, D. J.; Dvkstra. C. E. J . Chem. Phvs. 1985. 83. 6307. ( 7 ) Augspurger. J b , Dykstra, C E , Oldfiild, E J Am Chem SOC 1991. 113. 2441
0 1991 American Chemical Society
The Journal of Physical Chemistry, Vol. 95, No. 23, 1991 9231
Electromagnetic Properties of Molecules retical analysis of the corresponding response properties. With the growing need for property and response information, there has come a tremendous advance in the methodology for energy derivatives of a b initio electronic wave functions. Much of this advance has been in the past fifteen years or so [see the papers collected in ref 81. The development we report here exploits different aspects of available methods and uses a simple computational organization to achieve the ultimate goal, a completely open-ended, analytical method for all properties. Though the formal theory is straightforward and already understood, we include as Appendix a general derivation in a concise form that is the most suggestive of the computational implementation we have adopted. Among the earliest work on methods for energy derivatives was McWeeny's use of a perturbative expansion of the HartreeFock equations to obtain analytical derivatives9 Stevens, Pitzer, and Lipscomb developed a perturbative expansion for the basis set form of the Hartree-Fock equations to obtain certain electrical and magnetic properties in 1963.1° Langhoff, Karplus, and Hurst tested different forms of perturbed Hartree-Fock, leaving the "coupled-perturbed Hartree-Fock" (CPHF)" as the complete form used almost exclusively today for first and second derivatives. Gerratt and Mills12generalized the formulation of Stevens et a1.I0 so that geometrical derivatives could be obtained with CPHF. h l a y pioneered the contemporary treatments of energy gradients and their use in finding potential surface minima." More recently, the number of reports of methods, approaches, and computational schemes has grown sharply, and among the highlights are the explicit expressions for several low orders of energy derivatives given by McWeeny and c o - ~ o r k e r s , l ~the ~ ' recognition ~ of the uniformity of the first- and second-order equations by Dupuis,16 and the derivation and early use of "2"' and "2N 1" rules by Nee, Parr, and Bartlett." Pulay has carefutly reviewed and summarized the analytic differentiation formalism and methodology through the mid- 198Os.l8 In the past ten years, there has been considerable progress in achieving analytic differentiation of correlated electronic wave functions with respect to geometrical parameters, starting essentially with the configuration interaction (CI) derivative methodologies of Brooks et aI.l9 and of Krishnan, Schlegel, and Pople.20 Formal expressions for up to fifth-order energy derivatives were reported by King and Komornicki,21while Handy and Schaefer22demonstrated that 2N and 2N 1 rules apply to MCSCF and SDCI wave functions. First- and second-order differentiation capability exists for MCSCF, MBPT, CI, and coupled cluster (CC) wave functions (refs 8 and 23-26, and
+
+
references therein). It is fast becoming routine to find many energy derivatives of S C F wave functions at third and fourth order, and even analytically determined fifth derivatives are being reported.27-29 The derivative Hartree-Fock (DHF) method,2 which is open-ended with respect to the order of differentiation, has been applied to finding dipole hyperpolarizabilities all the way to the ninth d e r i ~ a t i v eand , ~ ~ this report presents a more global development and implementation of DHF.
General Approach for High-Order Differentiation Organization of the Differentiation with Respect to Multiple Parameters. Properties that are often of interest involve differentiation with respect to structural parameters of the molecule, to external fields and field gradients, to internal fields due to nuclei, or to parameters introduced via some model Hamiltonian. Clearly, a general treatment must provide for simultaneous differentiation with respect to multiple, independent parameters. Uniformity in the derivative Schrodinger equations can be exploited by first organizing the parameters as an ordered set, such as {a,6, c, d, ...),and then by using ordered integer lists to designate specific differentiation operations. An integer list associates one integer with each parameter, with the integers ordered the same as the parameters. Then, a given integer specifies the level of differentiation for the associated parameter. For instance, a mixed third-order differentiation operation represented by the list [ 1 0 2 0 ...I is taken to be differentiation with respect to the first parameter, a, and second differentiation with respect to the third parameter, c. Organization of the derivatives is essential for implementing any procedure that is open-ended with respect to the order of differentiation and which exploits the uniformity of the derivative equations. A workable scheme is to order the unique derivatives according to their integer lists, and in particular, an anticanonical arrangement (where the "odometer" is driven from the left) is convenient. Ordering of the lists for the case of four parameters illustrates the sequencing:
[oooO] [ 10001
(no differentiation-zero order in all parameters) (differentiation with respect to the first parameter)
[OlOO] [OOIO] [OOOI] [2000] [ 11001 [IOIO]
(second differentiation with respect to the first parameter) (mixed second differentiation)
[ IOOl]
[0200]
[OIIO] [OIOl] (8) Jorgensen, P., Simons, J., Eds. Geometrical Deriuatiues of Energy Surfaces and Molecular Properties; Reidel: Dordrecht, Holland, 1986. (9)McWeeny, R. Phys. Rev. 1961,126, 1028. (IO) Stevens, R. M.; Pitzcr, R. M.; Lipscomb, W. N. J . Chem. Phys. 1963, 38. 550. ( I I ) Langhoff, P. W.; Karplus, M.; Hurst, R. P. J . Chem. Phys. 1966,44, 505. (12) Gerratt. J.; Mills, 1. M. J . Chem. Phys. 1968,49, 1719. (13) Pulay, P. Mol. Phys. 1969,17, 197. (14)Dodds, J. L.; McWeeny, R.; Raynes, W. T.; Riley, J. P. Mol. Phys. 1977,33, 61 I . (15) Dodds, J. L.; McWeeny, R.; Sadlej, A. J. Mol. Phys. 1977,34,1791. (16)Dupuis, M. J . Chem. Phys. 1981,74,5758. (17)Nee,T.-S.; Parr, R. G.; Bartlett, R. J . J . Chem. Phys. 1976,64,2216. (18)Pulay, P. Adu. Chem. Phys. 1987,69,241. (19)Brooks, B. R.; Laidig, W. D.; Saxe, P.; Gcddard, J. D.; Yamaguchi, Y.; Schaefer, H. F. J . Chem. Phys. 1980,72,4652. (20)Krishnan, R.; Schlegel, H. B.; Pople, J. A. J. Chem. Phys. 1980, 72,
4654. (21)King, H.F.; Komornicki, A. J . Chem. Phys. 1986,84, 5645. (22)Handy, N. C.; Schaefer, H. F. J . Chem. Phys. 1984,81, 5031. (23) Watts, J. D.; Trucks, G. W.; Bartlett, R. J. Chem. Phys. Lett. 1989, 164,502. Stanton, J. F.;Watts, J. D.; Bartlett, R. J. J . Chem. Phys. 1991, 94. 404. (24) Koch, H.; Jensen, H. J. Aa.; Jorgensen, P.; Helgaker, T.; Scuseria, G. E.; Schaefer, H. F. J . Chem. Phys. 1990,92,4924.Scuseria, G . J . Chem. Phys. 1991, 94, 442. (25)Urban, M.; Diercksen, G. H. F.; Sadlej, A. J.; Noga. J. Theor. Chim. Acta 1990, 77,29.
[0020] [OOI11 [00021 [3000] (third differentiation with respect to the first parameter) [2 1001...
We will use & to designate a single arbitrary integer list and will use a (no arrow) to be the sequence number in the integer list. Thus, in the case above when a = 3, the corresponding integer list is Z = [0 1 0 01. The integer lists are used explicitly for the logic operations in a uniform derivative procedure. Logic operations require definitions of addition and comparison of integer lists. Addition is defined as numerical addition of corresponding elements, e.g., [ 2 1 0 01 + [0 1 1 01 = [ 2 2 1 01. For this reason, the integer lists may be thought of (and are designated as) vectors. Integer lists may be compared as in the Kronecker delta, Sa$, in the same way ~ _ _ _
~______
~~~________
~~~
(26) Pal, S.Phys. Reu. A 1990, 42,4385. (27)Sekino, H.;Bartlett, R. J. J. Chem. Phys. 1986,85, 976. (28) Karna, S. P.; Dupuis, M. IBM Tech. Rep. No. KGN-211, 1990. Karna, S. P.; Prasad, P. N.; Dupuis, M. J . Chem. Phys. 1991, 94, 1171. (29)Allen, W. D.; Yamaguchi, Y.; Csaszar, A. G.; Clabo,Jr., D. A,; Remington, R. 8.; Schaefer, H. F. Chem. Phys. 1990. 145, 427. (30)Liu, S.-Y.; Dykstra, C. E.; Malik, D. J. Chem. Phys. Lett. 1986, 130,
403.
9232 The Journal of Physical Chemistry, Vol. 95, No. 23, 1991 as two geometrical vectors. If_every element in G is the same as the corresponding element in 0, then the lists are the same. If there is a difference in any element, then Je,8 = 0. We define n ( a ) to be the order of the a derivative. n(a) =
eai
(1)
i
where G = [cu,, cy2, ,..I. The array n distinguishes first derivatives from second derivatives, and so on. Of course, in view of the anticanonical ordering of the integer lists, all first-derivative lists precede all second-derivative lists, and so on. Thus, instead of using the array n, it is convenient to use short arrays that point to the beginning and ending sequence number of the interger lists for a given order of differentiation. We call these arrays b and e, respectively, and continuing with the example given above, the elements of these arrays are b(l) = 2 e(1) = 5
Augspurger and Dykstra largest order of differentiation is sizable, perhaps greater than IO, eq 3 may involve ratios of much larger integers than will eq 4 and this may pose numerical complications. Overall, the integer arithmetic for organizing the differentiation operations is going to be minor in relation to the whole process of computations that yield the desired derivative values. Formal Differentiation of the Time-Independent Scbriidinger Equation. The notation we use here is that an a superscript means that differentiation has been carried out with respect to the set of parameters specified by a. A zero subscript means that the item has been evaluated for the equilibrium, reference, or zeroorder choice of the parameters. The Schriidinger equation with unspecified values for the embedded parameters is written here as ( H - E)* = 0, whereas the Schriidinger equation with specific values of the parameters being either zero or some reference value is (Ho - Eo)*, = 0. Formal differentiation with respect to the first parameter, a, and rearrangement yields ( H - E)*" = -(Ha - E"*
b(2) = 6 e ( 2 ) = 15 b(3) = 16 e ( 3 ) = 35
*;
and so on
For any particular differentiation operation a, all other operations of the same order are sequenced in the integer lists between the positions b(n(a))and e ( n ( a ) )inclusive. For computer implementation, an integer list is constructed explicitly for all derivatives sought. This list depends on the number of parameters, the order of differentiation, and on whether differentiation is to be carried to the same level for each parameter. The arrays b and e must be set up, but the array n is simple enough that it need not be set up explicitly. A minor routine is required that compares two vectors, Le., that finds bs,+ Another routine is needed for computing certain products of binomial coefficients. The following notation is helpful for these products. N
G!
= nai! i= I
N is the number of parameters. The factors that are needed in an open-ended differentiation procedure are designated t and are given by i= I
-
E& a first-derivative property, and are obtained by solving this inhomogeneous equation after the parameters have been set to zero or their reference values: Continued differentiation of the Schriidinger equation to arbitrarily high order yields
The summation indices p and y are independent, but only because of the presence of the 6a,8+y factor. This is a useful form for computation because it casts the expression as a loop-structured process, and is convenient relative to carrying out logic operations in a direct search for contributing (Hp - E@)@?terms. In expressions that follow, extension (or other change) of summation indices where there are suitable &factors is done freely. With parameters set to their equilibrium or reference values, the energy derivative, E:, can be found by multiplying eq 8 on the left by Po,integrating, and rearranging.
Et = (\kolHa190) +
(iGi)! f ( G & ,...,tin)=
I e(n(al-1) e(n(a)-1)
(3)
c
fiG>
i=1
i= I
Th_ese factors are used only when, for example, G, 8, T, and ( a + T) each show up in the set oQnteger lists; that is, the factor will not be needed if the sum (5 + p + T ) is not in the list. Thus, it is sufficient to compute and store a table of values, G!, for all a's, and then evaluate any t factor as needed by looking up the necessary values in that table and applying eq 3 . The alternative to a table look-up to find the t factors on a parameter-by-parameter basis. This is best done with a stored table of binomial coefficients:
+
(f+')
(6)
2
p=2
c
7=2
J e , ~ + ~ + ~ f ( ~ , p , ~ )-( *EglW ~ l H g (9)
This result expresses the N + 1 rule, that (A'+ 1)-order energy derivatives may be found from only N and lower order wave function derivatives. The trivial summation over t has been introduced in order to write out 2N and 2N + 1 rules. These develop by reducing the range of y at the expense of increasing the range of t . The following steps uniformly eliminate the y derivatives of highest order while adding the next highest order t derivatives. First, the summation over y in eq 9 is partitioned.
(k + l ) !
=
7
These are used to form the/ factors used in previous explanations of derivative Hartree-Fock theory,2 (4)
Equation 8 relates an N-order derivative of \k to a summation over all N - 1 and lower order derivatives, and it may be used31 to replace the second summation in eq loa:
--c
These factors can be multiplied together to yield equivalent t factors:
(=I
e(n(o)-l) e(n(a)-2)
E 8=2
~ ( 1 )e(n(a)-l) e ( n ( o ) - l )
c +c c i=2
7=2
812
c
y=l
(lob)
And then the summations may be recombined.
N
f(a,P*T)= IIjIaiybi,yi)
1
(5)
1=1
The vector designation for a, @,and y on the left will be suppressed from here on, because it is clear that t depends on the elements of the integer lists. Storing binomial coefficients and forming f factors according to eqs 4 and 5 is preferable to storing the factorials of eq 2 and using eq 3 in only one small way. When the
~~~~~
______
~~
(31) Substitutions employ eq 8 twice. Working with just the terms of eq 9 where y ranges between b(n(a)- 1 ) and e ( n ( a ) - l ) , which are the highest order of y , (3 will range only over first-order derivatives. Using the adjoint
form for the terms, substitution is made following eq 8, which increases e by I order, while reducing fl to be just I , as in (H- E ) 9 " = -(He- E')*. Next, eq 8 is used to increase B while decreasing the range on y . The process may be repeated for the next highest order y terms.
The Journal of Physical Chemistry, Vol. 95, No. 23, 1991 9233
Electromagnetic Properties of Molecules dI) r(n(a)-l) O W - 2 )
- c 812c c e=l
(10c)
1-1
At this point, E; may be evaluated by only knowing the wave function derivatives through 2 orders less than a (Le., n(a)- 2). The summation changes may be carried out N times where the order of a is 2 N (or 2N + 1) (Le., n(a) = 2 N or 2 N + 1). This will increase the limit on c to N a n d reduce the limit on y to N - 1 (or N). The result is the 2N or 2 N 1 rule for energy derivatives (or perturbative corrections).
+
n(a) = 2N: E; = (901H@'o) 4 0 r(2N-I)
e(N-I)
C 8'2C ('I
C y=l
+ 1:
Eg =
n(a) = 2N
+
b;+g+&B,y)(%I~
- E6198 (1 l a )
TABLE I: Molecular Geometries for the Property Cnlculrtionrr (A) molecule atom X V z H, H1. H2 h0.370424 0.0 0.0 H ~ N H -1.623 922 7 0.0 0.0 -0.558 022 7 0.0 C 0.0 N 0.595 077 4 0.0 0.0 N 0.067 661 0 0.0 0.0 HI -0.3133669 0.937 529 6 0.0 -0.468 764 8 h0.81 I 924 4 H2, H3 -0.3133669 0.0 0.0 0.0 C CH4 HCCH H2CO
co
H I , H2 H3, H4 c1. c 2 HI: H2 H I , H2 C
k0.629 677 7 h0.629 677 7 f0.601 k1.661 -0.541 8136
0
1.210 -0.644 483 5 0.483 516 5 -0.870615 5 0.046 184 5 0.0 k0.854 362 9 f0.854 362 9 0.0 k1.16 0.065 569 2 -0.5203130
C
0
The Appendix presents these equations for the specific problems of S C F wave functions in a basis set expansion and of vibrational wave functions. Complex Perturbations. It is possible to encounter problems where the operator associated with the perturbation is complex, and a key example is the interaction with an external magnetic field. In treating the response to these properties, it could prove costly in terms of disk space and memory to merely declare all variables as complex, though this would succeed. However, it should usually be possible to separate a complex perturbational operator into real and imaginary parts, and so the complication is to manage the treatment of the imaginary operators using a computational implementation with strictly real numbers. It has been pointed that this may be accomplished by factoring i out of the operator and combining it with the parameter of differentiation. Thus, if there exists some perturbation of the form AH,where H is real, then the procedure is to differentiate the energy with respect to X = iX rather than with respect to A. The energy and wave function derivatives with respect to X are found by scaling with the appropriate power of i.
That is, the desired derivatives are obtained by multiplying a purely real derivative value by i raised to the sum of the powers of the purely imaginary perturbations in the specific derivative. Sometimes a basis set is constructed with complex function^,^*^^ and usually these are perturbation dependent. For magnetic properties, the perturbation dependence is usually on the ga~ge,'~ and then the differentiation has the same complication as encountered for geometrical derivatives.
Results and Discussion We have implemented open-ended differentiation procedures (see Appendix) and have applied them to a number of small molecules. The basis sets used for all atoms except silicon and the hydrogen atoms in H2 are those described in ref 36. These basis sets, which have been termed ELP (electrical properties), are triple-zeta, triply polarized bases augmented with diffuse valence functions. Because we find that gauge insensitivity of the (32) Lazzeretti, P.;Zanasi, R. J . Chem. Phys. 1977, 67, 382. Schindler, M.; Kutzelnigg, W. J . Chem. Phys. 1982, 76, 191. (33) Pople, J. A.; Krishnan, R.; Schlegel, H.B.; Binkley, J. S. Int. J . Quantum Chem. 1979, S13.225. (34) Dodds, J. L.;McWeeny, R.; Sadlej, A. J. Mol. Phys. 1977, 34, 1779. Sadlej, A. J. Chem. Phys. Leu. 1977,47,50: Mol. Phys. 1977,34,731: Theor. Chim. Acra 1978,47, 205. Szalewicz, K.; Adamowicz, L.; Sadlej, A. J. Chem. Phys. Lett. 1979, 61, 548. (35) Ditchfield, R. J. Chem. Phys. 1972,56,5688; Mol. Phys. 1974,27, 789. Kutzelnigg, W. Isr. J . Chem. 1980, 19, 193. Schindler, M.; Kutzelnigg, W. J. Chem. Phys. 1982, 76, 1919. (36) Dykstra, C. E.;Liu, S.-Y.;Malik, D. J. Adu. Chem. Phys. 1989, 75, 37.
HF
H
SiHd
F Si
co2 H20
H I , H2 H3, H4 C 01,02
0 HI. H2
0.0
0.629 677 7 -0.629 677 7 0.0 0.0 h0.959 605 2
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.854 362 9 -0.854 362 9 0.0 0.0 0.0 k0.756 950 3
f0.629 677 7 70.629 677 7
0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 *0.854 362 9 70.854 362 9
0.0 0.0 0.0 0.0
chemical shifts can be increased by partly uncontracting the p functions on first-row atoms, a (6p/5p) contraction has been employed instead of the usual (6p/3p) contraction of a triple-zeta valence set. For the first-row atoms, the total composite bases may be described as (1 ls8p3d/7s7p3d), except for CO, where the basis of ref 7 was used which employed four sets of d functions (exponent values of 2.4,0.6,0.15, and 0.04). For silicon, the basis used was the valence (12s9p/6s5p) set of McLean and Chandler?' augmented with diffuse s and p functions (exponents 0.03 and 0.02),a set of 3d functions (1.6, 0.24, 0.035), and a single f function (exponent 0.4), for a final basis of (13sIOp3dlf/ 7s6p3dlf). The hydrogen basis for Hzwas that used in ref 38, which consisted of 19 functions (7s3pld/5s3pld). In all cases, the molecular axis or principal symmetry axis is the x axis. The geometries used are given in Table I. All chemical shieldings and related magnetic properties were computed with the gauge center of the external magnetic field fixed at the nucleus of interest. Table I1 gives calculated values for the isotropic chemical shieldings. The chemical shielding is obtained as the second derivative of the energy with respect to the nuclear magnetic moment and an external field, ua6 = dB,. Absolute paramagnetic chemical shieldings can be derived from spin-rotation constants, which must then be combined with theoretical diamagnetic shielding values to get absolute chemical shieldings which can be compared with our computed values. Alternatively, once a particular absolute shielding is known, it can be used to convert chemical shifts to absolute shieldings. Our computed shieldings are in good agreement with the experimental values, within a few ppm for hydrocarbons. However, for compounds containing nitrogen, the errors are a bit larger, but less than IO ppm. For oxygen and fluorine containing compounds the errors are about 10-20 ppm. The second-order dependence of the shielding on the external magnetic field represents the nonlinearity in the shielding; it is the second derivative of the chemical shielding tensor with respect to external magnetic field components. For a freely tumbling molecule, which experiences an external field such that B, = B = B, = '/'B, the isotropic shielding can be written a = + a2ii/aBZ(see footnote to Table 11). The largest nonlinearity in the shielding is for the oxygen in HzCO. Since field strengths used in N M R spectrometers range to around 10 T , which is equivalent to a magnetic field of 4.25 X in atomic units,3g
I/$''
(37) McLean, A. D.; Chandler, G. S. J . Chem. Phys. 1980, 72, 5639. (38) Augspurger, J. D.; Dykstra, C. E.J . Chem. Phys. 1988, 88, 3817.
9234
The Journal of Physical Chemistry, Vol. 95, No. 23, 1991
TABLE II: Chemical Shielding and Nonlinear Components for Several Small Molecules (in ppm, nu)' B a2a/as2 nucleus' molecule t+ BP 195.9 -0.1 1 "C CH, 296.3 -100.4 195.1' 0.05 HCCH 321.4 -201.6 1 19.8 117.2' 0.05 76.3 326.5 -250.2 HCN 82.1' -17.9 0.06 326.5 -344.4 co 3.2 f 0.3' -324' -1.9 338.8 -340.7 0.36 HZCO -7.7' I4N 261.8 354.1 -92.3 -0. IO NH3 264.59 -9O.lf 0.16 -45.0 HCN 378.8 -423.8 -3 7e -4 18' 416.0 -90.7 170 325.3 -0.07 H20 346 0.19 -79.2 co 444.8 -524.0 -42.3 f 17f -485 f 17f I .76 -469.0 452.3 -921.3 HJO -830 f 1 0 6 -375 f 1 0 6 -312 f 2 6 I9F HF 482.2 -70.6 411.6 -0.03 419.7 f 0.3d -63* O.OOh 33.7 'H HCCH 99.1 -65.4 32.Y 99.5 -66.7 32.7 O.OOh HCN 334 27.7 112.6 -84.9 0.01 H2CO 18.3 f 2.09 -93.2 f 1.28 108.4 -70.0 HF O.Oh 38.4 -79.7 f 0.3d 29.2 f 0.5d O.Oh 102.4 -65.7 36.7 H20 26.0 f 0.3d O.Oh 95.5 -60.4 35.1 N H3 31.2 f 1.09 -65.38 87.5 -53.8 O.Oh 33.7 CH4 30.9 f 0.3d 32.2 -5.6 O.Oh 26.6 H2 26.4 f 0.6'
'Experimental values are given on the second line of each entry. ba2a/aB2 = ( I / 9 ) ~ y a u m , / a B , 2 . cReference 41. The value for H 2 C 0 was derived from 6(C) = -1 with reference to CS2 by P. C. Lauterbur reported in ref 40 and o(C) = -6.7 for CS2.41 dReference 42. eReference 43. /Reference 44. #Reference 45. hValue less than 0.01. 'Reference 46. /Reference 47. Values derived by using the value of a ( H ) for CHI from ref 42. 'Gauge origin for each chemical
shielding is the nuclear center. the nonlinear effects on the shielding are negligible, amounting to - I 0-9 ppm. A more significant sensitivity of the chemical shielding is revealed by the results for several molecules in Table 111. It gives the chemical shielding derivative with respect to an applied electric field and field gradient. Buckingham formally expressed the chemical shielding as a power series in applied electric field4*
where d")is an n
Augspurger and Dykstra TABLE 111: Shielding Polarizabilities for Several Small Molecules (in ppm/au) molecule nucleus A B C H2 -50.3 -1 87.7 5.8 H2 HCCH H2 -67.2 -13.6 -109.5 HCN H -54.1 -173.2 89.5 H2CO 0.1 HI -474.5 5.0 HF H -81.5 328.5 61.8 HI -47.3 -80.4 3.5 H20 H1 -27.7 -278.3 -18.1 NH3 HI -45.1 -114.4 -2.1 CH4 H4 -35.8 SiH4 -1.7 -276.4 C -697.4 H2CO 3874.9 746.0 co C -374.5 -535.5 532.7 HCN C -422.6 1502.9 512.1 c2 -733.9 -1 106.9 HCCH -35 1.8 0.0 C -359.7 -54.1 CH4 N 50.8 1613.4 333.3 NH3 -1 1336.1 N 1910.1 HCN 603.6 0 401.1 -4339.3 -190.8 H2O co 0 1526.7 -5906.1 1044.0 H2CO 0 7019.0 -70843.0 1377.8 -741.7 HF F 636.5 -8486.1 Si 0.0 6939.5 SiH4 456.5 TABLE IV: Dipole Moments, Polarizabilities, and Hyperpolarizabilities and Their Geometrical Derivatives for HF, CO, and CO,"
co
co,
pxJJ.x
HF -0.766 -5.613 8.323 -264.0
-0.103 -14.42 3 1.03 -1 166.0
0.0 -23.63 0.0 -791.0
ar:iaR ac,xlaR apx.x.xlaR apx,.x,laR
-0.406 -5.543 27.09 -523.8
-1.025 -9.719 6.216 -596.2
0.0 -1 0.48 0.0
(apx,laR)lp,,x (apx.x,laR)lpx,, (apx.x,.xlaR)lpx,x,,
100% 325% 200%
67% 20% 51%
44% 0% 28%
r: px.x p x A X
-224.9
'Values are in au. The polarizabilities are reported using the P convention described in ref 31. P,, corresponds to -axx,P,,,, to &, and P,,,, to -yxxxx. For H F and CO, R is the internuclear separation coordinate, while for C 0 2 it is the symmetric stretch coordinate. nonzero elements of these tensors for various symmetries has been worked out by Raynes and R a t ~ l i f f e . Electric ~~ field gradient effects may be considered in the same manner. So, for a molecule being perturbed by an electric field and field gradient whose only nonzero elements remain fixed along the molecular axis (throughout taken to be the x axis) as it freely tumbles in the applied magnetic field, we can write the isotropic shielding as a function of the electric field and gradient =
8,
+ A E , + Y ~ B E -+, ~CE,, + ...
+ 2 rank tensor. The number of independent,
(39) Dykstra, C. E.; Augspurger, J. D.; Kirtman, B.; Malik, D. J. Rev. Comput. Chem. 1990, I , 83. (40) Neumann, D. B.; Moskowitz, J . W. J . Chem. Phys. 1969, 50, 2216. (41) Jameson, A. K.;Jameson, C. J. Chem. Phys. Lett. 1987, 134, 461. (42) Hinderman, D. K.;Cornwell, C. D. J . Chem. Phys. 1968, 48,4148.
(43) Jameson, C. J.; Jameson, A. K.; Oppusunggu, D.; Burrell, P. M.; Mason, J . J . Chem. Phys. 1981, 74, 81. (44) Wasylishen, R. E.; Mooibroek, S.; McDonald, J . B. J . Chem. Phys. 1984, 81. 1057. (45) Kukolich, S.G.J . Am. Chem. SOC.1975, 97, 5704. (46) Myint, T.; Kleppner, D.; Ramsey, N. F.; Robinson, H. G.Phys. Rec.
Lett. 1966, 17. 405. (47) Schneider, W. G.; Bernstein, H. J.; Pople, J. A. J . Chem. Phys. 1958, 28, 601. (48) Buckingham, A. D. Can. J . Chem. 1960. 38, 300.
Note that commas in the subscripts of CT are used to separate the elements of the electric field or field gradient from the first two subscripts which refer to the components of the magnetic moment is an element of a third-rank and magnetic field. Thus, CT&!,, tensor, the derivative of crXx with respect to Exx. Analogous to defining the polarizabilities as the various derivatives of the multipole moments with respect to elements of an electric field, we refer to A , B, and C as various shielding polarizabilities.' A ~
(49) Raynes, W T , Ratcliffe. R Mol Phys 1979, 37, 571
The Journal of Physical Chemistry, Vol. 95, No. 23, 1991 9235
Electromagnetic Properties of Molecules TABLE V: Geometrical Derivatives of Chemical Shieldings and Shielding Polarizabilities (in DPOI au) of HF and CO
hydrogen fluoride H F 8
aalaR A
BAfaR
38.37 35.81 -81.52 76.52
411.6 230.3 636.5 2689.2
carbon monoxide C 0 -17.93 305.9 -374.5 14.99
-79.23 -604.8 1526.8 2871.52
very important result from this body of results is that the shielding polarizabilities are highly dependent on the chemical environment. For instance, the magnitude of A values ranges from 375 to 734 for multiply bonded carbon. As well, we see that the B values for a proton connected to a triply bonded carbon depend sizably on whether that carbon is bonded to a nitrogen or a carbon (HCCH vs HCN). The proton chemical shifts, on the other hand, are similar for both molecules, 33.7 ppm vs 32.7 ppm. SO,the "signature" of a chemical environment that occurs for chemical shielding is not found for the shielding polarizability. It is sensitive to subtler features of the electronic structure. Table IV provides another illustration of the computational capability for derivatives. Here, the geometrical derivatives of the dipole moment and polarizabilities have been obtained for three small molecules. The ratio of the property derivatives to the equilibrium property values provides a relative measure of the sensitivity of the properties to the internuclear coordinate. For a given change of 0.1 au in the internuclear separation, the polarizability changes by 4-1 0'30,the hyperpolarizability by 0-30%, and the second hyperpolarizability by 3-20%. The HF single bond has a more sharply changing set of response properties than the multiply bonded CO or COS. In Table V, chemical shielding, shielding polarizabilities, and their geometrical derivatives are given for H F and CO. For linear diatomics, the primary effect of an axial electric field is to shift charge along the axis from one atom to the other. This shift in charge density leads to a diminishment of shielding on one atom and an enhancement on the other. For a homonuclear diatomic, the A values for the two centers would be equal, but opposite in sign. In the case of C O and HF,the 0 and F have positive A values, whereas H and C have negative values that reflect their orientation with respect to the applied field. The magnetic properties of Table V are more uniformly sensitive to internuclear separation than the electrical properties in Table IV, with the exception of the derivative of the shielding polarizability of C. However, this is a case of a rapidly changing derivative, as evidenced by the derivative values of 249.0 at +0.03 A from the equilibrium separation and -152.7 at -0.03 A from equilibrium. The general differentiation approach is applied to the electromagnetic properties of the first three vibrational states of C O by means of the DNC method (see Appendix). Results are given in Table VI. The potential energy was obtained from large basis set calculations using well-correlated, coupled cluster wave functionss0 a t 12 points spanning almost 0.5 A about the equilibrium separation, with the farthest out points 16000 cm-' above the minimum. The electromagnetic properties were also obtained at several points, and then the potential and properties fitted to low-order polynomials to be represented as functions of the internuclear separation. Comparing the equilibrium electrical properties of CO givei: in Table IV to the vibrational state properties, the nuclear contributions to the dipole moment and polarizability are small, yet are quite significant for the hyperpolarizability. These results are similar to those of Adamowicz and Bartlett's study of HFS5' Most properties change little for the first three vibrational states, with the exception of the off-axis elements of the chemical shieldings and their associated shielding polarizabilities. An axial component of a magnetic field may be thought of as inducing an electric current about the molecular axis. Likewise, an off-axial component of a field will act on the
TABLE VI: Electromagnetic Properties" of the VibnHonrl States
of co
v=o
property PXb PXS PYY PXSS PXYY . x m .yy(C)
W ) ~XX(0) .xxs(C) U Y Y m
axx,(O) g Y Y m gxx,(C) gyy.Zx(C) ~XXSX(0)
v = 2 -0.004 -1 5.39 -1 1.38 13.98 2.53 270.60 -189.39 -36.06 410.66 -384.16 -1 19.22 -41.93 -814.01 39.21 1933.51 -44.94 985.07 -29.26 1933.51
%rX(O) Electrical properties in atomic units, chemical shielding, and shielding polarizabilities in ppm atomic units. See footnote to Table IV for conventions describing electrical polarizabilities. b Apositive dipole moment has the sense of C-0'.
electron distribution along the axis. Since the electron distribution along the axis tends to depend more strongly on the internuclear separation than does the distribution along any off-axis slice, it would not surprising that the off-axis component of the shielding would be more sensitive to the separation distance as is seen. This leads to "chemical shifts" of 7 ppm for 13Cand 15 ppm for I7O for the first excited vibrational state of CO with respect to the ground state. Acknowledgment. This work was supported, in part, by a grant from the Physical Chemistry Program of the National Science Foundation (CHE-8721467).
Appendix Basis Set Representations. The general differentiation of the Schrijdinger equation may be carried out with a basis set representation such as that used for configuration interaction (CI) wave functions. The general equation is (H - E)C = 0 (18) where H is the matrix representation of the Hamiltonian in an orthonormal basis, E is the diagonal eigenvalue matrix, and C is the matrix of expansion coefficients for the states of the system. General differentiation yields a e(n(a)-l)
(H - E)C" = -E 82'
&$+,t(P,y)(H@ - E@)C7 (19)
7-1
This corresponds to eq 8, and in a like manner the energy expressions of eqs 9 and 11 may be written. The N 1 rule is given by E; = (Cot&Co) +
+
e(n(d-1) dn(a)-l)
I
c c
8=2
t=1
c
712
s,i+a+,t(t,8,r)(CV,+(Ho- Eo)@ca) (20)
The brackets indicate that the trace of the enclosed matrix must be evaluated. The 2N and 2N + 1 expressions are as follows. n ( a ) = 2N: E; = (Cot&Co) e(N) dZN-I)e(N-I)
c c c
t=l
n(a) = 2N
8'2
+ 1:
?=I
+
s,,~+a+,t(c,8,r)(C~+(Ho -
Eo)W(21)
E; = e ( N ) d2N) r ( N )
(50) Parish, C. A.; Augspurger, J. D.; Dykstra, C. E. J . Phys. Chem., submitted. (51) Adamowicz, L.; Bartlett, R. J. J . Chem. Phys. 1986, 84, 4988.
v=l 0.007 -15.19 -1 1.33 14.36 2.73 270.81 -179.36 -29.30 410.72 -361.55 -104.13 -41.62 -815.45 38.81 1777.06 -44.99 925.46 -29.38 1855.60
0.017 -15.00 -1 1.28 13.62 2.92 27 1.03 -169.56 -22.70 410.79 -339.55 -89.37 -41.29 -816.44 38.42 1779.63 -45.03 867.93 -29.49 1667.50
(CotHgCo) + C C €=I
8-2
7'1
6,L+a+~r(t,B,r)(C','(Ho- E0)Q) (22)
9236 The Journal of Physical Chemistry, Vol. 95, No. 23, 1991 Hartree-Fock Wave Functions. Derivatives of self-consistent field or Hartree-Fock wave functions given in a basis of orbital functions are more complicated because of the implicit basis dependence of the Fock operator. Beginning with the HartreeFock-Roothan equation FC - SCE = 0
(23)
it can be differentiated just as the Schrdinger equation. The first-derivative equation is
F"C + FC" - SOCE - SCOE - SCE" = 0
(24)
To solve this derivative equation, again the embedded parameters must be set to their reference values. The derivative orbital matrices can be represented in terms of the reference orbitals
cg = coug
(25)
and it is useful to introduce two new matrices
G" = CotF{Co and R{ = C$SgCo (26) so that, after multiplying the derivative equation by Cot and integrating, eq 24 becomes (since Ro = 1, and Go = Eo)
&Ut - U{Eo = R{Eo - Gg + E{
(27)
From this point, the subscripts will be suppressed, for in all cases we will be referring to solutions of the equations where the parameters are set to their zero or reference values. For a general derivative, a , eq 27 becomes
EU* - U ~ E = RUE - G* + E" +
ra
(28)
where I'" is completely determined by lower order derivative U matrices. e(n(a)-l)
c(n(a)-l)
=
6a,8+4+it(P,~,t)RBU'E'-
8,?~2
(29) Since the off-diagonal blocks of E are identically zero, eq 28 is used to find the off-diagonal-block elements of the U matrix by rewriting it for a specific element, where m and r refer to orbitals of different shells (doubly occupied, singly occupied, or empty)
=~
: ~ t- ,
The terms which depend on lower order U matrices, plus the diagonal blocks of U", can be assembled to form a fixed part of the derivative matrix, which in turn can be used as a first guess in assembling the derivative Fock operator to begin the iterative process for finding the off-diagonal U matrix elements. Those elements can then be added to the fixed derivative density matrix and the iterative process continued until the off-diagonal elements of Ua become constant. Finding the energy derivatives from the Hartree-Fock wave function is a separate step, since the eigenvalues of the Fock matrix are just the one-electron orbital energies. As in perturbation theory, the (2N 1)- and 2N-order energy derivatives can be found from just the N-order derivatives of the wave function. Using a 2N + 1 or 2 N evaluation is advantageous for two reasons. First, finding the U matrices is the computationally slow step, and so reducing the number of them required should reduce computation time. Second, whatever slight residual error that remains from finding the U matrices will be propagated, so that either more accurate energy derivatives can be found, or a lower convergence limit can be used for finding the U matrices to get equivalent accuracy in the energy derivatives. For the case of a closed-shell molecule, the Hartree-Fock energy is given by E = (D(h + F ) ) (35)
+
It must be differentiated to find the energy derivative.
The Fock operator consists of two pieces, the one-electron hamiltonian, h, and the two-electron operator, 2J(D) - K(D), which will be designated as Q(D). The derivative Fock operator then becomes
C &,,+,t(P,r)GQJ'
i3,7=2
u:,(t, - t,)
Augspurger and Dykstra
cp + r;,
(30)
where e, and t, are orbital energies. Equation 30 represents an iterative process, since the derivative Fock operator, G", depends on the U" matrix being found. Alternatively, this equation can be solved by writing it as a set of coupled linear equation^.^,^^ The diagonal blocks of U" can be found from the orthogonality condition.
U'RU = 1
(31)
Thus, the arbitrary derivative of the orthogonality condition is
Fa = h"
a
+ (Q(D))" = ha +
(37)
&,,j+&3,-y)@(DY) &7=l
where
[QB(D')I,, = CDZ,,[2(iiVk)O - ( i ~ l r j ) 8 1 k,l
(38)
After substituting for the derivative Fock operator in eq 36, and taking advantage of the identity (W@(DY)) = (D?@(D")), eq 36 can be written as e('v)
5
E" = 2 C
P I y=b(N+I)
6G3~+,t(P,r)l(D8ht) + (DYFB))+ e(M
a
C C ba,8+.i.+~t(P,r,t)(wQ(D7)) (39) B.-(=l t = 2 where n(a) = 2 N + 1. While the first and last terms of the summation only depend on U matrices of order N a n d lower, the second term depends on all orders. It is this second term which requires our attention. It can be shown that the dependence on N 1 and higher order U matrices can be eliminated, thus resulting in the 2 N + 1 energy evaluation. The first step is to replace the derivative density matrix, via eq 34, in the second term in eq 39.
+
and the individual elements of the diagonal blocks can be found by choosing these blocks to be Hermitian. Of course, this does not determine U" uniquely, since any antiHermitian matrix can be added and this condition will still be maintained. The important idea is that the diagonal blocks can thus be specified completely by just the lower order U matrices. The density matrix, D, depends on the U matrices and can be written as UTUt, where T is the zeroth-order idempotent density matrix represented in the molecular orbital basis. This can be differentiated as well.
e(.%)
X
8.7=1
&-j+.i.t(P,~)U@Ua'
(34)
~ ~ , ~ + + ( P , Y ) ( D=' F ~ ) e(M
a
c c
@=I y,e=1
6 a . ~ + T + : t ( ~ , ~(TU'+F~~UY ,t) (40)
The order of the matrices has been cyclically permuted, which can be freely done since the trace is invariant to cyclical permutation. Starting with differentiating the orthogonality condition (eq 31), it can be shown that e(n3
a
p=1
y,e=1
xc
a
D8 =
01
X X 8=I y = b ( N + I )
so that
~,,B+-;+~~(P,~,~)(TU~+(RU)YEB)
=
o
(41)
The Journal of Physical Chemistry, Vol. 95, No. 23, 1991 9237
Electromagnetic Properties of Molecules
with those of the third summation and rearranging, the 2N energy derivative can be cast in another form. e(N)
a
e(N)
(I
C C baB+~+it(P,Y,t)(TUet[FBU’- (RU)”’E@])(42) @=I y.t=l
E” = C
C
@=I y=b(N)
e(M
& , 8 + 4 P 9 y , l ) ( m 7 ) ( 2 - 6,(@,+,v)+ dN-1)
e(N)
E C C @=2 y=b(N+I-n(@)) e52
The key step now, as shown by Nee, Parr, and Bartlett,” is to partition these summations.
dN) dN) a
(RU)’E@I) ( 2 - bn(B),N) + C C C 6aj+T+zt(B,Y,t)
e(M a
cc-
@ = I 1‘’ t=2 e(N-1) dN-l-n(iW d 2 N - n ( @ + ~ ) ) dN)
B=l ?.#=I
(J-WYD’)) - 2
@=I
7-1
8-2 y=MN+I-n(@))
t=l
@=I
t=l
y=b(N+I)
The second term depends on only the N a n d lower U matrices, and so now the first and third terms must be dealt with. The first term can be shown to be identically zero by separating the = 1 term, and then substituting for the resultant F U 7 - (RU)’E by differentiating the Hartree-Fock equation with respect to y. The third term can be dealt with by relating it to the first term. Since F, R, and E are Hermitian, and E is block diagonal it can be shown that
(TUtt[FBUy - RWyE@]) (TU’+[F@U‘- R O E @ ] ) *
(44)
After expanding the (RU)Y, it can then be shown that the third term is equal to the first, except for terms which depend only on N a n d lower U matrices. So the energy derivative can be written as dN)
C
N
y=b(N+I)
@=I
&a+?t(P,Y,l)(m’)
+
C
C
C
@=I
C 6aJ+r+S+z x
(Pb(N-n(B+T)) e = l
ye1
e(N)
c
@-b(N)
e(N-1) @V-~(@+Y))
C
7’1
C
{=b(N-n(@+y))
I
C 6,*J+?+?+i‘(P,Y,S;4 x t=l (TU‘~RW~E@) (46b)
These two forms represent the “four-term” and “two-term” second-derivative formulas described by King and Kormornicki.21 They showed formally that the four-term formula (represented here in eq 46a) is more accurate, since errors in the wave function derivative enter quadratically, while the two-term formula (eq 46b) depends linearly on the error in the wave function derivative. Numerical Wave Functiorrs for Vibrational States. Vibrational problems with one degree of freedom can be solved numerically via the Numerov-Cooley method,52 where the Schradinger equation is converted to a difference equation, and the wave function is represented as a set of amplitudes. The derivative approach presented here has been applied to these types of problems as well by the derivative Numerov-Cooley (DNC) meth~d.~ Both the Numerov-Cooley and DNC methods can be outlined concisely. The Schradinger equation for a one-dimensional oscillator is 1 P”)(R) = [ U ( R ) - E ] P ( R ) 2m
dN-n(B)) d2N+I-n(@+t))
e(N)
2c c &y=l
til
c
f=b(N+I-n(@+t))
n
6a.a+?+i+F(Blr.c, ( T U ‘ + R W W
(45) When the basis is independent of any of the parameters, Le., Ra and Q“ are zero for all (Y except I , then all the terms in the last two summations are identically zero (for example, electromagnetic field derivatives when field-independent basis functions are employed). An analogous 2N energy derivative evaluation (Le. n(a) = 2 N ) is e(N-l)
E” = 2
C
C
-pb(N+I)
@=I
& d + ? t ( P , Y * l ) ( m y )+
el M
(DBQ’(D9)
-2
e(N-I) e(N) e(2N-n(@+y)) e(N-I-n(@)) ,911
I3
y = l f=b(N-n(8+7))
c
e-l
&i$+?+?+Z
x
f(P,yJ,c) (TUttRWyE@)(46a)
However, certain ambiguities exist in determining the 2N energy derivative expression which are not present for the 2N + 1 case. For the 2N + 1 derivatives, rearranging to arrive at eq 39 leaves all terms like depending on derivatives of the wave function beyond the Nth order. For the 2N derivatives though, the equivalent expression to eq 39 contains these types of terms where n(@)= n(Y) = N . In arriving at eq 46a, these terms were explicitly left in the final expression. However, by combining these terms
(m)
X
t(@,y,f,t) (TU“RWyE@)-
ell
(43)
Ea = 2
a,i6a,a+?+7f(P,y,t)(TU~t[FBUy -
(47)
where R is the displacement coordinate, P ( R ) the wave function, U ( R ) the potential, E the vibrational state energy, and f l 2 ) ( R ) the second derivative of the wave function with respect to R. Of course, there is a corresponding P(R) and E for every vibrational state, but for notational simplicity any subscripts denoting a particular vibrational state will be suppressed. U ( R ) consists of any interacting potential and Uo(R),the stretching potential. The vibrational coordinate, R , is divided into a set of M discrete points, R,, i = 1, M , of a small, common length. Then, the second derivative is approximated by the second difference. Making an initial guess to the vibrational state energy, E, it is then possible to find Ri-l from R , and Ri+, (or Ri+l from Ri and R,-l). Then, starting from the nonclassical regions, far out, and assuming that the farthest out point has a zero amplitude and the next point in has a very small amplitude, it is possible to find the amplitudes a t each point recursively. This process is also started from the near in nonclassical region until they meet. If E is the correct energy, the amplitudes will match. If not, a new guess to E is made, and the process repeated until E is found. The DNC method involves differentiating the vibrational Schradinger equation, eq 47, which for an arbitrary derivative a (n(a) = N + 1) is
1 --P”)“(R) 2m
= [ U ( R ) - E ] P ( R ) - EaP(R)
+ Xa(R)
(48)
where X”(R) depends on only N and lower order wave function derivatives. X“(R) = Ua(R) P ( R ) +
e(M &Y-2
6ag+lt(P,y)[I/s(R)- E @ l f l ( R )
(49) ( 5 2 ) Cooley, J. W.Math. Compur. 1961, I S , 363.
J . Phys. Chem. 1991, 95,9238-9242
9238
If eq 48 is multiplied by P*(R)on the left and integrated, and if the first term on the right is brought to the left, they will add to be zero, since P ( R ) and E are the wave function and energy of the zero-order Hamiltonian. Thus, the energy derivative is
Ea = (P(R)IX"(R))
(50)
(provided that the wave function is normalized), which represents an N 1 energy derivative evaluation. This can be substituted
+
back into eq 48 to then find P ( R ) , and since all other quantities are known, unlike finding the zeroth-order wave function, this process is not iterative. Since the derivative wave functions are found noniteratively, there is no advantage in computation speed or accuracy in developing a 2N/2N 1 derivative evaluation. Since this process consists almost entirely of vector adds and dot products, it is extremely well suited to take advantage of specialized vector processing hardware.
+
Ab Initio Study of the Nitronyl and Imino Nitroxides. Relation between Electronic Structure and Magnetic Properties in Metal-Nitroxide Complexes Andre Grand, Paul Rey, Robert %bra,*,+ Vincenzo Barone,$ and Camilla Minichinos Laboratoires de Chimie (U.A., C.N.R.S. I 194), Dgpartement de Recherche Fondamentale, Centre d'Etudes Nuclaires, 38041 -erenoble Cgdex, France, Istituto Chimico, Universith di Napoli, 4, Via Mezzocannone, 80134- Napoli. Italy, and Dipartimento di Chimica, Universith della Basilicata, 85100- Potenza, Italy (Received: February 28, 1991)
A complete ab initio optimization of the simplest nitronyl and imino nitroxides was performed at the UHF-SCF level, using a 3-21G basis set. The molecular electrostatic potential maps were drawn and correctly predict the coordinating ability of oxygen and nitrogen in both radicals. The largest fraction of the unpaired electron is always found on the oxygens, but the imino nitrogen bears a noticeable part of the spin population in agreement with experimental findings. All the trends of experimental EPR spectra are well reproduced by the calculated spin-projected densities at nuclei, the imino nitrogen being found to have the largest coupling constant at the SCF level.
Among the nitroxides, the nitronyl and imino nitroxides'-* (Figure 1) exhibit interesting properties in relation with their delocalized electronic structure. These stable free radicals have two sites of coordination and they have been extensively used for designing high-spin molecular based species3 Indeed, they can function as bridging ligands toward metal ions and afford extended arrays of alternating exchanged coupled organic and metal spins. We recently reported such unidimentional ferri- or ferromagnetic compounds using nitronyl n i t r o ~ i d e s . ~ , ~ On the other hand, the coordination chemistry of imino nitroxides has not been so much studied, although they exhibit similar delocalized structure and two sites of coordination. A preliminary investigation of the magnetic properties of some dirhodium derivatives of these paramagnetic species brought unusual r e s u h 6 While the oxygen-bonded nitronyl nitroxides exhibit large intramolecular nitroxide-nitroxide antiferromagnetic interactions mediated by the Rh-Rh bond, the nitrogen-bonded imino nitroxide complexes are only weakly exchange coupled. It was suggested that these different behaviors were related to structural differences: whereas in the 0-bonded derivatives an overlap between the a* magnetic orbital of the nitroxides and the u* of the dirhodium fragment is allowed, the bonding parameters of the imino nitrogen atom are such that these orbitals are orthogonal in the N-bonded adducts. Such a behavior would also have been observed if the unpaired electron was mainly localized on the uncoordinated NO group. However, extended Huckel calculations showed that the unpaired electron is actually delocalized on the imino coordination site in agreement with EPR results. Moreover, recent studies of copper(l1) imino nitroxide complexes' showed that the large observed metal-nitroxide ferromagnetic interaction is only compatible with a significant unpaired spin population on the coordinated imino nitrogen. Whereas the unpaired spin in the symmetrical nitronyl nitroxides is obviously equally distributed on the two N O groups, the situation is more complicated in the imino nitroxides in which the 'Centre d'Etudes Nuclaires * Universiti di Napoli. 8 Universitl della Basilicata.
two sites of coordination are different. Solution EPR spectra of the former exhibit a five-line pattern in agreement with two equivalent nitrogen coupling constants;' on the other hand, in the latter two differently coupled nitrogen atoms afford a seven-line spectrt" Concerning the magnetic properties of the metal complexes, EPR results are of no use since these magnetic properties are dependent on the overlap of the magnetic orbitals which may have a node on the nucleus. This is the actual situation in the nitroxides where the unpaired electron residues primarily in an orbital of a* symmetry. Therefore, the understanding of the magnetic interactions in metal-nitroxide complexes must rely on a more precise description of the electronic structure of both the metal ion and the free-radical ligand. Thus, we decided to perform ab initio calculations on the imino and the nitronyl nitroxides in order to get a deeper insight on the factors which govern the magnetic properties of the metal complexes of these free radicals. We describe hereunder the results of these calculations which show that, in the imino nitroxides, the unpaired electron is indeed substantially delocalized on the imino atom and that the old attribution of the larger nitrogen hyperfine splitting to the NO group is erroneous. In addition, we report a full investigation of the EPR coupling constants in both types of radicals. Theoretical Background
The unrestricted Hartree-Fock (UHF) formulation provides a method for computing magnetic hyperfine coupling constants in free radicals which is generally considered as a useful alternative ( 1 ) Ullman, E. F.; Osiecki, J. H.; Boocock, D. G. B.; Darcy, R. J. Am. Chem. SOC.1972, 94, 7049.
(2) Ullman, E. F.; Call, L.; Osiecki, J. H. J . Org. Chem. 1970, 35, 3623. (3) Caneschi, A.; Gatteschi, D.; Rey, P. Prog. Inorg. Chem., in press.
(4) Caneschi, A.; Gatteschi, D.; Renard, J.-P.; Rey, P.;Sessoli, R. Inorg. Chem. 1989, 428, 1976. ( 5 ) Caneschi, A.; Gatteschi, D.; Laugier, J.; Rey, P. J. Am. Chem. Soc. 1987, 109, 2191. (6) Come. - A.; Grand, A.; Rey, P.; Subra, R. J. Am. Chem. Soc. 1989, I l l ,
3230. (7),Luneau, D.; Rey, P.; Laugier, J.; Fries, P.; Caneschi, A.; Gatteschi, D.; Sessoli, R. J . Am. Chem. SOC.,in press.
0022-3654/9 1/2095-9238$02.50/0
0 199 1 American Chemical Society
