- Email: [email protected]
Gauge-dependent ket calculation of the magnetic properties of molecules
Journal of Molecular Structure (Theochem), 161 (1987) 39-60 Elsevier Science Publishers B.V., Amsterdam - Printed in The Netherlands
GAUGE-DEPENDENT RET CALCULATION OF THE MAGNETIC PROPERTIES OF MOLECULES Part I. Method and application to Hz, Hz0 and CO
J. P. FLAMENT, Groupe D.C.M.T.,
H. P. GERVAIS and M. RERAT Departement
de Chimie, Ecole Polytechnique,
91128 Palaiseau (France)
(Received 15 September 1986)
ABSTRACT A gauge-invariant ket function theory of magnetic properties is presented. The perturbed state, which depends explicitly on the strength of the external fields, is determined variationally as a linear combination of Slater’s determinants associated with a gauge transformation factor. This “mixed” procedure combines advantages of both the polynomial and the familiar SCF-CI methods; it is applied, with inclusion of electronic correlation, to H,, H,O and CO. INTRODUCTION
Considerable interest has been devoted to the quantum mechanical calculations of the magnetic properties of molecules. The magnetic field does not appear in the theoretical ex ressions through its strength ‘B, but through its associated vector potential i for which one has to choose a particular gauge origin, i.e., a position vector Luch that A = l/23 A (i’- A). In the case of a molecule, there is no natural choice; use of infinite basis sets should lead to gauge-invariant results, especially with respect to the origin; in practice, sufficiently large basis sets that make the gauge problem irrelevant are hardly attainable. To solve the problem of translational invariance, an interesting approach uses expansions of the molecular orbitals over a set of gaugeinvariant atomic orbitals (GIAO) [l] , but it needs calculation of additional two-electron integrals and, to our knowledge, correlation effects have never been considered; the IGLO method [2] has been proposed to overcome these difficulties. When correlation effects are considered, GIAO are not used and large basis sets are necessary; gauge invariance however is not obtained [3a, 4-61. In the GIAO or IGLO methods gauge invariance is guaranteed; both in fact involve a “constrained” gauge, in the sense that the origin of a correcting factor is set either at the nucleus of the various atomic orbitals (GIAO) [l] or at the centroid of charge of localized orbitals (IGLO) [2]. Some other methods optimize the position of the origin by minimizing the paramagnetic contribution, eventually with inclusion of correlation effects [ 7-91. 0166-1280/87/$03.50
0 1987 Elsevier Science Publishers B.V.
40
It has been shown that constrained gauge invariance through GIAO’s does not ensure current conservation [lob] ; in the same paper, it has been proposed to consider the gauge function itself as a variational parameter, as already suggested and applied by others [ 11-131. This last procedure usually leads to acceptable results for proton magnetic shielding, but fails in the case of heavier nuclei [ 141: one major reason for this is the difficulty to find a suitable gauge function which accounts for the existence of nodal surfaces in the wavefunction when heavy nuclei are present. In contrast, customary perturbation methods would be well-adapted to this situation, but they unfortunately do not ensure gauge invariance (which could only result in a complete basis). The computational procedure we propose here follows Epstein’s suggestion [ lOa] : to translate both the “geometrical” effect of a gauge transformation and the proper physical influence of the field B,,, one has to consider the gauge factor 0 as well as the ket 19) as variational parameters in the expression I$) = e’“l$)
(1)
of the perturbed
molecular
state vector.
METHOD
Unitary transformation of the states Let 10) be the normalized ground state of the system in the absence of any magnetic field (10) belongs to the orthonormal set { 1~))formed by the eigenstates of the field-free hamiltonian A,,). The gauge transformation operator 0’ takes the familiar form
where 6’ (resp. &) is hermitian [I51
and a function
of the electronic
coordinates
We may consider the ket I$), appearing in eqn. (1) and depending on-the vector potential 6, as resulting from the action of a unitary operator ezs on the field-free ket IO) [16] I\ir) = p(e)
= eieee’glO)
For the field-including E(a)
= (*IQ*)
(41 hamiltonian
A, this leads to an energy
= (Ole-‘S^e-‘BAe’6eis^10)
In the following,
we shall be concerned
(51 by expressions
of the form
=
C
a,y:
(6b)
P
By hermiticity of @, we have a,, = a;. The development of the hermitian operator, $ = q/a li Sl, IlXmI can be limited to the part operating on IO> s =;
=;
c (S,lm>(Ol + S~lO>W) m#O
c c%nfik+ WL)
(7)
mf0
It acts as mixing the states {Ii>} under the influence of the external state field. In the magnetic case, the perturbation operator contains an imaginary part, so that the S matrix is not reducible to its imaginary component, as in
[W. Variation of the energy Using the development [17] : e&fie’ = & - i[a, 81 + 1/2[a, [g, A] ] + . . . one gets from eqn. (5) for the total energy associated with the ket I$> E(A, 6, $) = (Ole-iS^e-ibAeiQei~lO) = + (Ol[&,
[A, &:I ] IO> + 2(OI[& [A, 611 IO>}+ . . . (8)
that is, in matrix form, by use of eqns. (6) and (7)
E=Eo--i(eP,S,,,,S$)~~+
1(a~,S~,S~)~~~~~~~~~~)
(9)
with E. =
(01 [fp, A] IO>
VI, = ; (01 [I?;,
(Ad,,
=.$
(A&m
P%
(01 [?p,
?“I 1 10)
=$OI[[?,,,
A],eJl0~
A] IO>
v,,
=;
(A13),m
(01 [II,,
=$
W[[?j,,
Al, %ll IO)
A] IO>
42
&dm,,
= $
(&dmn
(01 [fik,
(01 &n,
(&A,n
[A, R’,]] IO)
rmJl10) (&~hfi = $ (01[km [A,
=$
(A32hn
=$
colr&l,
(01 P%,
[A, R, ] ] IO) (A33L
[A, R; ] ] IO)
$11 10)
=$
=$
COl[ll,,
[A, I?“]] IO>
The A, matrices satisfy the following relations A22
V, = -v;k=
Vand
=&=fj=()
Ai2 = A,, = A,*, = A& = C
i A2,=Af2=D from which, putting AlI = A and S, = X, + iY, (X,,
Y, real), we get
E = E,, .- i(a, S,, Sg)
E = E0 - 2(a,, X,,
=
Y,)
E,, - 2iiv + i.iMu
with obvious definitions for U, v and M. The stationarity conditions for E with respect to a,,, X, and Y, is thus Mu=v
(11)
Order analysis of the matrices Since we are interested in E due to the changes induced into # by the presence of an external field, we now develop eqn. (11) in powers of the perturbations and solve it for each order separately. For this purpose, every essential quantity will be expanded as follows
v=v,+ v,+ a=a,+
v2
al+ a2
43
x=x,+x,+x* Y= Y,+
y,+
Yz
As already mentioned, 10) represents the ground state of the field-free hamiltonian A,-, and thus belongs to order 0. Furthermore, as in other gaugeinvariant calculations [l, 2, lo], we restrict the gauge operator &’ to be a linear function of the positions (?j) of the electrons, i.e., with eqns. (3) and (6) e’=
c al,$I,withZ+
m+ n=l
(12)
lmn
We now concentrate upon the magnetic case (but the forthcoming developments may be easily adapted to the case of an external electric field). Zero th-order MouO = u.
(13)
Using the common commutation rules, we have w; = ; (01EM, A01IO)= -iq(OI($
lgrr
+ar,*$)lo>
For 5, restricted by eqn. (12), 8’: reduces to WE= -Siq = 0 for a closed-shell ground state.
When 10) and Im>are exact solutions of the field-free Schrijdinger equation, V&,, disappears. In current calculations, however, where these kets are only approximations, this normally no longer holds; fortunately, if IO>approximates sufficiently well the exact ground state, V&,, should be small, so that we shall neglect this term in the following. It thus follows that u. = 0 and u. = 0. First-order The first-order equation is Maul +Mluo = ul. On the basis of the preceding approximation, this reduces to M,u, = u1 where
(14)
CL
i (P
=COl&lO>) ,
Equation (14) then factorizes into Y1 = 0 and
Up to now, no hypothesis has been made about the nature of the vector potential A; the foregoing development thus applies to any kind of static field. In the general case, d may be a sum (e.g., of an external magnetic field B. and of the field due to a nuclear spin). As long as the fields are independent, eqn. (14) can be split into as many equations. It should be noted, however, that a “coupling” appears between the gauge and the field through matrix CL so that we may split the gauge contribution a, into parts related to the relevant fields. This means that, when considering chemical shift, the gauge function G’ does not depend solely upon the external magnetic field (a line already suggested by Hameka [ 111). Second-order We now focus on magnetic susceptibility appear at second order. Let us define
and magnetic shielding, which
and
(17b) where 8 is a homogeneous magnetic field, 3 the magnetic moment of nucleus N centered at &, q the position of electron j about some fixed origin, F”j = i$ - Rt, the position of electron j relative to N. Substitution of
45
into the equations leads to the calculation of x and u. With u solution of eqn. (ll), and recalling that u. = u. = 0, we get to second order E = COlA,lO>+ (OIA1lO) +
since fi, is purely imaginary
and
E =
Et')=(01C &?lO>+ (01C x'210, + 2(01 C & i - (i+‘l
i
l
$‘lO, -i&u;--i?al’u;
i
+ ii:‘u\)
With the use of eqn. (14), we may now demonstrate
the exchange
c’lv;’ = ti:‘u; since we have (the first-order
theorem (19)
eqn. (14) is linear in the perturbation)
M,u; = u;
PO@
M,u'l= u;'
POb)
Multiplying eqn. (20a) on the left by ii;’ and eqn. (20b) by ii:, and using the fact that MO is symmetric, we get eqn. (19), which is an extension of the usual exchange theorem, including the gauge. to Vectors u’, and ~9, as well as vi and uy, are to first order proportional IIEfll = B or II311 = M. We thus write a’=+e~
(21)
where Pis a unit vector along g or $. With this notation we get from eqns. (20a and b) a new set of equations where the coefficients are the same for both cases with eqns. (20a) and (20b) differing only by the type of operator (ds or la,) appearing on the right hand side
Using the {(I,,, x,} solutions
lG2) =
eg-;
of eqn. (22), the second-order
ii’u’ +
*
energy is
-$-* mffi’u” + . , . (23)
=
-~r3(xd+ f)S
+ B(od+ CJP)~ + . ..
where suffixes d and p denote as usual the diamagnetic and paramagnetic parts of the susceptibility and screening tensors. Using the definitions of eqn. (21) for ds and A,, it follows that
(244 1 1 -Zx~=-.-$&=--,
1 P
ap WkBfi+ C Xe VzP m
(24b)
=--
(24~)
[q A (Gj7p)]plO) + 2x m
xf$tOlC(3
A
(24d)
j
WW
=-g(cCLaMP(OI w 1 (I$ A (Gj,))&lO) i
+ 21 X~‘(OIC(~ m j
A Gj),lm,> (25e)
Equations (25d+e) clearly show the opportunity given by the exchange theorem. In the former, eqn. (25d), ap and X”,” are obtained when solving
eqns. (14) or (20a) for the perturbation BO; in the latter, eqn. (25e) afP and Xzfi are obtained when solving eqns. (14) or (20b) for the perturbation Mk Knowledge of (allBar, Xf) alone leads to both x and (I.
Connection with other methods By making specific choices of G, S and IO>,the above procedure reduces to some previously proposed methods: (1) Setting G = 0 and IO>= IHfl, the states Im>appearing in S are simply the monoexcited determinants ii + p). Vh = 0 then is no longer an approximation, since it expresses Brillouin’s theorem. This situation is closely related to the Tamm-Dancoff approximation (TDA). The equations are readily written
C [6,,Sii(e,-eEi)+ ia
2(ipljq)-(iqIjp)]X,j=--V&
which is an approximation to the CHF equation [ 181 E
jq
[6p*6ij(fp--i)
+ (WA--(0
IPq)lxqj=-v&
obtained with &, = hii = UiUi instead of Fiji = e$ilHF)(HF(. The CHF procedure is related to the random phase approximation (RPA) in which the B’matrix is not null. The decoupling has been introduced when eliminating the terms Im>(Zl with I # 0 in S, an operation which suppre_sses_ter_ms like S om X S,, in the development of the double commutator [S, [H, S] ] (although they do not act on IO), they introduce a coupling between the operators Im>(01). (2) Setting G = 0, with IO>being a correlated wavefunction, one gets Fukui’s SCF-CI method [3a] (see also Svensen [3b] for the case of electric polarizability). (3) Setting 10) = IHF), S = 0, and varying G alone for some analytic forms leads to various developments of gauge transformations [ 12-14, 19-221. (4) Finally, if we expand the exponentials to first order we get a wavefunction of the type I@‘) = i(&lO> + &$,In>). This is closely related to Okninsky-Sadlej’s method [ 231 where 19) refer to orbitals, while in the present paper it symbolizes states or Slater’s determinants).
Practical implementation The correlated ground state IO) has been determined by means of the CIPSI algorithm [24], which builds the wavefunction by iterative selection of the most important determinants; the final step included all those which had a coefficient superior to ca. 0.01 in the M$llerPlesset first order wavefunction. Use of states Im>in S would require a lot of computation, so we generated them in the spirit of CIPSI and perturbation theory: once IO>has
48
been built as a combination of Slater determinants (mainly 1HF) and diexcitations), all mono-substitutions are introduced to generate all mono-, di- and tri-excited determinants orthogonal to 10) (i.e., that do not2ppear in its development) which interact with it b+y means of operators L (coefficient sB or sL) = l(OlLlm)/& - ,?&I) or LN/r& (coefficient sM (or S&$) = I(01k N/r&lm>/Eo - &I) (both criteria can be used indifferently because of the exchange theorem). We then make the approximation (Olfi,,lm) = 0, which is not a serious drawback, since the most important determinants of 10) have already been taken into account (moreover, in highly symmetrical molecules this matrix element is generally zero for symmetry reasons). y,, being limited to monomials of degree 1 (cf. eqn. (12); y = 3c,y, z) coefficients appearing in the two parts of eqn. (22) write in the case of perturbation B, (n, is the number of electrons): 0
n,
2 0
-m$o
.
-&IO)
. ..=
0
-(ml
0
-(mldy d IO>
. ..=
-_(ml$lO)
. ..=
%
(mlHO-E”ln) -m-go
0 -
cd 2 (Y) T
(26)
. . . = -my$--s&O .
The closed shell ground state IO>belongs to the totally symmetric representation of the molecular point group; for the perturbation B,, the kets In> belong to the irreducible representation l’(L,). If the molecular point group is such that I’(&) differs from I’(x, V,), I’(y, Q,) and I’(z, VZ), the coefficients a, of the gauge function vanish (a,@) = a,,(%)= a?) = 0) when the origin is the electronic centroid ((x) = (y) = cZ>= 0). When this holds also for the perturbation B, and B, (as in the case of D-h), the problem reduces to its SS part and the gauge function is zero [25]. For less symmetrical systems, however, the gauge function must be taken into account: in the water molecule, for example, taking Oz as the CZ axis, only B, obeys the preceding rules if the origin lies on that axis and a GS coupling appears for the B, and B, components; when we move the origin away from the z=axis (e.g., from the oxygen to one hydrogen), new determinants will be generated by g,(H) - L,(O) = XOH a /i3y and gauge invariance may not be guaranteed if one considers only the SS part of eqns. (26). In contrast, if the GG and GS parts are explicitly taken into account, the gauge function absorbs completely the translation of the origin and the coefficients of the extradeterminants disappear; thus suppose that X, = XL (t/m) when
49
translating the origin, the G lines give a:--ua,=
0
a;-aa,=Zon
(27a)
a: - a, = -Yen 1 whereas for one S line -_(mlV,lO>(a: -a,)
- (mlV,IO)(a~ - uy) - hlV,I0>(u:
- Z,,Qrw,lO)
-a,)
= Y,,GnlV,I0> Wb)
Equation (27b) still holds after substitution of eqn. (27a). Consequently one may put the origin at the point of maximum symmetry, in order to reduce the number of kets lm>generated by perturbation. APPLICATIONS
Three procedures have been compared, in the case of the small molecules Hz, Hz0 and CO (at a fixed geometry, that is without introduction of rovibronic corrections) to the “polynomial” method (a f 9, X = 0) with degree p superior to one for the gauge function G, the “SCF-CI” method (a = 0, X # 0) and the “mixed” method outlined above (a # 0 with p = 1, X f 0). Hz (dnn = 1.401 u.a.) Two basis sets were used: (A) that of Lie and Clementi [26] ; (B) that of Iwai and Saika [4] (i.e. (A) augmented by (Is, 2p, Id) on each hydrogen). The ground state wavefunction was obtained through SCF (HONDO) and CI (CIPSI) calculations, as described in the Practical implementation section. With our actual program (diagonalization below dimension 200), a complete CI was possible for basis (A); for basis (B), more than half the CI space (S = 196) was diagonalized, the second-order M$ller-Plesset barycentric correction to the energy [24] due to the remainder being less than lOA u.a. As concerns the polynomial method (a f 0 with p = 2, X = 0), Fig. 1 shows the behavior of a and acX obtained, while increasing the dimension of CIPSI’s S-space. In the case of basis (B), one notices a sudden rise of energy at S = 196 (account being taken of the few mono-excited configurations which appear at this stage). Results from the various procedures are given (basis (B), S = 196) in Table 1. In the case of the polynomial method, gauge invariance is respected when translating the origin from the center of mass Sl to nucleus H. Thus AG, = (0.596yH - 0.149ynzu) = (0.596y - O.l49y(z = (0.596
-
0.149Z)y
- (-0.149yz)
+ 2) + 0.149yz =
+ 0.149 $
>
y
Fig. 1. Polynomial method (H,). Evolution of 6~ and or(H) dimension of space S.
(Basis sets A and B) versus
Development of & to the third degree is useless, since in Dmhgroup no monomial of this degree belongs to the same irreducible representation as fl,, &, L,. Development of & to degree one implies that the paramagnetic contribution to magnetic shielding is zero when the origin is put at the electronic centroid (a here). Comparing the SCF-CI and the mixed methods (the threshold & or sL for selection of the excited states lm> being identical and small, s < 10e4 u-a.: nearly all kets lm> are introduced), one sees that both give exactly the same result when gauge origin is set in a. Gauge invariance still holds for the mixed method (the results depending, of course, on the basis chosen), but not for the SCF-CI, even with a large basis and CI calculation. Finally, values ii, = 26.64 ppm and u%“(H) = -8.61 ppm, as given by the mixed method, can be compared to other values mentioned in [ 271 (although a rovibrational correction would be necessary for comparison with experiment). By contrast, the polynomial method, even if ensuring gauge invariance, is far from likely to give good results; this conclusion will be more explicit below, where heavy atoms are present. Hz0 (doH = 0.958 A; HbH = 104.5”) Satisfactory values of magnetic properties need to be calculated with large basis sets. We thus choose basis (IV) employed by P. Lazzeretti et al. [28]
51 TABLE 1 Hydrogen molecule (A) Polynomial method Monomial basis (degree < 2)
Gauge origin: H
Gauge origin: n a
G,
G,
0 0.596 0 8 0 4.149 0
34.63 0 0 -9.04 0 0 25.59 0
“H
2i.35
GY
0 0 0 0
-0.596 0
0”
z
3: 63 0’ 0 -9.04 0
-:
0.
00 8 0.149
149
0 0
25.19 0 0
0 0
27.88
0 25.19 0
0.40
8 0
2: 59 0’
0 0
8
0.149 0 0
840 0’
: 25.59 0 0 26.35
8 27.88
GZ
GY
2: 59 0’
0 0
27.88 0 0 0 0 0 27.88
For a polynomial of degree 1, &@(sL) = 0 v (a, p) 4 6~ = 26.09 ppm For a polynomial of degree 3, GH = 26.35 ppm (B) SCF-CI method (threshold s = lo-+ u.a.) Determinants introduced by Lp uv = 2981146 a;” = 2981146 a? = o/o
Gauge origin: H -8.36 0 0
0 -8.36 0
OH
26.27 0 0
GH @pm)
26.81
(C) Mixed method (threshold s = lo+ gauge origin OH =
26.02 8
PXX CH (H) =
0 26.02 0
-8.61
0 26.27 0
u.a.; polynomial
Gauge origin: ~2 0 0 0 0 0 27.88
0.83 0 0
0 0.83 0
0 0 0
26.02 0 0
0 26.02 0
0 270.88
26.64 degree l), results independent of
0 0 27.88 >
ppm
&“(
n ) = 0.83 ppm (for this choice of gauge origin, all coefficients zero; the value of 0.83 ppm is purely SCF-CI) OH- 26.64 ppm
%enter of mass of the molecule.
in the polynomial m
52
and consisting in the [6s, 5p] contraction of the (lls, 7~) of Huzinaga and Sakai’ [29], increased by a 3d (STO-SG, 5 = 1.93) for the oxygen atom and [ 2s, 2p] A.O.‘s of Salez and Veillard [ 301 for the hydrogen atoms. Trying to improve gauge invariance in the case of the pure SCF-CI method, we added to this a set of four d shells (5 = 1.0) centered in the molecular plane, 0.5 u.a. away from the oxygen atom on the coordinate axes. This did not prove to bevery useful, but we choose first to present (in Tables 2 and 3) a comparison of the three methods via the results computed on this basis of 62 A.O. In Tables 2 and 3, the S space of CIPSI has a dimension unity (i.e., the diagonalized energy is that of the SCF calculation and is taken as zero energy, TABLE
2
H,G: ,,@ oru) uorp and Cr Compa%on Af results given by the SCF-CI 4d, dim. S = 1, E, = 0,threshold s = 0) (A) SCF-CI:
P
x
(d)
x
(p)
Y Y
Y
(B) SCF-CI:
0
0 0 0
-39.21 39.04 -0.17
0 0 0
101.05 -59.01 2.05
17.62 -7.35 10.26
415.21 -53.79 361.42
(d) (P)
00
(ppm)
0)
0 0 0 ; 0 325.85
z
X
0
21.14
0
- -11.19
0
9.95 0 0 0
8.38 12.82 21.19
0 0 0 417.18 -112.44 304.74 0 0 0
0 0 0 0 0 0 415.74 -104.35 311.38
-10.85
0
23.11 6.14 29.25
0 0 0
Y
44.80
basis +
39.08
2.97
-7.04
C(
1,1,l) z
1.04
0
38.19
40.69 (C) Polynomial
z
Y
(d)
s z t
(Lazzeretti’s
Gauge origin:
29.32
u,, (gauge origin:
(P) (d) (P)
31.67 5.83 37.51
44.06
X
0
Y
0 0 0
x
Y Y Y
X
130.32 -85.79 44.53
i
x
X
z -39.21 39.25 0.05
P
ct
Gauge origin:
0 0 0
(d)
(PPm)
H
Y
75.67 -30.08
(P)
z z z tin
x
45.59
X
methods
oH Gauge origin:
(Y
and polynomial
method
Degree of polynom.
OH @pm)
ho (ppm)
1 2 3
24.84 26.32 27.03
415.13 415.04 414.42
53 TABLE 3 H,O: c$$orn), 0~ and 6, results given by the mixed method (Lazzeretti’s basis t 4d, perturbation, dim. S = 1, E, = 0,threshold s = O)*
Gaugeorigin:
P
x
H
Gaugeorigin: 0 z
Y
x
Y
L,
Gauge origin:C(1.1.1) I
x
Y
z
b
OH 1 2 3 4 ~~(126
2.10 0 0 0.99 0 75.67 0 -37.32 det.) 38.35,0,8.92 (37.36,0. 13.35)
1 2 3 4 oyp(156 det.) 1 2 3 4 a&96
det.)
+(PPm) 00
0 9.36 -1.00 0 0 130.32 0 -108.38 0. 21.94.0 (0.14.06.0) -7.35 0 0 -1.43 -39.21 0 49.47 0 10.26,0. 29.25 (17.62,0. 23.11)
-5.95 0
-39.21 48.13
0 1.43 0 0
6.14 0 101.05 -71.80
29.85
0 2.10 -0.12 0 0 31.67 0 6.67 38.35.0.8.92 0 9.36 0.11 0 0 8.38 0 13.56 0, 21.94,0 0 -7.35 0 0 0 17.62 0 -7.35 10.26,0. 29.25
-5.95
0 21.44 -12.22
0
0 0 0
6.14 0 23.11 6.14
2.10 0 0 0.88 71.37 -54.45 -93.02 54.45 38.35,0. 8.92
0 1.00 0 0
-7.35 0 1.00 -1.00 -22.08 -39.70 32.34 39.70 10.26.0.29.25
6.14 0 77.56 -48.31
b
1 -47.12 2 0 3 404.59 4 -37.60 o,p(126 det.) 366.99,0,0 (413.84.0.0) 1 2 3 4 oYp(156 det.) 1 2 3 4 a,+(96 det.) oo(PPm)
0 0.99 0 0
0 0 0 0
0 -107.69 -1.00 0 0 406.56 0 -98.09 0,308.47,0 (0.415.81.0)
0 1.43 0 0
0 0 0 -1.43 13.71 0 -13.71 0 0. 0, 311.38 (0.0, 415.73)
-104.35 0 415.73 -104.35
328.85
-47.12 0 -0.12 0 415.21 0 -48.22 0 366.99, 0.0
0 0 0
0 -107.69 0.11 0 0 417.18 0 -108.70 0, 308.47. 0
0 0 0
0 0 0 0 0 0 0 0 0, 0, 311.38
-1.00 -33.31 42.23
0 9.36 -0.89 0 0 102.52 0 -80.58 0, 21.94.0
29.85
29.85
-5.95
0
0
-104.35 0 415.73 -104.35
328.85
*Values in parentheses are given under the same conditions by the polynomial method (degree 1). bl. SCF-CI paramagnetic contribution. 2. Coefficients of the polynomial. 3. Diamagnetic contribution. 4. Paramagnetic contribution. One can verify that a change 3 o f gauge origin induces variations in the polynomial coefficients, exactly equal to the corresponding coordinates of $, thus satisfying AG, = (3 A &. (The position of the hydrogens is (1.43; 0; ll.lll), in ua., with respect to the oxygen.)
54
E. = 0). In this case, the excited kets lm> generated by B, (resp. M,J are not too numerous for the resolution of the linear system eqn. (22), even if sBa = 0 (resp. sMa = 0). Gauge invariance (AG, = @A it),), as well as the exchange theorem, can then be verified in Table 3 for the mixed method. When the S space expands up to dimension 192 (E,-, then becomes the lowest eigenvalue of fi,, projected onto it) it appears that & generates too many kets lm> for a current resolution of the linear system of eqn. (22) if stcu = 0. If sLa # 0, however, gauge invariance is not automatically ensured, because perturbation L, depends on the origin, as well as the kets Irn) generated by it. To avoid this difficulty, we choose sM as a threshold for selecting Im),even in the computation of susceptibilities (sM then has to be as low as possible): for HzO, we calculate the magnetic shieldings cn of hydrogen and co of oxygen, after selecting excited kets Im>,respectively with the help of the values (denoted by xH and x0) of SLY/,+ and of SL~/'~, and then compare magnetic susceptibility obtained via these two basis sets. Table 4 reports the results obtained for various values of sM, in Lazzeretti’s original basis of 42 A.O. Figure 2 is a plot of cn, uo, xH and x0 versus colog s,. Convergence is relatively slow, which requires resolution of large systems TABLE 4 H,O: @(CH), @(Go) an d x‘@(x), results given by the mixed method (Lazzeretti’s basis, L&h perturbationa, dim. S = l-192, for various selection thresholds) x(10-"J T-2)
~Hbpm) S-l.s=0 36.37 0 0 21.77 10.35 0
8.93 0 29.23
29.79
-13.78 0 0
0 -14.01 0
0 0 -14.64
-14.38
;,=,‘,““9 ‘;H/& 21.38 0
0
0 311.36
-14.73 0 0
0 -16.62 0
q)(PPd 0
0 -15.59
-15.65
s = 192.q.J= 0.001,90 = 0.005b 38.31 0 9.57 -14.15 0 0 21.37 0 0 -14.35 11.49 0 28.79 0 0
29.51
0 0 -15.59
-15.65
29.49
0 309.01 0
~'(10-~J T-*)
xH(lO-" J T-*)
25.59
366.72 0 0 329.03
S= 192,s= IO((I~)}= O; polynomialonly)b 37.73 0 12.84 -14.73 0 0 15.08 0 0 -16.62 17.07 0 23.77 0 0
10.90
0 0 -14.02
-13.94
aH(PPm)
0
qJ(PPd
-14.19 0 0
413.88 0 0
0 415.66 0
0 0 415.71
0 303.61 0
0 0 312.01
0 311.17 0
0 0 323.02
415.08
0 -14.35 0
0 0
-14.39
-14.32
371.10 0 0 328.91
= o.oo5,sLo = O.OOlb 9.42 0 28.87
-14.13 0 0 -14.27
0 -14.30 0
0 0 -14.39
-14.11 0 0
0 -14.12 0
0
0
-14.28
-14.15
aExcept for the last case. bFor the definition of xH and x0, see text.
375.08 0 0 336.42
55
‘. ‘.
*.
“....................,...,*
341.55
CT c
UH
346.56
_.,’
g
%.,
337.40 336.42 .._,# . . . . . . . *n 328.91 ~..____........... .a , . 29.44 29.49 29.51
25%9
X
z :: ;;
-14,3&y,. -14.54 .y:.. / , ,-’
-Mll, ..
-15.65...-.“’ q. .
-15.33 .=
I 0
I -I
’ ..
-15:65 I I
-14.15 .m -1y.27 ..*-. -14.38
XT, . ‘,.i$O XH.
I 2
I 3
colog
5
Fig. 2. H,O (dim. S 0 192). Convergence of the calculations as function of the selection threshold s of the excited kets Im).
of linear equations; a good procedure for this, especially when a lot of matrix elements are zero, is the Gauss-Seidel method (where only non-zero elements are stored, thus needing less memqry and less I-O operations). More, over, the c-algorithm [31 j is well suited to this method and can be added to accelerate convergence of the iterative computation; in this way, systems of the order 10 000 containing sometimes 99% of zeros can be solved. Other calculation results, as in refs. If or 32, can be referred, as well as the experimental values mentioned in ref. If (an = 30.2 ppm; o. = 334 ppm) or refs. 33-35. We must emphasize that the value u. = 334 ppm is obtained from a private communication mentioned in ref. 36 giving the quantity (%O
-
aH,O)*
For H20, it appears that the S space (i.e., electronic correlation) does not contribute very much, so that one can be content with the results obtained for S = 1 and s = 0. co (dco = 1.13 A) For the oxygen atom, the basis chosen was that already used in the case of HzO; for the carbon atom, it was a [6s, 5p] contraction of the (115, 7~) basis used by Huzinaga and Sakai’ [29], increased by a d shell ({ = 0.7) (see ref. 38). With this set of 52 A.O.‘s, we obtain Es,, = -112.777876 hartree. In this demonstration calculation (comparison with experiment would need rovibronic correction and will be described in a further paper), dimension of the S space was limited to 128. As in the case of HzO, we had to choose sM instead of sB as a threshold for selecting the excited kets Im>, when dim S > 1; corresponding values of x therefore are denoted xc and x0 in Table 6.
56
Table 5 compares results given by the SCF-CI and mixed methods for dim 5’ = 1 and s = 0. As expected, both are close when gauge origin is put as the center of electronic charge G. Table 6 shows the influence of selection threshold s for dim S = 1 and 128 (see also Figs. 3 and 4). In contrast to the case of HzO, the S space (i.e., introduction of electronic correlation) deeply transforms the results, e.g., for the dipole moment, the magnetic shielding constants bc and Do are almost zero. Calculations with dim S 9 1 and s 2 0, implying resolution of large systems of linear equations, thus seem necessary in the present case. Comparison with experiment (By = -16 ppm, ~c = +3 ppm) [35] still appears unrealistic, because it requires introduction of rovibronic corrections, which can be important (see other theoretical calculations in ref. 37). CONCLUSION
The mixed method, as it was outlined and applied above, originates in Epstein’s suggestion [lob] to ensure current conservation together with gauge invariance. Perturbation due to the field is treated as in SCF-CI procedures (thus introducing electron correlation, in a manner similar to that of CIPSI algorithm), whereas a polynomial part “absorbs” the effect of the gauge. (as in Guy’s and Stephen’s procedures). In the SCF-CI aspect, the perturbed ket should normally have the form I$> = Z I, ,,,Slm In (m I; in the present treatment, we limited ourselves to 14) = Z,,,S,lm)(01 + c.c., thus neglecting terms like (ml&lO) (m f 0). The results show that this approximation has a negligible effect. Programmation aspects of the mixed method are relatively simple, because they only involve adjonction of three columns and three lines (for a polynomial of degree one, which is the most convenient case) to the system of TABLE
5
CO: oF(CC), @(CO), xw(x), comparison of results methods (dim. S = 1, s = 0) (U in ppm, x in J T-*)
given by the
SCF-CI and mixed
SCF-CZ
C!J(X = -93.14 DC OXX = -65.44 =C G”X = -77.47 UC
GE = 27.93 Czz = 271.28 OC
* Oc = 37.26
X xx = -14.20
x 1O-5
X *= = -18.09
x lo-$
u8
= 411.00
=+ Gc” = 46.80 52 = 38.78
Mj$ed method ax =-79.75 u”c”= 271.28
o”o” = -141.05
= og”
* x = -15.50
==+Trg = 42.97
x lo-”
57 TABLE 6 co: “?(Cc), c$j%o), Xoa(X), results given by the mixed method (L.&k influence of dim. S and of excitation threshold sa
perturbation),
XC
S = 1, 8 = 10 ( {Im)} = 0, polynomial only) x 297.06 -38.43 -38.43 430.41 z 271.28 -18.09 -18.09 411.00 288.47 -31.65 -31.65 423.94
s = 1,s = 0.001 -81.66 -14.52 271.28 -18.09 35.99 -15.71
-14.68 -18.09 -15.82
-141.26 411.00 42.83
S=l,s=l x 38.91 z 271.28 116.37
-26.35 -18.09 -23.60
-28.08 -18.09 -24.75
42.13 411.00 165.09
s = 1,s = 0.0001 -79.76 -14.20 271.28 -18.09 37.25 -15.50
-14.20 -18.09 -15.50
-141.06 411.00 42.96
S = 1,s = 0.1 -58.62 2 27 1.28 51.35
-21.39 -18.09 -20.29
-21.44 -18.09 -20.32
-37.36 411.00 112.09
S=l,s=O -79.75 271.28 37.26
-14.20 -18.09 -15.50
-141.05 411.00 42.97
s = 1,s = 0.01 -84.77 -17.08 z 271.28 -18.09 33.91 -17.42
-16.71 -18.09 -17.17
-147.08 411.00 38.95
x
X
-14.20 -18.09 -15.50
S = 128, s = 10 ( {Im)} = 0, polynomial x 297.62 -38.51 -38.51 a 272.10 -18.21 -18.21 289.11 -31.74 -31.74
only) 430.50 410.66 423.89
S = 128, s = 0.01 (380 det.11642 det.) -39.10 -19.58 -17.29 -150.44 272.10 -18.21 -18.21 410.66 64.63 -19.12 -17.60 36.59
S = 128, s = 1 (4 det.112 det.)b 71.65 -27.99 -30.93 z 272.10 -18.21 -18.21 137.47 -24.73 -26.69
198.61 410.66 269.29
S = 128, s = 0.001 (4366 detJ5978 det.) -112.86 -13.00 -12.67 -213.61 272.10 -18.21 -18.21 410.66 15.46 -14.74 -14.52 -5.52
S = 128, s = 0.1 (50 det./66 det.) x -4.18 -23.77 -23.79 z 272.10 -18.21 -18.21 87.91 -21.92 -21.93
-59.54 410.66 97.19
S= 128,s= 0.0005 (6028 det./7510det.) -118.62 -12.31 -11.92 -224.20 272.10 -18.21 -18.21 410.66 11.62 -14.28 -14.02 -12.58
X
sFor the definition of xc and x0, see text. bNumbers in parentheses are the numbers of determinants generated by Lc/r& and Lo/r& for threshold s.
linear equations appearing in the SCF-CI treatment. As in the other methods, however, large basis sets of A.O.‘s are necessary if one wants to get acceptable values of magnetic properties, although the polynomial part may correct in some way the weakness of the basis. This leads to large systems of linear equations, because treatment of electron correlation (e.g., in the case of CO) converges slowly, implying the need to choose a large S-space and to lower as much as possible the threshold s.
58
3.50 3.00 P
_ "0 2 T B b
2.50 2.00 1.50
1.00 .50
: N ‘\.
-dim S=l28 ----dim SE 1
:
.
.
.
*. .
.
.
.
.
.
.
.
b
.
.
.
-t.
*.
‘.
‘bI.
*.
-.
‘.
w ‘.
-.
‘b...
Fig. 3. CO (dim. S = 1 and 128). Screening constants 6~ and 50: convergence of the calculations as functions of the selection threshold s.
-1.40 -1.60
-1.80 _
-2.00
"0 2
-2.20
-2.60 -2.80
-3.20
Fig. 4. CO (dim. S * 1 and 128). Susceptibilities fl functions of the selection threshold s.
and X0 (see definition in the text) as
59
Direct comparison with experiment needs introduction rections and will be published at a later date.
of rovibronic cor-
ACKNOWLEDGMENTS
The calculations were performed on the CRAY-1S computer of the Centre de Calcul Vectoriel pour la Recherche (C.C.V.R.). This work was supported by grant No. 3592 from the C.C.V.R. Scientific Council which is gratefully acknowledged. REFERENCES 1 (a) J. A. Pople, J. Chem. Phys., 37 (1962) 53. (b) H. F. Hameka, Mol. Phys., l(l958) 203. (c) F. London, J. Phys. Rad., 8 (1937) 397. (d) R. Ditchfield, Mol. Phys., 27 (1974) 789. (e) E. Vauthier, S. Odiot and F. Tonnard, Can. J. Chem., 60 (1982) 957. (f) H. Fukui, K. Miura, H. Yamasaki and T. Nosaka, J. Chem. Phys., 82 (1985) 1410. 2 W. Kutzelnigg, Isr. J. Chem., 19 (1980) 193; M. Schindler and W. Kutzelnigg, J. Chem. Phys., 76 (1982) 1919. 3 (a) H. Fukui, Int. J. Quantum Chem., 23 (1983) 633; H. Fukui, K. Miura and F. Tada, J. Chem. Phys., 79 (1983) 6112. (b) E. N. Svendren and T. Stroyer-Hansen, Theor. Chim. Acta, 45 (1977) 53. 4M. Iwai and A. Saika, J. Chem. Phys., 77 (1982) 1951. 5 A. J. Sadlej, Mol. Phys., 26 (1973) 1445. 6 D. F. Tuan, T. S. Epstein and J. 0. Hirschfelder, J. Chem. Phys., 44 (1966) 431. 7 A. J. Sadlej, Acta Phys. Pol. A, 49 (1976) 667. 8 B. Levy and J. Ridard, Mol. Phys., 44 (1981) 1099. 9 M. Maestro and R. Moccia, Mol. Phys., 29 (1975) 81. 10 (a) S. T. Epstein, J. Chem. Phys., 42 (1965) 2897. (b) S. T. Epstein, J. Chem. Phys., 58 (1973) 1952. 11 H. F. Hameka, Rev. Mod. Phys., 34 (1962) 87. 12 J. Tillieu, Ann. Phys., 2 (1957) 471, 631; J. Guy and J. Tillieu, J. Chem. Phys., 24 (1956) 1117. 13 M. Karplus and H. J. Kolker, J. Chem. Phys., 38 (1963) 1263. 14U. Sternberg, K. Salzer, H. Pfeifer and W. Haberditzl, Monatsch. Chem., 111 (1980) 505. 15 C. Cohen-Tannoudji, B. Diu and F. LaloE, Mecanique Quantique, I, Hermann, 1973, p. 314. 16 J. Olsen, D. L. Yeager and P. Jorgensen, Adv. Chem. Phys., LIV (1983) 1. 17 W. H. Louisell, Quantum Statistical Theory of Radiation, Wiley, New York, p. 136. 18 P. W. Langhoff, M. Karplus and R. P. Hurst, J. Chem. Phys., 44 (1966) 505. 19 T. P. Das and R. Bersohn, Phys. Rev., 104 (1956) 849. 20 U. Sternberg, W. Haberditzl, Mona&h. Chem., 106 (1975) 701; 109 (1978) 735. 21 U. Sternberg and H. Rosenberger, Chem. Phys., 55 (1981) 275. 22 A. Y. Sadykova and R. M. Aminova, J. Struct. Chem., 21 (1980) 142. 23 A. Okninsky and A. J. Sadlej, Acta Phys. Pol. A, 42 (1972) 709. 24 B. Huron, J. P. Mahieu and P. Rancurel, J. Chem. Phys., 58 (1973) 5745. 25 A. Okninsky and A. J. Sadlej, Mol. Phys., 26 (1973) 1545. 26 J. C. Lie and E. Clementi, J. Chem. Phys., 60 (1974) 1275. 27 A. J. Sadlej and W. T. Raynes, Mol. Phys., 35 (1978) 101. 28 P. Lazzeretti, R. Zanasi and B. Cadioli, J. Chem. Phys., 67 (1977) 382.
60 29 S. Huzinaga and Y. Sakai’, J. Chem. Phys., 50 (1969) 1371. 30 C. Salez and A. Veillard, Theor. Chim. Acta, 11(1968) 441. 31 C. Brezinski, Algorithmes d’accel&ation de la convergence. Etude numbique, Editions Technip., 1978. 32 P. W. Fowler, G. Riley and W. T. Raynes, Mol. Phys., 42 (1981) 1463. 33 J. A. Pople, W. G. Schneider and H. J. Bernstein, High Resolution Nuclear Magnetic Resonance, McGraw-Hill, New York, 59. 34 A. Carrington and A. D. McLachlan, Introduction to Magnetic Resonance, Harper and Row, New York, 1967, p. 55. 35 T. D. Gierke, W. H. Flygare, J. Am. Chem. Sot., 94 (1972) 21.7277. 36 A. Velenik and R. M. Lynder-Bell, Mol. Phys., 19 (1970) 371. 37 T. Weller, W. Meiler, H. J. Kohler, H. Lischka and R. HSller, Chem. Phys. Lett., 98 (1983) 541. 38 R. Jaquet, W. Kutzelnigg and V. Staemmler, Theor. Chim. Acta, 54 (1980) 205.
GAUGE-DEPENDENT RET CALCULATION OF THE MAGNETIC PROPERTIES OF MOLECULES Part I. Method and application to Hz, Hz0 and CO
J. P. FLAMENT, Groupe D.C.M.T.,
H. P. GERVAIS and M. RERAT Departement
de Chimie, Ecole Polytechnique,
91128 Palaiseau (France)
(Received 15 September 1986)
ABSTRACT A gauge-invariant ket function theory of magnetic properties is presented. The perturbed state, which depends explicitly on the strength of the external fields, is determined variationally as a linear combination of Slater’s determinants associated with a gauge transformation factor. This “mixed” procedure combines advantages of both the polynomial and the familiar SCF-CI methods; it is applied, with inclusion of electronic correlation, to H,, H,O and CO. INTRODUCTION
Considerable interest has been devoted to the quantum mechanical calculations of the magnetic properties of molecules. The magnetic field does not appear in the theoretical ex ressions through its strength ‘B, but through its associated vector potential i for which one has to choose a particular gauge origin, i.e., a position vector Luch that A = l/23 A (i’- A). In the case of a molecule, there is no natural choice; use of infinite basis sets should lead to gauge-invariant results, especially with respect to the origin; in practice, sufficiently large basis sets that make the gauge problem irrelevant are hardly attainable. To solve the problem of translational invariance, an interesting approach uses expansions of the molecular orbitals over a set of gaugeinvariant atomic orbitals (GIAO) [l] , but it needs calculation of additional two-electron integrals and, to our knowledge, correlation effects have never been considered; the IGLO method [2] has been proposed to overcome these difficulties. When correlation effects are considered, GIAO are not used and large basis sets are necessary; gauge invariance however is not obtained [3a, 4-61. In the GIAO or IGLO methods gauge invariance is guaranteed; both in fact involve a “constrained” gauge, in the sense that the origin of a correcting factor is set either at the nucleus of the various atomic orbitals (GIAO) [l] or at the centroid of charge of localized orbitals (IGLO) [2]. Some other methods optimize the position of the origin by minimizing the paramagnetic contribution, eventually with inclusion of correlation effects [ 7-91. 0166-1280/87/$03.50
0 1987 Elsevier Science Publishers B.V.
40
It has been shown that constrained gauge invariance through GIAO’s does not ensure current conservation [lob] ; in the same paper, it has been proposed to consider the gauge function itself as a variational parameter, as already suggested and applied by others [ 11-131. This last procedure usually leads to acceptable results for proton magnetic shielding, but fails in the case of heavier nuclei [ 141: one major reason for this is the difficulty to find a suitable gauge function which accounts for the existence of nodal surfaces in the wavefunction when heavy nuclei are present. In contrast, customary perturbation methods would be well-adapted to this situation, but they unfortunately do not ensure gauge invariance (which could only result in a complete basis). The computational procedure we propose here follows Epstein’s suggestion [ lOa] : to translate both the “geometrical” effect of a gauge transformation and the proper physical influence of the field B,,, one has to consider the gauge factor 0 as well as the ket 19) as variational parameters in the expression I$) = e’“l$)
(1)
of the perturbed
molecular
state vector.
METHOD
Unitary transformation of the states Let 10) be the normalized ground state of the system in the absence of any magnetic field (10) belongs to the orthonormal set { 1~))formed by the eigenstates of the field-free hamiltonian A,,). The gauge transformation operator 0’ takes the familiar form
where 6’ (resp. &) is hermitian [I51
and a function
of the electronic
coordinates
We may consider the ket I$), appearing in eqn. (1) and depending on-the vector potential 6, as resulting from the action of a unitary operator ezs on the field-free ket IO) [16] I\ir) = p(e)
= eieee’glO)
For the field-including E(a)
= (*IQ*)
(41 hamiltonian
A, this leads to an energy
= (Ole-‘S^e-‘BAe’6eis^10)
In the following,
we shall be concerned
(51 by expressions
of the form
=
C
a,y:
(6b)
P
By hermiticity of @, we have a,, = a;. The development of the hermitian operator, $ = q/a li Sl, IlXmI can be limited to the part operating on IO> s =;
=;
c (S,lm>(Ol + S~lO>W) m#O
c c%nfik+ WL)
(7)
mf0
It acts as mixing the states {Ii>} under the influence of the external state field. In the magnetic case, the perturbation operator contains an imaginary part, so that the S matrix is not reducible to its imaginary component, as in
[W. Variation of the energy Using the development [17] : e&fie’ = & - i[a, 81 + 1/2[a, [g, A] ] + . . . one gets from eqn. (5) for the total energy associated with the ket I$> E(A, 6, $) = (Ole-iS^e-ibAeiQei~lO) =
[A, &:I ] IO> + 2(OI[& [A, 611 IO>}+ . . . (8)
that is, in matrix form, by use of eqns. (6) and (7)
E=Eo--i(eP,S,,,,S$)~~+
1(a~,S~,S~)~~~~~~~~~~)
(9)
with E. =
(01 [fp, A] IO>
VI, = ; (01 [I?;,
(Ad,,
=.$
(A&m
P%
(01 [?p,
?“I 1 10)
=$OI[[?,,,
A],eJl0~
A] IO>
v,,
=;
(A13),m
(01 [II,,
=$
W[[?j,,
Al, %ll IO)
A] IO>
42
&dm,,
= $
(&dmn
(01 [fik,
(01 &n,
(&A,n
[A, R’,]] IO)
rmJl10) (&~hfi = $ (01[km [A,
=$
(A32hn
=$
colr&l,
(01 P%,
[A, R, ] ] IO) (A33L
[A, R; ] ] IO)
$11 10)
=$
=$
COl[ll,,
[A, I?“]] IO>
The A, matrices satisfy the following relations A22
V, = -v;k=
Vand
=&=fj=()
Ai2 = A,, = A,*, = A& = C
i A2,=Af2=D from which, putting AlI = A and S, = X, + iY, (X,,
Y, real), we get
E = E,, .- i(a, S,, Sg)
E = E0 - 2(a,, X,,
=
Y,)
E,, - 2iiv + i.iMu
with obvious definitions for U, v and M. The stationarity conditions for E with respect to a,,, X, and Y, is thus Mu=v
(11)
Order analysis of the matrices Since we are interested in E due to the changes induced into # by the presence of an external field, we now develop eqn. (11) in powers of the perturbations and solve it for each order separately. For this purpose, every essential quantity will be expanded as follows
v=v,+ v,+ a=a,+
v2
al+ a2
43
x=x,+x,+x* Y= Y,+
y,+
Yz
As already mentioned, 10) represents the ground state of the field-free hamiltonian A,-, and thus belongs to order 0. Furthermore, as in other gaugeinvariant calculations [l, 2, lo], we restrict the gauge operator &’ to be a linear function of the positions (?j) of the electrons, i.e., with eqns. (3) and (6) e’=
c al,$I,withZ+
m+ n=l
(12)
lmn
We now concentrate upon the magnetic case (but the forthcoming developments may be easily adapted to the case of an external electric field). Zero th-order MouO = u.
(13)
Using the common commutation rules, we have w; = ; (01EM, A01IO)= -iq(OI($
lgrr
+ar,*$)lo>
For 5, restricted by eqn. (12), 8’: reduces to WE= -Siq
When 10) and Im>are exact solutions of the field-free Schrijdinger equation, V&,, disappears. In current calculations, however, where these kets are only approximations, this normally no longer holds; fortunately, if IO>approximates sufficiently well the exact ground state, V&,, should be small, so that we shall neglect this term in the following. It thus follows that u. = 0 and u. = 0. First-order The first-order equation is Maul +Mluo = ul. On the basis of the preceding approximation, this reduces to M,u, = u1 where
(14)
CL
i (P
=COl&lO>) ,
Equation (14) then factorizes into Y1 = 0 and
Up to now, no hypothesis has been made about the nature of the vector potential A; the foregoing development thus applies to any kind of static field. In the general case, d may be a sum (e.g., of an external magnetic field B. and of the field due to a nuclear spin). As long as the fields are independent, eqn. (14) can be split into as many equations. It should be noted, however, that a “coupling” appears between the gauge and the field through matrix CL so that we may split the gauge contribution a, into parts related to the relevant fields. This means that, when considering chemical shift, the gauge function G’ does not depend solely upon the external magnetic field (a line already suggested by Hameka [ 111). Second-order We now focus on magnetic susceptibility appear at second order. Let us define
and magnetic shielding, which
and
(17b) where 8 is a homogeneous magnetic field, 3 the magnetic moment of nucleus N centered at &, q the position of electron j about some fixed origin, F”j = i$ - Rt, the position of electron j relative to N. Substitution of
45
into the equations leads to the calculation of x and u. With u solution of eqn. (ll), and recalling that u. = u. = 0, we get to second order E = COlA,lO>+ (OIA1lO) +
since fi, is purely imaginary
and
E =
Et')=(01C &?lO>+ (01C x'210, + 2(01 C & i - (i+‘l
i
l
$‘lO, -i&u;--i?al’u;
i
+ ii:‘u\)
With the use of eqn. (14), we may now demonstrate
the exchange
c’lv;’ = ti:‘u; since we have (the first-order
theorem (19)
eqn. (14) is linear in the perturbation)
M,u; = u;
PO@
M,u'l= u;'
POb)
Multiplying eqn. (20a) on the left by ii;’ and eqn. (20b) by ii:, and using the fact that MO is symmetric, we get eqn. (19), which is an extension of the usual exchange theorem, including the gauge. to Vectors u’, and ~9, as well as vi and uy, are to first order proportional IIEfll = B or II311 = M. We thus write a’=+e~
(21)
where Pis a unit vector along g or $. With this notation we get from eqns. (20a and b) a new set of equations where the coefficients are the same for both cases with eqns. (20a) and (20b) differing only by the type of operator (ds or la,) appearing on the right hand side
Using the {(I,,, x,} solutions
lG2) =
eg
of eqn. (22), the second-order
ii’u’ +
*
energy is
-$-* mffi’u” + . , . (23)
=
-~r3(xd+ f)S
+ B(od+ CJP)~ + . ..
where suffixes d and p denote as usual the diamagnetic and paramagnetic parts of the susceptibility and screening tensors. Using the definitions of eqn. (21) for ds and A,, it follows that
(244 1 1 -Zx~=-.-$&=--,
1 P
ap WkBfi+ C Xe VzP m
(24b)
=--
(24~)
[q A (Gj7p)]plO) + 2x m
xf$tOlC(3
A
(24d)
j
WW
=-g(cCLaMP(OI w 1 (I$ A (Gj,))&lO) i
+ 21 X~‘(OIC(~ m j
A Gj),lm,> (25e)
Equations (25d+e) clearly show the opportunity given by the exchange theorem. In the former, eqn. (25d), ap and X”,” are obtained when solving
eqns. (14) or (20a) for the perturbation BO; in the latter, eqn. (25e) afP and Xzfi are obtained when solving eqns. (14) or (20b) for the perturbation Mk Knowledge of (allBar, Xf) alone leads to both x and (I.
Connection with other methods By making specific choices of G, S and IO>,the above procedure reduces to some previously proposed methods: (1) Setting G = 0 and IO>= IHfl, the states Im>appearing in S are simply the monoexcited determinants ii + p). Vh = 0 then is no longer an approximation, since it expresses Brillouin’s theorem. This situation is closely related to the Tamm-Dancoff approximation (TDA). The equations are readily written
C [6,,Sii(e,-eEi)+ ia
2(ipljq)-(iqIjp)]X,j=--V&
which is an approximation to the CHF equation [ 181 E
jq
[6p*6ij(fp--i)
+ (WA--(0
IPq)lxqj=-v&
obtained with &, = hii = UiUi instead of Fiji = e$ilHF)(HF(. The CHF procedure is related to the random phase approximation (RPA) in which the B’matrix is not null. The decoupling has been introduced when eliminating the terms Im>(Zl with I # 0 in S, an operation which suppre_sses_ter_ms like S om X S,, in the development of the double commutator [S, [H, S] ] (although they do not act on IO), they introduce a coupling between the operators Im>(01). (2) Setting G = 0, with IO>being a correlated wavefunction, one gets Fukui’s SCF-CI method [3a] (see also Svensen [3b] for the case of electric polarizability). (3) Setting 10) = IHF), S = 0, and varying G alone for some analytic forms leads to various developments of gauge transformations [ 12-14, 19-221. (4) Finally, if we expand the exponentials to first order we get a wavefunction of the type I@‘) = i(&lO> + &$,In>). This is closely related to Okninsky-Sadlej’s method [ 231 where 19) refer to orbitals, while in the present paper it symbolizes states or Slater’s determinants).
Practical implementation The correlated ground state IO) has been determined by means of the CIPSI algorithm [24], which builds the wavefunction by iterative selection of the most important determinants; the final step included all those which had a coefficient superior to ca. 0.01 in the M$llerPlesset first order wavefunction. Use of states Im>in S would require a lot of computation, so we generated them in the spirit of CIPSI and perturbation theory: once IO>has
48
been built as a combination of Slater determinants (mainly 1HF) and diexcitations), all mono-substitutions are introduced to generate all mono-, di- and tri-excited determinants orthogonal to 10) (i.e., that do not2ppear in its development) which interact with it b+y means of operators L (coefficient sB or sL) = l(OlLlm)/& - ,?&I) or LN/r& (coefficient sM (or S&$) = I(01k N/r&lm>/Eo - &I) (both criteria can be used indifferently because of the exchange theorem). We then make the approximation (Olfi,,lm) = 0, which is not a serious drawback, since the most important determinants of 10) have already been taken into account (moreover, in highly symmetrical molecules this matrix element is generally zero for symmetry reasons). y,, being limited to monomials of degree 1 (cf. eqn. (12); y = 3c,y, z) coefficients appearing in the two parts of eqn. (22) write in the case of perturbation B, (n, is the number of electrons): 0
n,
2 0
-m$o
.
-&IO)
. ..=
0
-(ml
0
-(mldy d IO>
. ..=
-_(ml$lO)
. ..=
%
(mlHO-E”ln) -m-go
0 -
cd 2 (Y) T
(26)
. . . = -my$--s&O .
The closed shell ground state IO>belongs to the totally symmetric representation of the molecular point group; for the perturbation B,, the kets In> belong to the irreducible representation l’(L,). If the molecular point group is such that I’(&) differs from I’(x, V,), I’(y, Q,) and I’(z, VZ), the coefficients a, of the gauge function vanish (a,@) = a,,(%)= a?) = 0) when the origin is the electronic centroid ((x) = (y) = cZ>= 0). When this holds also for the perturbation B, and B, (as in the case of D-h), the problem reduces to its SS part and the gauge function is zero [25]. For less symmetrical systems, however, the gauge function must be taken into account: in the water molecule, for example, taking Oz as the CZ axis, only B, obeys the preceding rules if the origin lies on that axis and a GS coupling appears for the B, and B, components; when we move the origin away from the z=axis (e.g., from the oxygen to one hydrogen), new determinants will be generated by g,(H) - L,(O) = XOH a /i3y and gauge invariance may not be guaranteed if one considers only the SS part of eqns. (26). In contrast, if the GG and GS parts are explicitly taken into account, the gauge function absorbs completely the translation of the origin and the coefficients of the extradeterminants disappear; thus suppose that X, = XL (t/m) when
49
translating the origin, the G lines give a:--ua,=
0
a;-aa,=Zon
(27a)
a: - a, = -Yen 1 whereas for one S line -_(mlV,lO>(a: -a,)
- (mlV,IO)(a~ - uy) - hlV,I0>(u:
- Z,,Qrw,lO)
-a,)
= Y,,GnlV,I0> Wb)
Equation (27b) still holds after substitution of eqn. (27a). Consequently one may put the origin at the point of maximum symmetry, in order to reduce the number of kets lm>generated by perturbation. APPLICATIONS
Three procedures have been compared, in the case of the small molecules Hz, Hz0 and CO (at a fixed geometry, that is without introduction of rovibronic corrections) to the “polynomial” method (a f 9, X = 0) with degree p superior to one for the gauge function G, the “SCF-CI” method (a = 0, X # 0) and the “mixed” method outlined above (a # 0 with p = 1, X f 0). Hz (dnn = 1.401 u.a.) Two basis sets were used: (A) that of Lie and Clementi [26] ; (B) that of Iwai and Saika [4] (i.e. (A) augmented by (Is, 2p, Id) on each hydrogen). The ground state wavefunction was obtained through SCF (HONDO) and CI (CIPSI) calculations, as described in the Practical implementation section. With our actual program (diagonalization below dimension 200), a complete CI was possible for basis (A); for basis (B), more than half the CI space (S = 196) was diagonalized, the second-order M$ller-Plesset barycentric correction to the energy [24] due to the remainder being less than lOA u.a. As concerns the polynomial method (a f 0 with p = 2, X = 0), Fig. 1 shows the behavior of a and acX obtained, while increasing the dimension of CIPSI’s S-space. In the case of basis (B), one notices a sudden rise of energy at S = 196 (account being taken of the few mono-excited configurations which appear at this stage). Results from the various procedures are given (basis (B), S = 196) in Table 1. In the case of the polynomial method, gauge invariance is respected when translating the origin from the center of mass Sl to nucleus H. Thus AG, = (0.596yH - 0.149ynzu) = (0.596y - O.l49y(z = (0.596
-
0.149Z)y
- (-0.149yz)
+ 2) + 0.149yz =
+ 0.149 $
>
y
Fig. 1. Polynomial method (H,). Evolution of 6~ and or(H) dimension of space S.
(Basis sets A and B) versus
Development of & to the third degree is useless, since in Dmhgroup no monomial of this degree belongs to the same irreducible representation as fl,, &, L,. Development of & to degree one implies that the paramagnetic contribution to magnetic shielding is zero when the origin is put at the electronic centroid (a here). Comparing the SCF-CI and the mixed methods (the threshold & or sL for selection of the excited states lm> being identical and small, s < 10e4 u-a.: nearly all kets lm> are introduced), one sees that both give exactly the same result when gauge origin is set in a. Gauge invariance still holds for the mixed method (the results depending, of course, on the basis chosen), but not for the SCF-CI, even with a large basis and CI calculation. Finally, values ii, = 26.64 ppm and u%“(H) = -8.61 ppm, as given by the mixed method, can be compared to other values mentioned in [ 271 (although a rovibrational correction would be necessary for comparison with experiment). By contrast, the polynomial method, even if ensuring gauge invariance, is far from likely to give good results; this conclusion will be more explicit below, where heavy atoms are present. Hz0 (doH = 0.958 A; HbH = 104.5”) Satisfactory values of magnetic properties need to be calculated with large basis sets. We thus choose basis (IV) employed by P. Lazzeretti et al. [28]
51 TABLE 1 Hydrogen molecule (A) Polynomial method Monomial basis (degree < 2)
Gauge origin: H
Gauge origin: n a
G,
G,
0 0.596 0 8 0 4.149 0
34.63 0 0 -9.04 0 0 25.59 0
“H
2i.35
GY
0 0 0 0
-0.596 0
0”
z
3: 63 0’ 0 -9.04 0
-:
0.
00 8 0.149
149
0 0
25.19 0 0
0 0
27.88
0 25.19 0
0.40
8 0
2: 59 0’
0 0
8
0.149 0 0
840 0’
: 25.59 0 0 26.35
8 27.88
GZ
GY
2: 59 0’
0 0
27.88 0 0 0 0 0 27.88
For a polynomial of degree 1, &@(sL) = 0 v (a, p) 4 6~ = 26.09 ppm For a polynomial of degree 3, GH = 26.35 ppm (B) SCF-CI method (threshold s = lo-+ u.a.) Determinants introduced by Lp uv = 2981146 a;” = 2981146 a? = o/o
Gauge origin: H -8.36 0 0
0 -8.36 0
OH
26.27 0 0
GH @pm)
26.81
(C) Mixed method (threshold s = lo+ gauge origin OH =
26.02 8
PXX CH (H) =
0 26.02 0
-8.61
0 26.27 0
u.a.; polynomial
Gauge origin: ~2 0 0 0 0 0 27.88
0.83 0 0
0 0.83 0
0 0 0
26.02 0 0
0 26.02 0
0 270.88
26.64 degree l), results independent of
0 0 27.88 >
ppm
&“(
n ) = 0.83 ppm (for this choice of gauge origin, all coefficients zero; the value of 0.83 ppm is purely SCF-CI) OH- 26.64 ppm
%enter of mass of the molecule.
in the polynomial m
52
and consisting in the [6s, 5p] contraction of the (lls, 7~) of Huzinaga and Sakai’ [29], increased by a 3d (STO-SG, 5 = 1.93) for the oxygen atom and [ 2s, 2p] A.O.‘s of Salez and Veillard [ 301 for the hydrogen atoms. Trying to improve gauge invariance in the case of the pure SCF-CI method, we added to this a set of four d shells (5 = 1.0) centered in the molecular plane, 0.5 u.a. away from the oxygen atom on the coordinate axes. This did not prove to bevery useful, but we choose first to present (in Tables 2 and 3) a comparison of the three methods via the results computed on this basis of 62 A.O. In Tables 2 and 3, the S space of CIPSI has a dimension unity (i.e., the diagonalized energy is that of the SCF calculation and is taken as zero energy, TABLE
2
H,G: ,,@ oru) uorp and Cr Compa%on Af results given by the SCF-CI 4d, dim. S = 1, E, = 0,threshold s = 0) (A) SCF-CI:
P
x
(d)
x
(p)
Y Y
Y
(B) SCF-CI:
0
0 0 0
-39.21 39.04 -0.17
0 0 0
101.05 -59.01 2.05
17.62 -7.35 10.26
415.21 -53.79 361.42
(d) (P)
00
(ppm)
0)
0 0 0 ; 0 325.85
z
X
0
21.14
0
- -11.19
0
9.95 0 0 0
8.38 12.82 21.19
0 0 0 417.18 -112.44 304.74 0 0 0
0 0 0 0 0 0 415.74 -104.35 311.38
-10.85
0
23.11 6.14 29.25
0 0 0
Y
44.80
basis +
39.08
2.97
-7.04
C(
1,1,l) z
1.04
0
38.19
40.69 (C) Polynomial
z
Y
(d)
s z t
(Lazzeretti’s
Gauge origin:
29.32
u,, (gauge origin:
(P) (d) (P)
31.67 5.83 37.51
44.06
X
0
Y
0 0 0
x
Y Y Y
X
130.32 -85.79 44.53
i
x
X
z -39.21 39.25 0.05
P
ct
Gauge origin:
0 0 0
(d)
(PPm)
H
Y
75.67 -30.08
(P)
z z z tin
x
45.59
X
methods
oH Gauge origin:
(Y
and polynomial
method
Degree of polynom.
OH @pm)
ho (ppm)
1 2 3
24.84 26.32 27.03
415.13 415.04 414.42
53 TABLE 3 H,O: c$$orn), 0~ and 6, results given by the mixed method (Lazzeretti’s basis t 4d, perturbation, dim. S = 1, E, = 0,threshold s = O)*
Gaugeorigin:
P
x
H
Gaugeorigin: 0 z
Y
x
Y
L,
Gauge origin:C(1.1.1) I
x
Y
z
b
OH 1 2 3 4 ~~(126
2.10 0 0 0.99 0 75.67 0 -37.32 det.) 38.35,0,8.92 (37.36,0. 13.35)
1 2 3 4 oyp(156 det.) 1 2 3 4 a&96
det.)
+(PPm) 00
0 9.36 -1.00 0 0 130.32 0 -108.38 0. 21.94.0 (0.14.06.0) -7.35 0 0 -1.43 -39.21 0 49.47 0 10.26,0. 29.25 (17.62,0. 23.11)
-5.95 0
-39.21 48.13
0 1.43 0 0
6.14 0 101.05 -71.80
29.85
0 2.10 -0.12 0 0 31.67 0 6.67 38.35.0.8.92 0 9.36 0.11 0 0 8.38 0 13.56 0, 21.94,0 0 -7.35 0 0 0 17.62 0 -7.35 10.26,0. 29.25
-5.95
0 21.44 -12.22
0
0 0 0
6.14 0 23.11 6.14
2.10 0 0 0.88 71.37 -54.45 -93.02 54.45 38.35,0. 8.92
0 1.00 0 0
-7.35 0 1.00 -1.00 -22.08 -39.70 32.34 39.70 10.26.0.29.25
6.14 0 77.56 -48.31
b
1 -47.12 2 0 3 404.59 4 -37.60 o,p(126 det.) 366.99,0,0 (413.84.0.0) 1 2 3 4 oYp(156 det.) 1 2 3 4 a,+(96 det.) oo(PPm)
0 0.99 0 0
0 0 0 0
0 -107.69 -1.00 0 0 406.56 0 -98.09 0,308.47,0 (0.415.81.0)
0 1.43 0 0
0 0 0 -1.43 13.71 0 -13.71 0 0. 0, 311.38 (0.0, 415.73)
-104.35 0 415.73 -104.35
328.85
-47.12 0 -0.12 0 415.21 0 -48.22 0 366.99, 0.0
0 0 0
0 -107.69 0.11 0 0 417.18 0 -108.70 0, 308.47. 0
0 0 0
0 0 0 0 0 0 0 0 0, 0, 311.38
-1.00 -33.31 42.23
0 9.36 -0.89 0 0 102.52 0 -80.58 0, 21.94.0
29.85
29.85
-5.95
0
0
-104.35 0 415.73 -104.35
328.85
*Values in parentheses are given under the same conditions by the polynomial method (degree 1). bl. SCF-CI paramagnetic contribution. 2. Coefficients of the polynomial. 3. Diamagnetic contribution. 4. Paramagnetic contribution. One can verify that a change 3 o f gauge origin induces variations in the polynomial coefficients, exactly equal to the corresponding coordinates of $, thus satisfying AG, = (3 A &. (The position of the hydrogens is (1.43; 0; ll.lll), in ua., with respect to the oxygen.)
54
E. = 0). In this case, the excited kets lm> generated by B, (resp. M,J are not too numerous for the resolution of the linear system eqn. (22), even if sBa = 0 (resp. sMa = 0). Gauge invariance (AG, = @A it),), as well as the exchange theorem, can then be verified in Table 3 for the mixed method. When the S space expands up to dimension 192 (E,-, then becomes the lowest eigenvalue of fi,, projected onto it) it appears that & generates too many kets lm> for a current resolution of the linear system of eqn. (22) if stcu = 0. If sLa # 0, however, gauge invariance is not automatically ensured, because perturbation L, depends on the origin, as well as the kets Irn) generated by it. To avoid this difficulty, we choose sM as a threshold for selecting Im),even in the computation of susceptibilities (sM then has to be as low as possible): for HzO, we calculate the magnetic shieldings cn of hydrogen and co of oxygen, after selecting excited kets Im>,respectively with the help of the values (denoted by xH and x0) of SLY/,+ and of SL~/'~, and then compare magnetic susceptibility obtained via these two basis sets. Table 4 reports the results obtained for various values of sM, in Lazzeretti’s original basis of 42 A.O. Figure 2 is a plot of cn, uo, xH and x0 versus colog s,. Convergence is relatively slow, which requires resolution of large systems TABLE 4 H,O: @(CH), @(Go) an d x‘@(x), results given by the mixed method (Lazzeretti’s basis, L&h perturbationa, dim. S = l-192, for various selection thresholds) x(10-"J T-2)
~Hbpm) S-l.s=0 36.37 0 0 21.77 10.35 0
8.93 0 29.23
29.79
-13.78 0 0
0 -14.01 0
0 0 -14.64
-14.38
;,=,‘,““9 ‘;H/& 21.38 0
0
0 311.36
-14.73 0 0
0 -16.62 0
q)(PPd 0
0 -15.59
-15.65
s = 192.q.J= 0.001,90 = 0.005b 38.31 0 9.57 -14.15 0 0 21.37 0 0 -14.35 11.49 0 28.79 0 0
29.51
0 0 -15.59
-15.65
29.49
0 309.01 0
~'(10-~J T-*)
xH(lO-" J T-*)
25.59
366.72 0 0 329.03
S= 192,s= IO((I~)}= O; polynomialonly)b 37.73 0 12.84 -14.73 0 0 15.08 0 0 -16.62 17.07 0 23.77 0 0
10.90
0 0 -14.02
-13.94
aH(PPm)
0
qJ(PPd
-14.19 0 0
413.88 0 0
0 415.66 0
0 0 415.71
0 303.61 0
0 0 312.01
0 311.17 0
0 0 323.02
415.08
0 -14.35 0
0 0
-14.39
-14.32
371.10 0 0 328.91
= o.oo5,sLo = O.OOlb 9.42 0 28.87
-14.13 0 0 -14.27
0 -14.30 0
0 0 -14.39
-14.11 0 0
0 -14.12 0
0
0
-14.28
-14.15
aExcept for the last case. bFor the definition of xH and x0, see text.
375.08 0 0 336.42
55
‘. ‘.
*.
“....................,...,*
341.55
CT c
UH
346.56
_.,’
g
%.,
337.40 336.42 .._,# . . . . . . . *n 328.91 ~..____........... .a , . 29.44 29.49 29.51
25%9
X
z :: ;;
-14,3&y,. -14.54 .y:.. / , ,-’
-Mll, ..
-15.65...-.“’ q. .
-15.33 .=
I 0
I -I
’ ..
-15:65 I I
-14.15 .m -1y.27 ..*-. -14.38
XT, . ‘,.i$O XH.
I 2
I 3
colog
5
Fig. 2. H,O (dim. S 0 192). Convergence of the calculations as function of the selection threshold s of the excited kets Im).
of linear equations; a good procedure for this, especially when a lot of matrix elements are zero, is the Gauss-Seidel method (where only non-zero elements are stored, thus needing less memqry and less I-O operations). More, over, the c-algorithm [31 j is well suited to this method and can be added to accelerate convergence of the iterative computation; in this way, systems of the order 10 000 containing sometimes 99% of zeros can be solved. Other calculation results, as in refs. If or 32, can be referred, as well as the experimental values mentioned in ref. If (an = 30.2 ppm; o. = 334 ppm) or refs. 33-35. We must emphasize that the value u. = 334 ppm is obtained from a private communication mentioned in ref. 36 giving the quantity (%O
-
aH,O)*
For H20, it appears that the S space (i.e., electronic correlation) does not contribute very much, so that one can be content with the results obtained for S = 1 and s = 0. co (dco = 1.13 A) For the oxygen atom, the basis chosen was that already used in the case of HzO; for the carbon atom, it was a [6s, 5p] contraction of the (115, 7~) basis used by Huzinaga and Sakai’ [29], increased by a d shell ({ = 0.7) (see ref. 38). With this set of 52 A.O.‘s, we obtain Es,, = -112.777876 hartree. In this demonstration calculation (comparison with experiment would need rovibronic correction and will be described in a further paper), dimension of the S space was limited to 128. As in the case of HzO, we had to choose sM instead of sB as a threshold for selecting the excited kets Im>, when dim S > 1; corresponding values of x therefore are denoted xc and x0 in Table 6.
56
Table 5 compares results given by the SCF-CI and mixed methods for dim 5’ = 1 and s = 0. As expected, both are close when gauge origin is put as the center of electronic charge G. Table 6 shows the influence of selection threshold s for dim S = 1 and 128 (see also Figs. 3 and 4). In contrast to the case of HzO, the S space (i.e., introduction of electronic correlation) deeply transforms the results, e.g., for the dipole moment, the magnetic shielding constants bc and Do are almost zero. Calculations with dim S 9 1 and s 2 0, implying resolution of large systems of linear equations, thus seem necessary in the present case. Comparison with experiment (By = -16 ppm, ~c = +3 ppm) [35] still appears unrealistic, because it requires introduction of rovibronic corrections, which can be important (see other theoretical calculations in ref. 37). CONCLUSION
The mixed method, as it was outlined and applied above, originates in Epstein’s suggestion [lob] to ensure current conservation together with gauge invariance. Perturbation due to the field is treated as in SCF-CI procedures (thus introducing electron correlation, in a manner similar to that of CIPSI algorithm), whereas a polynomial part “absorbs” the effect of the gauge. (as in Guy’s and Stephen’s procedures). In the SCF-CI aspect, the perturbed ket should normally have the form I$> = Z I, ,,,Slm In (m I; in the present treatment, we limited ourselves to 14) = Z,,,S,lm)(01 + c.c., thus neglecting terms like (ml&lO) (m f 0). The results show that this approximation has a negligible effect. Programmation aspects of the mixed method are relatively simple, because they only involve adjonction of three columns and three lines (for a polynomial of degree one, which is the most convenient case) to the system of TABLE
5
CO: oF(CC), @(CO), xw(x), comparison of results methods (dim. S = 1, s = 0) (U in ppm, x in J T-*)
given by the
SCF-CI and mixed
SCF-CZ
C!J(X = -93.14 DC OXX = -65.44 =C G”X = -77.47 UC
GE = 27.93 Czz = 271.28 OC
* Oc = 37.26
X xx = -14.20
x 1O-5
X *= = -18.09
x lo-$
u8
= 411.00
=+ Gc” = 46.80 52 = 38.78
Mj$ed method ax =-79.75 u”c”= 271.28
o”o” = -141.05
= og”
* x = -15.50
==+Trg = 42.97
x lo-”
57 TABLE 6 co: “?(Cc), c$j%o), Xoa(X), results given by the mixed method (L.&k influence of dim. S and of excitation threshold sa
perturbation),
XC
S = 1, 8 = 10 ( {Im)} = 0, polynomial only) x 297.06 -38.43 -38.43 430.41 z 271.28 -18.09 -18.09 411.00 288.47 -31.65 -31.65 423.94
s = 1,s = 0.001 -81.66 -14.52 271.28 -18.09 35.99 -15.71
-14.68 -18.09 -15.82
-141.26 411.00 42.83
S=l,s=l x 38.91 z 271.28 116.37
-26.35 -18.09 -23.60
-28.08 -18.09 -24.75
42.13 411.00 165.09
s = 1,s = 0.0001 -79.76 -14.20 271.28 -18.09 37.25 -15.50
-14.20 -18.09 -15.50
-141.06 411.00 42.96
S = 1,s = 0.1 -58.62 2 27 1.28 51.35
-21.39 -18.09 -20.29
-21.44 -18.09 -20.32
-37.36 411.00 112.09
S=l,s=O -79.75 271.28 37.26
-14.20 -18.09 -15.50
-141.05 411.00 42.97
s = 1,s = 0.01 -84.77 -17.08 z 271.28 -18.09 33.91 -17.42
-16.71 -18.09 -17.17
-147.08 411.00 38.95
x
X
-14.20 -18.09 -15.50
S = 128, s = 10 ( {Im)} = 0, polynomial x 297.62 -38.51 -38.51 a 272.10 -18.21 -18.21 289.11 -31.74 -31.74
only) 430.50 410.66 423.89
S = 128, s = 0.01 (380 det.11642 det.) -39.10 -19.58 -17.29 -150.44 272.10 -18.21 -18.21 410.66 64.63 -19.12 -17.60 36.59
S = 128, s = 1 (4 det.112 det.)b 71.65 -27.99 -30.93 z 272.10 -18.21 -18.21 137.47 -24.73 -26.69
198.61 410.66 269.29
S = 128, s = 0.001 (4366 detJ5978 det.) -112.86 -13.00 -12.67 -213.61 272.10 -18.21 -18.21 410.66 15.46 -14.74 -14.52 -5.52
S = 128, s = 0.1 (50 det./66 det.) x -4.18 -23.77 -23.79 z 272.10 -18.21 -18.21 87.91 -21.92 -21.93
-59.54 410.66 97.19
S= 128,s= 0.0005 (6028 det./7510det.) -118.62 -12.31 -11.92 -224.20 272.10 -18.21 -18.21 410.66 11.62 -14.28 -14.02 -12.58
X
sFor the definition of xc and x0, see text. bNumbers in parentheses are the numbers of determinants generated by Lc/r& and Lo/r& for threshold s.
linear equations appearing in the SCF-CI treatment. As in the other methods, however, large basis sets of A.O.‘s are necessary if one wants to get acceptable values of magnetic properties, although the polynomial part may correct in some way the weakness of the basis. This leads to large systems of linear equations, because treatment of electron correlation (e.g., in the case of CO) converges slowly, implying the need to choose a large S-space and to lower as much as possible the threshold s.
58
3.50 3.00 P
_ "0 2 T B b
2.50 2.00 1.50
1.00 .50
: N ‘\.
-dim S=l28 ----dim SE 1
:
.
.
.
*. .
.
.
.
.
.
.
.
b
.
.
.
-t.
*.
‘.
‘bI.
*.
-.
‘.
w ‘.
-.
‘b...
Fig. 3. CO (dim. S = 1 and 128). Screening constants 6~ and 50: convergence of the calculations as functions of the selection threshold s.
-1.40 -1.60
-1.80 _
-2.00
"0 2
-2.20
-2.60 -2.80
-3.20
Fig. 4. CO (dim. S * 1 and 128). Susceptibilities fl functions of the selection threshold s.
and X0 (see definition in the text) as
59
Direct comparison with experiment needs introduction rections and will be published at a later date.
of rovibronic cor-
ACKNOWLEDGMENTS
The calculations were performed on the CRAY-1S computer of the Centre de Calcul Vectoriel pour la Recherche (C.C.V.R.). This work was supported by grant No. 3592 from the C.C.V.R. Scientific Council which is gratefully acknowledged. REFERENCES 1 (a) J. A. Pople, J. Chem. Phys., 37 (1962) 53. (b) H. F. Hameka, Mol. Phys., l(l958) 203. (c) F. London, J. Phys. Rad., 8 (1937) 397. (d) R. Ditchfield, Mol. Phys., 27 (1974) 789. (e) E. Vauthier, S. Odiot and F. Tonnard, Can. J. Chem., 60 (1982) 957. (f) H. Fukui, K. Miura, H. Yamasaki and T. Nosaka, J. Chem. Phys., 82 (1985) 1410. 2 W. Kutzelnigg, Isr. J. Chem., 19 (1980) 193; M. Schindler and W. Kutzelnigg, J. Chem. Phys., 76 (1982) 1919. 3 (a) H. Fukui, Int. J. Quantum Chem., 23 (1983) 633; H. Fukui, K. Miura and F. Tada, J. Chem. Phys., 79 (1983) 6112. (b) E. N. Svendren and T. Stroyer-Hansen, Theor. Chim. Acta, 45 (1977) 53. 4M. Iwai and A. Saika, J. Chem. Phys., 77 (1982) 1951. 5 A. J. Sadlej, Mol. Phys., 26 (1973) 1445. 6 D. F. Tuan, T. S. Epstein and J. 0. Hirschfelder, J. Chem. Phys., 44 (1966) 431. 7 A. J. Sadlej, Acta Phys. Pol. A, 49 (1976) 667. 8 B. Levy and J. Ridard, Mol. Phys., 44 (1981) 1099. 9 M. Maestro and R. Moccia, Mol. Phys., 29 (1975) 81. 10 (a) S. T. Epstein, J. Chem. Phys., 42 (1965) 2897. (b) S. T. Epstein, J. Chem. Phys., 58 (1973) 1952. 11 H. F. Hameka, Rev. Mod. Phys., 34 (1962) 87. 12 J. Tillieu, Ann. Phys., 2 (1957) 471, 631; J. Guy and J. Tillieu, J. Chem. Phys., 24 (1956) 1117. 13 M. Karplus and H. J. Kolker, J. Chem. Phys., 38 (1963) 1263. 14U. Sternberg, K. Salzer, H. Pfeifer and W. Haberditzl, Monatsch. Chem., 111 (1980) 505. 15 C. Cohen-Tannoudji, B. Diu and F. LaloE, Mecanique Quantique, I, Hermann, 1973, p. 314. 16 J. Olsen, D. L. Yeager and P. Jorgensen, Adv. Chem. Phys., LIV (1983) 1. 17 W. H. Louisell, Quantum Statistical Theory of Radiation, Wiley, New York, p. 136. 18 P. W. Langhoff, M. Karplus and R. P. Hurst, J. Chem. Phys., 44 (1966) 505. 19 T. P. Das and R. Bersohn, Phys. Rev., 104 (1956) 849. 20 U. Sternberg, W. Haberditzl, Mona&h. Chem., 106 (1975) 701; 109 (1978) 735. 21 U. Sternberg and H. Rosenberger, Chem. Phys., 55 (1981) 275. 22 A. Y. Sadykova and R. M. Aminova, J. Struct. Chem., 21 (1980) 142. 23 A. Okninsky and A. J. Sadlej, Acta Phys. Pol. A, 42 (1972) 709. 24 B. Huron, J. P. Mahieu and P. Rancurel, J. Chem. Phys., 58 (1973) 5745. 25 A. Okninsky and A. J. Sadlej, Mol. Phys., 26 (1973) 1545. 26 J. C. Lie and E. Clementi, J. Chem. Phys., 60 (1974) 1275. 27 A. J. Sadlej and W. T. Raynes, Mol. Phys., 35 (1978) 101. 28 P. Lazzeretti, R. Zanasi and B. Cadioli, J. Chem. Phys., 67 (1977) 382.
60 29 S. Huzinaga and Y. Sakai’, J. Chem. Phys., 50 (1969) 1371. 30 C. Salez and A. Veillard, Theor. Chim. Acta, 11(1968) 441. 31 C. Brezinski, Algorithmes d’accel&ation de la convergence. Etude numbique, Editions Technip., 1978. 32 P. W. Fowler, G. Riley and W. T. Raynes, Mol. Phys., 42 (1981) 1463. 33 J. A. Pople, W. G. Schneider and H. J. Bernstein, High Resolution Nuclear Magnetic Resonance, McGraw-Hill, New York, 59. 34 A. Carrington and A. D. McLachlan, Introduction to Magnetic Resonance, Harper and Row, New York, 1967, p. 55. 35 T. D. Gierke, W. H. Flygare, J. Am. Chem. Sot., 94 (1972) 21.7277. 36 A. Velenik and R. M. Lynder-Bell, Mol. Phys., 19 (1970) 371. 37 T. Weller, W. Meiler, H. J. Kohler, H. Lischka and R. HSller, Chem. Phys. Lett., 98 (1983) 541. 38 R. Jaquet, W. Kutzelnigg and V. Staemmler, Theor. Chim. Acta, 54 (1980) 205.