- Email: [email protected]
Generalized valence bond solutions from a constrained coupled cluster method
Chemical Physics ELSEVIER
Chemical Physics 202 (1996) 217-229
Generalized valence bond solutions from a constrained coupled cluster method John Cullen Department of Chemistry, Universityof Manitoba, Winnipeg,Manitoba, CanadaR3T 2N2 Received 31 May 1995
Abstract
The GVB-PP wavefunction is cast into a coupled cluster form with the coupled cluster operator constrained to intrabond double excitations. Following the coupled duster ansatz, where the trial wave function is assumed to satisfy the Schr6dinger equation, projections onto the reference and doubly excited configurations are used to determine the energy and coefficients respectively. A decoupling of these equations results, allowing analytical solutions. The active orbital space is simultaneously optimized to produce the lowest energy. This is carried out efficiently using a procedure previously developed by Head-Gordon and Pople for the direct optimization of a Hartree-Fock wave function. Preliminary calculations for singlet/triplet states and bond dissociations show that although this method is strictly nonvariational, local minima are found which lie within a few tenths of a millihartree of the true GVB-PP energy.
I. Introduction
Multiconfiguration self-consistent field calculations (MCSCF) when coupled with a modest configuration interaction expansion of single and double excitations (CISD) have recently been shown to predict spectroscopic constants remarkably well when extrapolated to the complete basis set limit [1]. In comparison, similar CISD calculations from a Hartree-Fock single reference state produced appreciatively much larger errors even for very small closed shell molecules. The MCSCF accurate results were also found to hold true for the simple multireference generalized valence bond (GVB) method. The wave function for this computational model in it most elemental form is represented by an antisymmetized product of doubly occupied core orbitals, ~bi, singly occupied orbitals of the same spin, q~j, and spin paired geminals, gk,
~OGVB= A[ ~bl ~1 q~2 ~2""" ~°1~°2 ''"
glg2g3""gn]
(la)
Each of the two electron functions, gk, is in turn constructed from local ground and doubly excited configurations of a pair of optimized orbitals
gk = ciA[ ¢~i~i] -~- c;A[ ¢~i*~i* ] 0301-0104/96/$15.00 © 1996 Elsevier Science B.V. All rights reserved SSDI 0301-0104(95)00321-5
(lb)
218
J.
Cullen / Chemical Physics 202 (1996) 217-229
For computational economy, a further constraint of orthogonality between geminals is usually imposed [2]. Known as the perfect pairing model or GVB-PP, resulting calculations without CISD predict the ground and low lying excited state properties and potential energy surfaces semiquantitatively. The set of localized bonding and antibonding natural orbitals, {thi, ~bi*} produced, when decomposed into overlapping hybrids provides a useful interpretive tool for rationalizing stabilities, geometries and reactivities of various closed and open shell species [3] and has been especially useful for transition metal chemistry [4] where Hartree-Fock methods often fail. The GVB method owes much of its continuous development to the work of the Goddard group. Recently a pseudospectral version of the program has been implemented [5] and the DIIS algorithm used in iterating to a converged solutions has been generalized to cases of multiple GVB-PP geminals [6]. The GVB-PP wavefunction in Eq. (1) can also be expanded out into a set of excited configurations from a single reference state of core orbitals and singly or doubly occupied bonding orbitals. When renormalized, the resultant wave function becomes a constrained cluster type A
I~CVB)=er°v~[~o),
(~val~o)
--1,
(2a)
where GEM
TGvB =
E
t i a ~ * a;+* a i a i
(2b)
i q-,
,
The above summation is over the set of GVB-PP geminals, a i , a-,+ denote creation operators for electrons in antibonding orbitals ~bi*, ~* respectively, and a i, a~ are annihilation operators for electrons in bonding orbitals 6i, ~ . If the coupled cluster coefficients, t~, are variationally optimized along with the orbitals the result is identical to the standard GVB-PP solution. However, unlike the ordinary coupled cluster doubles method based on a canonical Hartree-Fock orbital space, this single reference coupled cluster method will perform very well for bond dissociation/formation processes since it is identical to the GVB-PP result. This point was recognized long ago by Ukrainskii in a study of the Peierls instability with the Hubbard model for one-dimensional metals [7]. Despite this, much of the more recent developmental work of coupled cluster theory for multireference problems has focused on either extensions to triple and quadruple excitations from a Hartree-Fock single reference state [8] or developing theories from MCSCF starting points [9]. In this paper instead of using the strictly variational Rayleigh-Ritz procedure to obtain the coupled cluster coefficients as is implicitly done in the ordinary GVB-PP procedure, the more usual coupled cluster ansatz is followed. That is, one assumes the trial wavefunction satisfies the Schrfdinger equation and then projects onto the reference and doubly excited configurations to solve for the energy and coupled cluster coefficients respectively. Our general objective is to examine to what extent a coupled cluster wave function can capture the nondynamic electron correlation when orbitals are properly selected. Here both occupied and unoccupied orbitals within the core + active space are also optimized to produce the lowest energy. Surprisingly, even though this method is nonvariational, results are found to lie within a few tenths of a millihartree of the true GVB-PP energy. A distinctive feature of this approach is that a decoupling of the constrained coupled cluster equations occurs allowing one to solve for the coupled cluster coefficients analytically. A similar decoupling was also noted earlier by Takahashi and Paldus in their one-dimensional model studies of metallic systems with a localized Wannier basis [10]. Here the decoupling allows us to focus our attention exclusively on the orbital optimization step and produce a simpler algorithm for the generation of GVB-PP solutions. For the calculation of single points on the potential surface and energy gradients this feature does not provide a significant advantage over the conventional GVB-PP method. Both treatments spend most of their computational effort within the orbital optimization steps. However, for the efficient determination of vibration frequencies and infrared intensities, when analytical second derivatives are required [11], the standard GVB-PP method will entail the additional calculation of derivatives of the CI coefficients. Although this can be accomplished within the pair excitation coupled perturbed MCSCF procedures developed by Duran et al. [12,13], the resulting
J. Cullen / Chemical Physics 202 (1996) 217-229
219
number of CI coefficient derivatives grows factorially with the number of GVB-PP pairs. In comparison, the development of the analytical second derivatives for the generalized valence bond coupled cluster (GVBCC) method outlined below will require only derivatives of orbital coefficients and integrals. In this initial work we develop the GVBCC method for single point energy calculations using an efficient orbital optimization algorithm which is similar to the procedure developed by Head-Gordon and Pople [14] for the simultaneous optimization of a Hartree-Fock wave function and geometry. Although a Broyden-FletcherGoldfarb-Shanno (BFGS) method [15] is employed, the recursion update procedure of Fischer and Alml6f [16] is also used removing the need to store large matrices such as the Hessian. This treatment does require a two electron integral transformation step at each iteration and when ab initio integrals are generated and transformed in standard fashion a rate limiting computational speed proportional to nN 4 is obtained where n is the number of GVB-PP geminals and N is the basis set size. In comparison, the conventional GVB-PP method avoids integral transformations by using a one- and two-electron density matrix formalism [2]. However, in its place (2n + 1) generalized Fock matrices must be constructed and diagonalized [5]. This results in a scaling of n N 4 for conventional generation of ab initio integrals and nN 25 for pseudospectral formulated methods [17]. Pseudospectrally generated integrals also undergo quicker integral transformations and recently Murphy, Beachy, Friesner and Ringnalda [18] have shown for a localized Moiler Plesset second order perturbation treatment that the integral transformation step executes with a speed proportional to N 275. Since a similar scaling should in principle hold for the localized integral transformation steps within the GVBCC procedure, this method promises to be competitive with conventional GVB-PP and potentially faster for pseudospectral calculations involving a large number of GVB-PP pairs.
2. Theory We begin by expanding out the coupled cluster representation in Eqs. (2), GEM GEM I I~GVB > = I I~0 > "-I- E ti I ~/ii*i "> "~ ½ E E t i t j l i i j
Ibi'i*j°3">Ti i j3 + "",
(3)
where summations are over all GVB-PP geminals (GEM). Substituting (3) back into the Schr6_dinger equation and taking dot products with the reference I~0) and doubly excited configurations I~0~"k" ~ ) produces a modified set of coupled cluster doubles equations GEM
= <~00 IHI qJ0> + ~
ti<~bOIHI
~Jii*i*>
i GEM
=E0+
y' ti(i*ili*i ) i
(4)
= E( ~bo I qJova > = E,
<
I t-/I 0ovB > = < q, 'r
GEM GEM ( 't'k'-k" i'i" k*-k* ib.i*!']i*)_" ) Inlq,0> + E t ix~Yk -k I n l q , , ~ )+½~]~_,titj(~9 * -k IHI T, , , j i i j
=(k'klk*k)
+ t k [ E o - 2 e k + 2e~, + ( k ' k *
Ik'k*)+(kklkk)
GEM
-4(k*k*lkk)+Z(k*klk*k)]
+ Y'~ t k t j ( j * j l j * j ) j~k
k'k*
=E(~Ok 7, I~GVI3)
=Et~.
(5)
J. Cullen/ ChemicalPhysics202 (1996)217-229
220
To eliminate the energy, E, Eq. (4) is substituted into (5). A set of decoupled equations for the coupled cluster coefficients, tk, then results (k*kl k ' k ) - (2(e k - ek. ) + 4 ( k ' k * I kk) - 2 ( k * k l k ' k ) - (kkl kk) - ( k ' k * [ k * k * ) ) t k
-(k*klk*k)t2=O.
(6)
Solving for coefficients and the total energy yields the analytical solutions
tk = l +
2a----k-
1 + 2ak
+1
,
(7a)
GEM
E=E0+7
•
b k + 2 a k - ( b ~ + g a k b k + 8 a ~ ) 1/2
k
where ak = (k*kl k'k),
(7c)
bk=2(8 k.-ek)-a(k*k*lkk)+(kklkk)+(k*k*
Ik*k').
(7d)
For n doubly occupied orbitals, including bonding orbitals of GVB pairs, followed by m singly occupied orbitals of same spin, the uncorrelated part of the energy, E0, is given by Ot~C
O C t OCC
Eo=En,(klhlk)+½Y'.Y'~nknl((kklll)-½(kllkl))-¼ k
k
l
OCC
E k=n+l
OCC
Y'~ (kllkl),
(7e)
l=n+l
where (kl hi k) is a diagonal element of the core Hamiltonian matrix and the summations in the first two terms are over the entire occupied n + m orbital space. Summations in the last exchange term are restricted to the unpaired singly occupied orbitals. Orbital occupancies, n k are defined to be two for k = 1, 2 , . . . , n and one for k = n + 1, n + 2 . . . . . n + m. Corresponding orbital energies ek., s k in Eq. (7d) are simply diagonal elements of the Fock matrix, f, with occ
fkt= ( k l h l l ) + E n r [ ( k l l r r ) - ½(kr I/r)].
(7f)
r
The optimized orbital solutions to Eqs. (7) are found by directly minimizing the energy through unitary transformations from an initial orthonormal reference set. To determine the energy gradient with respect to the rotation angles, ~9i~, of the unitary matrix the chain rule is applied, i.e. dE =p~q
dOij
dE
dCpq
dCpq dOij"
(8)
The coefficient matrix, C, is a unitary matrix where the GVB-PP orbitals thq are constructed from the orthonormal orbitals, Xp which form the initial reference set. To determine the dCpq/d~)ij derivatives the Givens representation of the unitary matrix, C, is used following the Head-Gordon and Pople 'zipper algorithm' [14]. That is
C = HG(~ij),
iy
(9a)
J. Cullen / Chemical Physics 202 (1996) 217-229
221
where
COS Oij
G(e.)
O~j
sin
=
(9b)
Oil
-sin
cos
Oij
The corresponding derivatives of C are then given by dC = U LEFTdG(Oij ) U amrtT, dOij dOi:
(10a)
where '0 0 -
sin
Oij
cos O~j 0
dG(O.) b
(lOb)
d O~j -- COS ~)ij
-- sin O~j
.LE (H rs < ij
(10c,
uRIGHT ----"( U O ( O r s ) ) . \ rs>ij /
(lOd)
Differentiating the energy expression from Eqs. (7) with respect to the coefficients of the orbitals results in dE
dE o
dCpq
dCpq
+
daq(1 4,+,)
dCpo
d Ae q- d f p----~+
Ib2q + 4aqbq + 8aZq
1 db'q ( 2aq+bq ) 2 dCp----qq1 - tb 2 + 4aqbq -I- 8a 2 '
(lla)
222
J.
Cullen/ ChemicalPhysics202 (1996)217-229
where
dAe
d GEM
(
2a~+b, ) dCpq Zr ( 8r" - ~?r) 1 -- ~b2 -[- 4arbr -]- 8 a 2 ,
dCpq
bq = bq -- 2(eq.
-
(llb) (llc)
OOq).
For bonding orbitals, ~bq, the derivatives
dEo/dCpq , dAe/dCpq, daq/dCpq, db'q/dCpq are given by
dE0
(12a)
dCpq = 2nqfpq -t- (nq - 2) r=n+E1( pr l qr), dAB
dCpq
=
(~ ECt~pCuq
PtrA[(/.ty [OrA) -- ½(/zA]tru)])
tzu
( -4fpq 1 -
2aq+bq ~b2 + 4aqbq + 8a2q '
daq = 2(pq* Iqq'), dCpq
--
db'q
(12c)
8(qp I q ' q * ) + 4(qp I qq),
dCpq
(12b)
(12d)
where nq is the orbital occupancy and P is a modified density matrix defined by
(
2a, + b, ) P~A= Y'~r [ C°'r*CAr* C°'rChr] 1 -- CbZ + 4a, b, + SaZ~ " GEM
-
-
(12e)
In the case of antibonding orbitals, ~q., the derivatives dEo/dCpq, , d As/dCpq. both vanish while and
daq//dfpq,, dbq//dfpq, become daq 2(qp] qq*), dCpq,
db.
- -
dCpq .
= 4fq.p - 8 ( q ' p ]
qq) + 4(q*pl q'q* ).
(13a)
(13b)
In addition to forming first derivatives, it is also useful to determine the diagonal matrix elements of the Hessian since the convergence of the BFGS algorithm is found to be substantially improved [14]. As previously noted by Head-Gordon and Pople this can be accomplished by scaling coordinates so that the diagonal second derivatives become unity. For orbital optimizations the scale factors, sij, are defined by
sZ=d2E/dO 2,
(laa)
and the new angular variables, oJij, are formed from
OOij= Sij~)ij.
(14b)
In general, differentiating Eq. (8) results in the complex equation
dfpq dCrs dE d2fpq d~2 =" ~rs Zpq dCrs dfpq dlgi-~j dl~i--~j-t- ~pq dCpq d{~i2 " d2E
d2E
(15)
J. Cullen/ Chemical Physics 202 (1996) 217-229
223
However, since all rotation angles are initially set to zero, the corresponding G(lgij) matrices in the Givens representation become unity while the corresponding first and second derivatives become null matrices with the exception of two off-diagonal or two on-diagonal terms respectively. This results in
dCpq
dOij = t~pq,i) -
6pq,ji,
(16a)
d2Cpq d~)i 2 = - 3pq,i i - 6pq,jj.
(16b)
Eq. (15) then simplifies to dZE
d~)i2
dZE
dC 2
dZE 2
dfijdfji
d2E + --
dffii
dE
dE
dfii
dfjj
(17)
In the case of the GVB-PP coupled cluster model
d2E
da i db i (b 2 + 4aib i+ 8aZi) -3/2
2ai( b i - ai) dCji dCij
dCij dCji
daj dbj . 2\-3/2 + 2aj(bj - at) dfij dfJ i ( b 2 @ 4ajbj + ~aj } 1 d2bj 2a j + bj ) 2ai + bi ) + 1 d2bi ( 1 - - dCji 1 - ~ # + 4a~bj + 8a~ ' 2 dCij dCji ~/bZi + 4aib i + 8a 2 + 2 dCij (18a)
daE0 d2oj(1-
) ÷ d2br1-
2b2( daj
+
J
+
2[ db~ ~2 ~ 2 `.3/2 - Y'.2arl--] (b2+4a~bj+Saj) .
r
2G + br l ~b2 + 4arb r -I- 8a~
~ dCij]
(18b)
A similar expression to (18b) with indices interchanged results for dZE/dC 2. In the case of the cross term (18a) a nonvanishing result occurs only when both orbitals j and i lie in the active space with one being occupied and the other being an antibonding orbital. The corresponding derivatives of b then are db i - = 4(i*i* ] i ' j ) - 4(iili*j) - 2(i*i* l i ' j ) - (ii* l ij), dCi* j
d2bi dCi* j dfji"
= 6 ( i * j l i * j ) - 2(i*i* I jj).
(19a)
(19b)
J. Cullen / Chemical Physics 202 (1996) 217-229
224
Similarly, in the case of the uncorrelated contribution from E 0, produced from Eq. (17), the only nonvanishing rotations occur when ~bj is an occupied core or bonding orbital and ~bi is an antibond or unoccupied orbital, ~bp, then d2Eo -
-
dO~.
d2Eo =
-
dE 0
-
dC2j
dC#
= 2nj( sp - gj) + n213( pj[ pj) - ( pp{ jj)] +(nj-2)
]~
[(pilpi)-(jilji)]+(nj-2)[(ppljj)+(pjlpj)].
(20)
i = n + l
For the remaining correlation energy contribution generated by the GVB-PP method, nonvanishing rotations occur when ~bj is in the active space as either an occupied bond or an unoccupied antibond, thf, and q~i is any other orbital, thp. In the case that ~bj is a bond, the corresponding derivatives of the 'a' and 'b' terms are d2aj
(21a)
dC2j = 2( pj* ] pj* ), d2br
dC2j
[4ep + S ( p p l j * j * ) - 4 ( p p l j j )
-8(pjlpj)]6r,j+
Ar,pj,
(21b)
where Ar,pj vanishes unless p represents an unoccupied orbital, then
Ar,p j = -- [24( pj ] pj ) - 8( pp ] jj ) ] 6~,j + 8( r* r* ] pp ) -- 4 ( r * p l r ' p ) -- 8( rr ] pp) + 4(rp] rp).
(21c)
Similarly for antibond ~bj*
d2aj dC2j '
2( pjl pj),
d2b~ dC2j" = [4ep + 8(pj* [p j * ) - 8(pp[ jj) + 4(pp[ j ' j * )] 6~,j.
(21d)
(21e)
3. Implementation As an initial guess a constrained Hartree-Fock calculation is first performed to produce a set of nonorthogonal strictly localized molecular orbitals (SLMOs) [19]. The SLMOs are formed from atomic orbitals belonging to either one or two atomic centers corresponding respectively to inner shell and lone pair orbitals or two center bonds. Each of the two center bonds is in turn decomposed into a pair of hybrids and an orthogonal antibond is constructed. Once this is done, the original set of occupied SLMOs is symmetrically orthogonalized and its components are projected out of the antibond space. The resulting antibonding orbitals are then symmetrically orthogonalized to complete the construction of the GVB pairs. The remaining orbital space is built by projecting both sets of occupied and antibonding orbitals out from the initial set of atomic orbitals and then carrying out a canonical orthogonalization to remove linear dependencies. In the initial cycle both first and diagonal second derivatives of the energy are formed, the latter being used to set the scale factors for the rotational angles used in the unitary transformations. The energy gradient computed in following cycles is found using the aptly named zipper algorithm of Head-Gordon and Pople [14]. This procedure begins with the old unitary matrix, C, set to U RI°HT and U LEFT set to unity. U Rmr~r is then
J. Cullen/ Chemical Physics 202 (1996) 217-229
225
unzipped by annihilating successive components G(~gij) and replacing it with the corresponding derivative dG(~gij)/d~gij. Multiplying through with U LEFTcompletes the formation of dC/d~gii in Eq. (10a). Combining these results with d E / d C terms produces the gradient, dE/d~ggj, found in Eq. (8). As U RIGHT is unzipped, U LE~T is zipped up by successive multiplications of the G(O~j)'s removed from U RIGHT. For a basis set of N atomic orbitals and an core + active space of m orbitals, the slowest step of the algorithm scales as mN 3. However since only matrix add/multiply operations are involved the procedure executes relatively rapidly in comparison to the other steps within the constrained coupled cluster method. If two electron integrals are generated pseudospectrally then the unmodified zipper algorithm outlined above becomes the rate limiting step. However, because the resulting GVB-PP orbitals are localized, energy gradients are dominated by a small subset of rotation angles linking bonds and antibonds of neighbors and next nearest neighbors. By iteratively rotating angles only within these subspaces a modified algorithm can be obtained which executes in an N 3 step. We have recently implemented such a procedure at the semiempirical level [20]. For standard ab initio two electron integrals, the rate determining step is the transformation of the atomic two electron integrals into the intrabond types (kk I kk), (k * k I k* k), (k * k* I kk), (k * k* I k* k* ) of Eqs. (7c) and (7d) and their corresponding derivatives found in Eqs. (12c), (12d) and (13a), (13b). Formally this scales as m N 4 for conventional generation of two electron integrals and mN 3 for pseudospectral coded methods. As mentioned in the introduction this is as computationally as fast as ordinary GVB-PP methods which instead of performing integral transformations use a generalized Fock matrix formulation requiring different density matrices for each GVB pair. One potential advantage of our approach is the possibility of speeding up the local integral transformations by prescreening integrals with the Schwarz inequality [21]. Optimization of the core + active space orbitals is carried out by minimizing the total energy. This is accomplished using the BFGS algorithm and recursively constructing the dot product of the inverse Hessian with the gradient at each cycle. The latter procedure is due to Fischer and Alml6f [16] who also combined the direct inversion of the iterative subspace (DIIS) technique with the BFGS method. The principal advantage of this approach is that the potentially large inverse Hessian matrix is not stored but instead three small vectors used in the reconstruction are stored at each cycle resulting in a memory requirement of only m(3n + 1) elements for an optimization problem of m degrees of freedom which converges in n cycles.
4. Results of preliminary applications In Table 1 we present comparative SCF and GVB results using a standard 6-31G * * basis for both the singlet and triplet states of the ethylene molecule. All bond lengths were fixed at their equilibrium values and rotations about the CC bond in 30 ° increments were performed. For each method, absolute energies are reported for the planar geometry along with the relative energy increases (in millihartrees) for twisted 30 °, 60° and 90 ° configurations. Singlet/triplet energy splittings, T~'s, are however reported as the differences between the absolute energies of the two states for each geometry. For all six valence bond/antibond pairs found in ethylene both the absolute and relative energies as well as the singlet/triplet splittings show that our generalized valence bond constrained coupled cluster method, which we denote as GVBCC (6PR), to lie within 0.02 millihartrees of the GVB-PP result. For the one pair pi bond/antibond case when both GVB-PP and GVBCC are equivalent, relative energies also closely parallel the six pair results, though at 90 ° there is a difference of 4.2 mH. Relative to the GVB results, the singlet restricted Hartree-Fock (RHF) calculations show a marked deterioration as the broken pi bond geometry of 90 ° is approached. In the triplet state however, the restricted open shell Hartree-Fock (ROHF) calculations are not nearly as bad showing an error of only 4.2 mH at 90 °. A cancellation of errors produces a GVB-PP (1PR) singlet/ROHF triplet splitting that almost reproduces the GVB-PP (6PR) results. In Tables 2 and 3, bond breaking by stretching is examined for a variety of closed and open shell species. Results are reported in the same format as used in Table 1. Literature reference values [22-25] for the complete
J. CuUen / Chemical Physics 202 (1996) 217-229
226 Table 1 Rotational energies of ethylene Method a
0 = 0.0 °
0 = 30.0 °
0 = 60.0 °
0 = 90.0 °
- 78.047240 - 77.917994
19.829 - 11.886
78.568 - 32.729
202.761 - 42.034
97.531
17.949
- 115.549
-78.075154 1.57160
18.141 127.134
69.303 55.128
118.502 -3.376
- 78.142501 - 77.985123 157.378
18.247 - 10.543 128.588
70.381 - 29.309 57.688
122.768 - 37.780 --3.170
SCF IA I 3A 2 T¢ b
129.246
G V B - P P (1PR) IA 1 Te c G V B - P P (6PR) IA i 3A 2 Te b G V B C C (6PR) 1A l
- 78.142521
18.246
70.370
122.781
3A 2 Te b
- 77.985129 157.392
- 10.543 128.603
- 29.308 57.713
- 37.779 -3.169
a Absolute energies are reported for the planar geometry along with the relative energy increases (in millihartrees) for twisted 30 °, 60 ° and 90 ° configurations. b Singlet/triplet energy splittings are reported as the differences between the absolute energies of the two states for each geometry in millihartrees. c GVB(1PR) s i n g l e t / R O H F triplet splitting.
Table 2 Stretching energies a of HF, H 2 0 , and N H 2 Molecule
r = r0
r = 1.5r 0
r = 2.Oro
HF 1X+ SCF MP2 b
- 100.047620 - 100.243698
113.820 93.764
229.378 186.736
G V B - P P (1PR)
- 100.070700
83.752
149.283
G V B - P P (4PR) G V B C C (4PR) CASSCF c
- 100.103008 - 100.103014 - 100.071996
85.037 85.018 82.984
151.120 151.081 149.411
FCI c
- 100.250969
90.576
169.861
H 2 0 1A l SCF MP2 b
- 76.040749 - 76.243660
240.013 195.565
458.118 345.057
GVB-PP (2PR) GVB-PP (4PR) G V B C C (4PR) CASSCF d FCI d
- 76.085365 - 76.106469 - 76.106469 - 76.129876 -- 76.256624
168.455 172.195 172.139 176.735 185.219
270.824 274.666 274.387 289.960 304.355
N H 2 2B 1 SCF MP2 e GVB-PP (3PR) G V B C C (3PR) CASSCF ~ FCI c
- 55.573146 - 55.722105 - 55.604636 - 55.604632 - 55.620751 -55.742620
185.676 156.907 161.797 161.845 122.626 137.411
384.352 305.931 303.844 304.157 209.683 237.096
a Absolute energies are reported for the equilibrium geometry along with the relative energy increases (in millihartrees) for stretched the r = 1.5r 0 and r = 2 . 0 r o configurations. b Ref. [25]. c Ref. [23]. d Ref. [22]. e Computed using G A M E S S , Ref. [26].
J. Cullen / Chemical Physics 202 (1996) 217-229
227
Table 3 Stretching energies a and singlet/triplet splitting b of silylene Molecule
r = r0
r = 1.5r 0
r = 2.0r 0
Sill 2 1A 1 SCF
- 289.997774
142.868
310.374
MP2 c
-290.084770
134.486
267.565
G V B C C (2PR)
- 290.037486
101.042
186.410
G V B - P P (3PR)
- 290.042869
105.631
191.731
G V B C C (3PR)
- 290.042875
105.573
191.337
CASSCF d
- 290.042911
103.129
189.993
FCI d
--290.110207
115.823
202.136
S i l l 2 3B 2 SCF
- 289.990049
150.102
272.429
MP2 c
- 290.063167
135.615
246.330
G V B - P P (2PR)
- 290.010613
137.233
263.269
G V B C C (2PR)
- 290.010613
137.221
263.270
CASSCF d
-- 290.016813
115.502
174.360
FCI d
--290.082313
115.302
181.023
- 30.220
ro SCF
7.725
14.959
MP2
21.603
22.732
0.369
G V B C C (2PR-2PR)
26.873
63.064
103.733
G V B C C (3PR-2PR)
32.262
63.910
104.195
CASSCF
26.098
38.471
10.465
FCI
27.894
27.373
6.781
a Absolute energies are reported for the equilibrium geometry along with the relative energy increases (in millihartrees) for stretched the configurations. b Singlet/triplet energy splittings are reported as the differences between the absolute energies of the two states for each geometry in millihartrees. c Computed using G A M E S S , Ref. [26]. r = 1.5r 0 and r = 2 . 0 r 0
d Ref. [24].
active space self consistent field (CASSCF) method as well as full CI results are also given. The double zeta + polarization basis sets used in these calculations were also employed in our work. However since our programs are implemented within GAMESS [26], six Cartesian d functions were used rather than the normal five. The additional resulting s function though should have only a minor effect (less than a millihartree) on the relative energy curves. Results for closed shell species show once again the SCF energies rapidly deteriorating in bond breaking regions. The M¢ller-Plesset second order (MP2) perturbation calculations partially correct these large errors but cannot overcome the inadequacy of a single reference RHF starting point. The GVB-PP and GVBCC results on the other hand closely parallel one another and are not markedly different from the more computationally intensive CASSCF values. Semiquantitative accuracy with errors of 10-30 mH is found for the relative energies which systematically underestimate the full CI results. Maximum differences between GVB-PP and GVBCC of 0.4-0.1 mH occur at twice the equilibrium bond length. Further calculations at more extended bond lengths show a decrease in this error. For example in the case of H 2 0 at r = 20.0r 0, this error is of the order of a microhartree (see Table 4). For the open shell 2B 1 NH 2 a n d 3 B 2 Sill 2 cases, similar trends as seen in the closed shell results are found for the SCF, MP2 and CASSCF methods. For the singlet/triplet splittings in Sill 2, there is some cancellation of errors with CASSCF results displaying semiquantitative accuracy and the MP2 values showing systematic errors between 5.3 and 6.4 mH. In contrast the GVB-PP and GVBCC calculations overestimate open shell relative
228
J. Cullen / Chemical Physics 202 (1996) 217-229
Table 4 Comparison of total energies and orbital coefficients between GVB-PP and GVBCC Molecule H20 (2.0r 0) (20.0r 0)
GVB-PP energy
GVBCC energy
Coefficient RMS error
# active pairs
- 75.831803 - 75.731796
- 75.832082 - 75.731797
2.90 X 10 4 6.89 x 10 5
4 4
- 290.042869 - 290.010613
- 290.042875 - 290.010613
9.75 X 10- 5 4.53 X 10- 7
3 2
C2H4 09 = 60°)
- 78.072197
- 78.072150
1.63 X 10- 4
6
benzene (-n) (~r + ~cc)
- 230.752709 -230.832772
- 230.752723 -230.832797
4.39 X 10- 5 1.44 X 10-4
3 9
glycine
-283.019065
- 283.019104
6.96 X 10 4
10
Sill 2 (1A1) (3B2)
energies with much larger errors than found in the singlet states where relative energies were underestimated. These errors add together in the resulting s i n g l e t / t r i p l e t splittings producing with the exception of the equilibrium geometry, nonphysical results as seen in Table 3. Carter and Goddard [27] have advocated using a balanced approach in which the singlet nonbonding electrons are allowed to occupy two orbitals just as is the case for the nonbonding electrons in the triplet state. For the perfect pairing model this leads to an energy difference between the GVB-PP (3PR) singlet energy and the G V B - P P (2PR) triplet energy. However our results for S i l l 2 show if the nonbonding pair is not correlated at all (i.e. using a GVB-PP (2PR) singlet energy) the resulting Te is slightly more accurate. In addition to examining energy differences between G V B - P P and G V B C C results we have also compared wave functions. In Table 4 the root mean square differences between the two methods for orbital coefficients of the correlated pairs are presented. Errors are typically of the order 1 0 - 4 - 1 0 -5 indicating the closeness of the wave functions between the two methods.
5. Conclusions The constrained coupled cluster algorithm developed here produces results in excellent agreement with GVB-PP calculations. This demonstrates that the optimization of the orbital space within a single reference coupled cluster procedure is quite capable of picking up the nondynamic electron correlation. This method is competitive with conventional GVB-PP treatment and should execute more rapidly for pseudospectral calculations involving a large number of GVB-PP pairs. The Head-Gordon, Pople algorithm used, also allows within the multireference calculation for the simultaneous optimization of molecular geometry with orbital coefficients. If analytical second derivatives are required then the G V B C C method provides a computationally attractive route since derivatives of CI coefficients are not required.
Acknowledgement W e thank Joe Paldus for providing references to the early work done in this area of research. J.C. gratefully acknowledges the financial support and the very stimulating environment provided by the Goddard group during
J. Cullen / Chemical Physics 202 (1996) 217-229
229
his sabbatical. Special thanks to Jean-Marc Langlois and Richard Muller for their enlightening e-mail discussions on GVB theory.
References [1] D.E. Woon and T.H. Dunning Jr., J. Chem. Phys. 99 (1993) 1914; K.A. Peterson, R.A. Kendall and T.H. Dunning Jr., J. Chem. Phys. 99 (1993) 1920, 9790. [2] F.W. Bobrowicz and W.A. Goddard III, in: Modem theoretical chemistry: methods of electronic structure theory, Vol. 3, ed. H.F. Schaefer III (Plenum Press, New York 1977) p. 79. [3] E.V. Ansyln, M.J. Brusich and W.A. Goddard III, Organometallics 7 (1988) 98. [4] W.A. Goddard III and L.B. Harding, Annu. Rev. Phys. Chem 29 (1978) 363. [5] J.M. Langlois, R.P. Muller, T.R. Coley, W.A. Goddard III, M.N. Ringnalda, R.A. Friesner and W.A. Goddard III, J. Chem. Phys. 92 (1990) 7488. [6] R.P. Muller, J.M. Langlois, M.N. Ringnalda, R.A. Friesner and W.A. Goddard III, J. Chem. Phys. 100 (1994) 1226. [7] I.I. Ukrainskii, Theoret. Math. Phys. 32 (1978) 816; Soviet Phys. JETP 49 (1979) 381; Phys. Stat. Sol. 106 (1981) 55. [8] J.D. Watts and R.J. Bartlett, Intern. J. Quantum Chem. $27 (1993) 51 [9] J. Paldus, in: Methods in computational molecular physics, eds. S. Wilson and G.H.F. Diercksen (Plenum Press, New York, 1992) p. 99. [10] M. Takahashi and J. Paldus, Phys. Rev. B 31 (1984) 5121. [11] Y. Yamaguchi, M.J. Frisch, T.J. Lee, H.F. Schaefer III and J.S. Binkley, Theoret. Chim. Acta 69 (1986) 337. [12] M. Duran, Y. Yamaguchi, R.B. Remington and H.F. Schaefer III, Chem. Phys. 122 (1988) 201. [13] M. Duran, Y. Yamaguchi, R.B. Remington, Y. Osamura and H.F. Schaefer III, J. Chem. Phys. 90 (1989) 334. [14] M. Head-Gordon and J.A. Pople, J. Phys. Chem. 92 (1988) 3063. [15] W.H. Press, B.P. Flannery, S.A. Teukolsky and W.T. Vetterling, in: Numerical recipes: the art of scientific computing (Cambridge Univ. Press, New York, 1986) p. 308. [16] T.H. Fischer and J. Alml6f, J. Phys. Chem. 96 (1992) 9768. [17] R.B. Murphy, R.A. Friesner, M.N. Ringnalda and W.A. Goddard Ili, J. Chem. Phys. 101 (1994) 2986 [18] R.B. Murphy, M.D. Beachy, R.A. Friesner and M.N. Ringnalda, J. Chem. Phys. 103 (1995) 1481. [19] J.M. Cullen, Intern. J. Quantum Chem. $25 (1991) 193. [20] J.M. Cullen, Intern. J. Quantum Chem. 46 (1995) 97. [21] M. Haser and R. Alhrichs, J. Comput. Chem. 10 (1989) 104. [22] C.W. Baushlicher Jr, S.R. Langhoff, P.R. Taylor, N.C. Handy and P.J. Knowles, J. Chem. Phys. 85 (1986) 1469. [23] C.W. Baushlicher Jr. and P.R. Taylor, J. Chem. Phys. 85 (1986) 2779. [24] C.W. Baushchlicher Jr. and P.R. Taylor, J. Chem. Phys. 86 (1987) 1420. [25] S.A. Kucharski, J. Noga and R.J. Bartlett, J. Chem. Phys. 90 (1989) 7282. [26] M.W. Schmidt, K.K. Baldridge, J.A. Boatz, J.H. Jensen, S. Koseki, M.S. Gordon, K.A. Nguyen, T.L. Windus and S.T. Elbert, QCPE Bull. 10 (1990) 52. [27] E.A. Carter and W.A. Goddard Ill, J. Chem. Phys. 86 (1987) 862.
Chemical Physics 202 (1996) 217-229
Generalized valence bond solutions from a constrained coupled cluster method John Cullen Department of Chemistry, Universityof Manitoba, Winnipeg,Manitoba, CanadaR3T 2N2 Received 31 May 1995
Abstract
The GVB-PP wavefunction is cast into a coupled cluster form with the coupled cluster operator constrained to intrabond double excitations. Following the coupled duster ansatz, where the trial wave function is assumed to satisfy the Schr6dinger equation, projections onto the reference and doubly excited configurations are used to determine the energy and coefficients respectively. A decoupling of these equations results, allowing analytical solutions. The active orbital space is simultaneously optimized to produce the lowest energy. This is carried out efficiently using a procedure previously developed by Head-Gordon and Pople for the direct optimization of a Hartree-Fock wave function. Preliminary calculations for singlet/triplet states and bond dissociations show that although this method is strictly nonvariational, local minima are found which lie within a few tenths of a millihartree of the true GVB-PP energy.
I. Introduction
Multiconfiguration self-consistent field calculations (MCSCF) when coupled with a modest configuration interaction expansion of single and double excitations (CISD) have recently been shown to predict spectroscopic constants remarkably well when extrapolated to the complete basis set limit [1]. In comparison, similar CISD calculations from a Hartree-Fock single reference state produced appreciatively much larger errors even for very small closed shell molecules. The MCSCF accurate results were also found to hold true for the simple multireference generalized valence bond (GVB) method. The wave function for this computational model in it most elemental form is represented by an antisymmetized product of doubly occupied core orbitals, ~bi, singly occupied orbitals of the same spin, q~j, and spin paired geminals, gk,
~OGVB= A[ ~bl ~1 q~2 ~2""" ~°1~°2 ''"
glg2g3""gn]
(la)
Each of the two electron functions, gk, is in turn constructed from local ground and doubly excited configurations of a pair of optimized orbitals
gk = ciA[ ¢~i~i] -~- c;A[ ¢~i*~i* ] 0301-0104/96/$15.00 © 1996 Elsevier Science B.V. All rights reserved SSDI 0301-0104(95)00321-5
(lb)
218
J.
Cullen / Chemical Physics 202 (1996) 217-229
For computational economy, a further constraint of orthogonality between geminals is usually imposed [2]. Known as the perfect pairing model or GVB-PP, resulting calculations without CISD predict the ground and low lying excited state properties and potential energy surfaces semiquantitatively. The set of localized bonding and antibonding natural orbitals, {thi, ~bi*} produced, when decomposed into overlapping hybrids provides a useful interpretive tool for rationalizing stabilities, geometries and reactivities of various closed and open shell species [3] and has been especially useful for transition metal chemistry [4] where Hartree-Fock methods often fail. The GVB method owes much of its continuous development to the work of the Goddard group. Recently a pseudospectral version of the program has been implemented [5] and the DIIS algorithm used in iterating to a converged solutions has been generalized to cases of multiple GVB-PP geminals [6]. The GVB-PP wavefunction in Eq. (1) can also be expanded out into a set of excited configurations from a single reference state of core orbitals and singly or doubly occupied bonding orbitals. When renormalized, the resultant wave function becomes a constrained cluster type A
I~CVB)=er°v~[~o),
(~val~o)
--1,
(2a)
where GEM
TGvB =
E
t i a ~ * a;+* a i a i
(2b)
i q-,
,
The above summation is over the set of GVB-PP geminals, a i , a-,+ denote creation operators for electrons in antibonding orbitals ~bi*, ~* respectively, and a i, a~ are annihilation operators for electrons in bonding orbitals 6i, ~ . If the coupled cluster coefficients, t~, are variationally optimized along with the orbitals the result is identical to the standard GVB-PP solution. However, unlike the ordinary coupled cluster doubles method based on a canonical Hartree-Fock orbital space, this single reference coupled cluster method will perform very well for bond dissociation/formation processes since it is identical to the GVB-PP result. This point was recognized long ago by Ukrainskii in a study of the Peierls instability with the Hubbard model for one-dimensional metals [7]. Despite this, much of the more recent developmental work of coupled cluster theory for multireference problems has focused on either extensions to triple and quadruple excitations from a Hartree-Fock single reference state [8] or developing theories from MCSCF starting points [9]. In this paper instead of using the strictly variational Rayleigh-Ritz procedure to obtain the coupled cluster coefficients as is implicitly done in the ordinary GVB-PP procedure, the more usual coupled cluster ansatz is followed. That is, one assumes the trial wavefunction satisfies the Schrfdinger equation and then projects onto the reference and doubly excited configurations to solve for the energy and coupled cluster coefficients respectively. Our general objective is to examine to what extent a coupled cluster wave function can capture the nondynamic electron correlation when orbitals are properly selected. Here both occupied and unoccupied orbitals within the core + active space are also optimized to produce the lowest energy. Surprisingly, even though this method is nonvariational, results are found to lie within a few tenths of a millihartree of the true GVB-PP energy. A distinctive feature of this approach is that a decoupling of the constrained coupled cluster equations occurs allowing one to solve for the coupled cluster coefficients analytically. A similar decoupling was also noted earlier by Takahashi and Paldus in their one-dimensional model studies of metallic systems with a localized Wannier basis [10]. Here the decoupling allows us to focus our attention exclusively on the orbital optimization step and produce a simpler algorithm for the generation of GVB-PP solutions. For the calculation of single points on the potential surface and energy gradients this feature does not provide a significant advantage over the conventional GVB-PP method. Both treatments spend most of their computational effort within the orbital optimization steps. However, for the efficient determination of vibration frequencies and infrared intensities, when analytical second derivatives are required [11], the standard GVB-PP method will entail the additional calculation of derivatives of the CI coefficients. Although this can be accomplished within the pair excitation coupled perturbed MCSCF procedures developed by Duran et al. [12,13], the resulting
J. Cullen / Chemical Physics 202 (1996) 217-229
219
number of CI coefficient derivatives grows factorially with the number of GVB-PP pairs. In comparison, the development of the analytical second derivatives for the generalized valence bond coupled cluster (GVBCC) method outlined below will require only derivatives of orbital coefficients and integrals. In this initial work we develop the GVBCC method for single point energy calculations using an efficient orbital optimization algorithm which is similar to the procedure developed by Head-Gordon and Pople [14] for the simultaneous optimization of a Hartree-Fock wave function and geometry. Although a Broyden-FletcherGoldfarb-Shanno (BFGS) method [15] is employed, the recursion update procedure of Fischer and Alml6f [16] is also used removing the need to store large matrices such as the Hessian. This treatment does require a two electron integral transformation step at each iteration and when ab initio integrals are generated and transformed in standard fashion a rate limiting computational speed proportional to nN 4 is obtained where n is the number of GVB-PP geminals and N is the basis set size. In comparison, the conventional GVB-PP method avoids integral transformations by using a one- and two-electron density matrix formalism [2]. However, in its place (2n + 1) generalized Fock matrices must be constructed and diagonalized [5]. This results in a scaling of n N 4 for conventional generation of ab initio integrals and nN 25 for pseudospectral formulated methods [17]. Pseudospectrally generated integrals also undergo quicker integral transformations and recently Murphy, Beachy, Friesner and Ringnalda [18] have shown for a localized Moiler Plesset second order perturbation treatment that the integral transformation step executes with a speed proportional to N 275. Since a similar scaling should in principle hold for the localized integral transformation steps within the GVBCC procedure, this method promises to be competitive with conventional GVB-PP and potentially faster for pseudospectral calculations involving a large number of GVB-PP pairs.
2. Theory We begin by expanding out the coupled cluster representation in Eqs. (2), GEM GEM I I~GVB > = I I~0 > "-I- E ti I ~/ii*i "> "~ ½ E E t i t j l i i j
Ibi'i*j°3">Ti i j3 + "",
(3)
where summations are over all GVB-PP geminals (GEM). Substituting (3) back into the Schr6_dinger equation and taking dot products with the reference I~0) and doubly excited configurations I~0~"k" ~ ) produces a modified set of coupled cluster doubles equations GEM
ti<~bOIHI
~Jii*i*>
i GEM
=E0+
y' ti(i*ili*i ) i
(4)
= E( ~bo I qJova > = E,
<
I t-/I 0ovB > = < q, 'r
GEM GEM ( 't'k'-k" i'i" k*-k* ib.i*!']i*)_" ) Inlq,0> + E t ix~Yk -k I n l q , , ~ )+½~]~_,titj(~9 * -k IHI T, , , j i i j
=(k'klk*k)
+ t k [ E o - 2 e k + 2e~, + ( k ' k *
Ik'k*)+(kklkk)
GEM
-4(k*k*lkk)+Z(k*klk*k)]
+ Y'~ t k t j ( j * j l j * j ) j~k
k'k*
=E(~Ok 7, I~GVI3)
=Et~.
(5)
J. Cullen/ ChemicalPhysics202 (1996)217-229
220
To eliminate the energy, E, Eq. (4) is substituted into (5). A set of decoupled equations for the coupled cluster coefficients, tk, then results (k*kl k ' k ) - (2(e k - ek. ) + 4 ( k ' k * I kk) - 2 ( k * k l k ' k ) - (kkl kk) - ( k ' k * [ k * k * ) ) t k
-(k*klk*k)t2=O.
(6)
Solving for coefficients and the total energy yields the analytical solutions
tk = l +
2a----k-
1 + 2ak
+1
,
(7a)
GEM
E=E0+7
•
b k + 2 a k - ( b ~ + g a k b k + 8 a ~ ) 1/2
k
where ak = (k*kl k'k),
(7c)
bk=2(8 k.-ek)-a(k*k*lkk)+(kklkk)+(k*k*
Ik*k').
(7d)
For n doubly occupied orbitals, including bonding orbitals of GVB pairs, followed by m singly occupied orbitals of same spin, the uncorrelated part of the energy, E0, is given by Ot~C
O C t OCC
Eo=En,(klhlk)+½Y'.Y'~nknl((kklll)-½(kllkl))-¼ k
k
l
OCC
E k=n+l
OCC
Y'~ (kllkl),
(7e)
l=n+l
where (kl hi k) is a diagonal element of the core Hamiltonian matrix and the summations in the first two terms are over the entire occupied n + m orbital space. Summations in the last exchange term are restricted to the unpaired singly occupied orbitals. Orbital occupancies, n k are defined to be two for k = 1, 2 , . . . , n and one for k = n + 1, n + 2 . . . . . n + m. Corresponding orbital energies ek., s k in Eq. (7d) are simply diagonal elements of the Fock matrix, f, with occ
fkt= ( k l h l l ) + E n r [ ( k l l r r ) - ½(kr I/r)].
(7f)
r
The optimized orbital solutions to Eqs. (7) are found by directly minimizing the energy through unitary transformations from an initial orthonormal reference set. To determine the energy gradient with respect to the rotation angles, ~9i~, of the unitary matrix the chain rule is applied, i.e. dE =p~q
dOij
dE
dCpq
dCpq dOij"
(8)
The coefficient matrix, C, is a unitary matrix where the GVB-PP orbitals thq are constructed from the orthonormal orbitals, Xp which form the initial reference set. To determine the dCpq/d~)ij derivatives the Givens representation of the unitary matrix, C, is used following the Head-Gordon and Pople 'zipper algorithm' [14]. That is
C = HG(~ij),
iy
(9a)
J. Cullen / Chemical Physics 202 (1996) 217-229
221
where
COS Oij
G(e.)
O~j
sin
=
(9b)
Oil
-sin
cos
Oij
The corresponding derivatives of C are then given by dC = U LEFTdG(Oij ) U amrtT, dOij dOi:
(10a)
where '0 0 -
sin
Oij
cos O~j 0
dG(O.) b
(lOb)
d O~j -- COS ~)ij
-- sin O~j
.LE (H rs < ij
(10c,
uRIGHT ----"( U O ( O r s ) ) . \ rs>ij /
(lOd)
Differentiating the energy expression from Eqs. (7) with respect to the coefficients of the orbitals results in dE
dE o
dCpq
dCpq
+
daq(1 4,+,)
dCpo
d Ae q- d f p----~+
Ib2q + 4aqbq + 8aZq
1 db'q ( 2aq+bq ) 2 dCp----qq1 - tb 2 + 4aqbq -I- 8a 2 '
(lla)
222
J.
Cullen/ ChemicalPhysics202 (1996)217-229
where
dAe
d GEM
(
2a~+b, ) dCpq Zr ( 8r" - ~?r) 1 -- ~b2 -[- 4arbr -]- 8 a 2 ,
dCpq
bq = bq -- 2(eq.
-
(llb) (llc)
OOq).
For bonding orbitals, ~bq, the derivatives
dEo/dCpq , dAe/dCpq, daq/dCpq, db'q/dCpq are given by
dE0
(12a)
dCpq = 2nqfpq -t- (nq - 2) r=n+E1( pr l qr), dAB
dCpq
=
(~ ECt~pCuq
PtrA[(/.ty [OrA) -- ½(/zA]tru)])
tzu
( -4fpq 1 -
2aq+bq ~b2 + 4aqbq + 8a2q '
daq = 2(pq* Iqq'), dCpq
--
db'q
(12c)
8(qp I q ' q * ) + 4(qp I qq),
dCpq
(12b)
(12d)
where nq is the orbital occupancy and P is a modified density matrix defined by
(
2a, + b, ) P~A= Y'~r [ C°'r*CAr* C°'rChr] 1 -- CbZ + 4a, b, + SaZ~ " GEM
-
-
(12e)
In the case of antibonding orbitals, ~q., the derivatives dEo/dCpq, , d As/dCpq. both vanish while and
daq//dfpq,, dbq//dfpq, become daq 2(qp] qq*), dCpq,
db.
- -
dCpq .
= 4fq.p - 8 ( q ' p ]
qq) + 4(q*pl q'q* ).
(13a)
(13b)
In addition to forming first derivatives, it is also useful to determine the diagonal matrix elements of the Hessian since the convergence of the BFGS algorithm is found to be substantially improved [14]. As previously noted by Head-Gordon and Pople this can be accomplished by scaling coordinates so that the diagonal second derivatives become unity. For orbital optimizations the scale factors, sij, are defined by
sZ=d2E/dO 2,
(laa)
and the new angular variables, oJij, are formed from
OOij= Sij~)ij.
(14b)
In general, differentiating Eq. (8) results in the complex equation
dfpq dCrs dE d2fpq d~2 =" ~rs Zpq dCrs dfpq dlgi-~j dl~i--~j-t- ~pq dCpq d{~i2 " d2E
d2E
(15)
J. Cullen/ Chemical Physics 202 (1996) 217-229
223
However, since all rotation angles are initially set to zero, the corresponding G(lgij) matrices in the Givens representation become unity while the corresponding first and second derivatives become null matrices with the exception of two off-diagonal or two on-diagonal terms respectively. This results in
dCpq
dOij = t~pq,i) -
6pq,ji,
(16a)
d2Cpq d~)i 2 = - 3pq,i i - 6pq,jj.
(16b)
Eq. (15) then simplifies to dZE
d~)i2
dZE
dC 2
dZE 2
dfijdfji
d2E + --
dffii
dE
dE
dfii
dfjj
(17)
In the case of the GVB-PP coupled cluster model
d2E
da i db i (b 2 + 4aib i+ 8aZi) -3/2
2ai( b i - ai) dCji dCij
dCij dCji
daj dbj . 2\-3/2 + 2aj(bj - at) dfij dfJ i ( b 2 @ 4ajbj + ~aj } 1 d2bj 2a j + bj ) 2ai + bi ) + 1 d2bi ( 1 - - dCji 1 - ~ # + 4a~bj + 8a~ ' 2 dCij dCji ~/bZi + 4aib i + 8a 2 + 2 dCij (18a)
daE0 d2oj(1-
) ÷ d2br1-
2b2( daj
+
J
+
2[ db~ ~2 ~ 2 `.3/2 - Y'.2arl--] (b2+4a~bj+Saj) .
r
2G + br l ~b2 + 4arb r -I- 8a~
~ dCij]
(18b)
A similar expression to (18b) with indices interchanged results for dZE/dC 2. In the case of the cross term (18a) a nonvanishing result occurs only when both orbitals j and i lie in the active space with one being occupied and the other being an antibonding orbital. The corresponding derivatives of b then are db i - = 4(i*i* ] i ' j ) - 4(iili*j) - 2(i*i* l i ' j ) - (ii* l ij), dCi* j
d2bi dCi* j dfji"
= 6 ( i * j l i * j ) - 2(i*i* I jj).
(19a)
(19b)
J. Cullen / Chemical Physics 202 (1996) 217-229
224
Similarly, in the case of the uncorrelated contribution from E 0, produced from Eq. (17), the only nonvanishing rotations occur when ~bj is an occupied core or bonding orbital and ~bi is an antibond or unoccupied orbital, ~bp, then d2Eo -
-
dO~.
d2Eo =
-
dE 0
-
dC2j
dC#
= 2nj( sp - gj) + n213( pj[ pj) - ( pp{ jj)] +(nj-2)
]~
[(pilpi)-(jilji)]+(nj-2)[(ppljj)+(pjlpj)].
(20)
i = n + l
For the remaining correlation energy contribution generated by the GVB-PP method, nonvanishing rotations occur when ~bj is in the active space as either an occupied bond or an unoccupied antibond, thf, and q~i is any other orbital, thp. In the case that ~bj is a bond, the corresponding derivatives of the 'a' and 'b' terms are d2aj
(21a)
dC2j = 2( pj* ] pj* ), d2br
dC2j
[4ep + S ( p p l j * j * ) - 4 ( p p l j j )
-8(pjlpj)]6r,j+
Ar,pj,
(21b)
where Ar,pj vanishes unless p represents an unoccupied orbital, then
Ar,p j = -- [24( pj ] pj ) - 8( pp ] jj ) ] 6~,j + 8( r* r* ] pp ) -- 4 ( r * p l r ' p ) -- 8( rr ] pp) + 4(rp] rp).
(21c)
Similarly for antibond ~bj*
d2aj dC2j '
2( pjl pj),
d2b~ dC2j" = [4ep + 8(pj* [p j * ) - 8(pp[ jj) + 4(pp[ j ' j * )] 6~,j.
(21d)
(21e)
3. Implementation As an initial guess a constrained Hartree-Fock calculation is first performed to produce a set of nonorthogonal strictly localized molecular orbitals (SLMOs) [19]. The SLMOs are formed from atomic orbitals belonging to either one or two atomic centers corresponding respectively to inner shell and lone pair orbitals or two center bonds. Each of the two center bonds is in turn decomposed into a pair of hybrids and an orthogonal antibond is constructed. Once this is done, the original set of occupied SLMOs is symmetrically orthogonalized and its components are projected out of the antibond space. The resulting antibonding orbitals are then symmetrically orthogonalized to complete the construction of the GVB pairs. The remaining orbital space is built by projecting both sets of occupied and antibonding orbitals out from the initial set of atomic orbitals and then carrying out a canonical orthogonalization to remove linear dependencies. In the initial cycle both first and diagonal second derivatives of the energy are formed, the latter being used to set the scale factors for the rotational angles used in the unitary transformations. The energy gradient computed in following cycles is found using the aptly named zipper algorithm of Head-Gordon and Pople [14]. This procedure begins with the old unitary matrix, C, set to U RI°HT and U LEFT set to unity. U Rmr~r is then
J. Cullen/ Chemical Physics 202 (1996) 217-229
225
unzipped by annihilating successive components G(~gij) and replacing it with the corresponding derivative dG(~gij)/d~gij. Multiplying through with U LEFTcompletes the formation of dC/d~gii in Eq. (10a). Combining these results with d E / d C terms produces the gradient, dE/d~ggj, found in Eq. (8). As U RIGHT is unzipped, U LE~T is zipped up by successive multiplications of the G(O~j)'s removed from U RIGHT. For a basis set of N atomic orbitals and an core + active space of m orbitals, the slowest step of the algorithm scales as mN 3. However since only matrix add/multiply operations are involved the procedure executes relatively rapidly in comparison to the other steps within the constrained coupled cluster method. If two electron integrals are generated pseudospectrally then the unmodified zipper algorithm outlined above becomes the rate limiting step. However, because the resulting GVB-PP orbitals are localized, energy gradients are dominated by a small subset of rotation angles linking bonds and antibonds of neighbors and next nearest neighbors. By iteratively rotating angles only within these subspaces a modified algorithm can be obtained which executes in an N 3 step. We have recently implemented such a procedure at the semiempirical level [20]. For standard ab initio two electron integrals, the rate determining step is the transformation of the atomic two electron integrals into the intrabond types (kk I kk), (k * k I k* k), (k * k* I kk), (k * k* I k* k* ) of Eqs. (7c) and (7d) and their corresponding derivatives found in Eqs. (12c), (12d) and (13a), (13b). Formally this scales as m N 4 for conventional generation of two electron integrals and mN 3 for pseudospectral coded methods. As mentioned in the introduction this is as computationally as fast as ordinary GVB-PP methods which instead of performing integral transformations use a generalized Fock matrix formulation requiring different density matrices for each GVB pair. One potential advantage of our approach is the possibility of speeding up the local integral transformations by prescreening integrals with the Schwarz inequality [21]. Optimization of the core + active space orbitals is carried out by minimizing the total energy. This is accomplished using the BFGS algorithm and recursively constructing the dot product of the inverse Hessian with the gradient at each cycle. The latter procedure is due to Fischer and Alml6f [16] who also combined the direct inversion of the iterative subspace (DIIS) technique with the BFGS method. The principal advantage of this approach is that the potentially large inverse Hessian matrix is not stored but instead three small vectors used in the reconstruction are stored at each cycle resulting in a memory requirement of only m(3n + 1) elements for an optimization problem of m degrees of freedom which converges in n cycles.
4. Results of preliminary applications In Table 1 we present comparative SCF and GVB results using a standard 6-31G * * basis for both the singlet and triplet states of the ethylene molecule. All bond lengths were fixed at their equilibrium values and rotations about the CC bond in 30 ° increments were performed. For each method, absolute energies are reported for the planar geometry along with the relative energy increases (in millihartrees) for twisted 30 °, 60° and 90 ° configurations. Singlet/triplet energy splittings, T~'s, are however reported as the differences between the absolute energies of the two states for each geometry. For all six valence bond/antibond pairs found in ethylene both the absolute and relative energies as well as the singlet/triplet splittings show that our generalized valence bond constrained coupled cluster method, which we denote as GVBCC (6PR), to lie within 0.02 millihartrees of the GVB-PP result. For the one pair pi bond/antibond case when both GVB-PP and GVBCC are equivalent, relative energies also closely parallel the six pair results, though at 90 ° there is a difference of 4.2 mH. Relative to the GVB results, the singlet restricted Hartree-Fock (RHF) calculations show a marked deterioration as the broken pi bond geometry of 90 ° is approached. In the triplet state however, the restricted open shell Hartree-Fock (ROHF) calculations are not nearly as bad showing an error of only 4.2 mH at 90 °. A cancellation of errors produces a GVB-PP (1PR) singlet/ROHF triplet splitting that almost reproduces the GVB-PP (6PR) results. In Tables 2 and 3, bond breaking by stretching is examined for a variety of closed and open shell species. Results are reported in the same format as used in Table 1. Literature reference values [22-25] for the complete
J. CuUen / Chemical Physics 202 (1996) 217-229
226 Table 1 Rotational energies of ethylene Method a
0 = 0.0 °
0 = 30.0 °
0 = 60.0 °
0 = 90.0 °
- 78.047240 - 77.917994
19.829 - 11.886
78.568 - 32.729
202.761 - 42.034
97.531
17.949
- 115.549
-78.075154 1.57160
18.141 127.134
69.303 55.128
118.502 -3.376
- 78.142501 - 77.985123 157.378
18.247 - 10.543 128.588
70.381 - 29.309 57.688
122.768 - 37.780 --3.170
SCF IA I 3A 2 T¢ b
129.246
G V B - P P (1PR) IA 1 Te c G V B - P P (6PR) IA i 3A 2 Te b G V B C C (6PR) 1A l
- 78.142521
18.246
70.370
122.781
3A 2 Te b
- 77.985129 157.392
- 10.543 128.603
- 29.308 57.713
- 37.779 -3.169
a Absolute energies are reported for the planar geometry along with the relative energy increases (in millihartrees) for twisted 30 °, 60 ° and 90 ° configurations. b Singlet/triplet energy splittings are reported as the differences between the absolute energies of the two states for each geometry in millihartrees. c GVB(1PR) s i n g l e t / R O H F triplet splitting.
Table 2 Stretching energies a of HF, H 2 0 , and N H 2 Molecule
r = r0
r = 1.5r 0
r = 2.Oro
HF 1X+ SCF MP2 b
- 100.047620 - 100.243698
113.820 93.764
229.378 186.736
G V B - P P (1PR)
- 100.070700
83.752
149.283
G V B - P P (4PR) G V B C C (4PR) CASSCF c
- 100.103008 - 100.103014 - 100.071996
85.037 85.018 82.984
151.120 151.081 149.411
FCI c
- 100.250969
90.576
169.861
H 2 0 1A l SCF MP2 b
- 76.040749 - 76.243660
240.013 195.565
458.118 345.057
GVB-PP (2PR) GVB-PP (4PR) G V B C C (4PR) CASSCF d FCI d
- 76.085365 - 76.106469 - 76.106469 - 76.129876 -- 76.256624
168.455 172.195 172.139 176.735 185.219
270.824 274.666 274.387 289.960 304.355
N H 2 2B 1 SCF MP2 e GVB-PP (3PR) G V B C C (3PR) CASSCF ~ FCI c
- 55.573146 - 55.722105 - 55.604636 - 55.604632 - 55.620751 -55.742620
185.676 156.907 161.797 161.845 122.626 137.411
384.352 305.931 303.844 304.157 209.683 237.096
a Absolute energies are reported for the equilibrium geometry along with the relative energy increases (in millihartrees) for stretched the r = 1.5r 0 and r = 2 . 0 r o configurations. b Ref. [25]. c Ref. [23]. d Ref. [22]. e Computed using G A M E S S , Ref. [26].
J. Cullen / Chemical Physics 202 (1996) 217-229
227
Table 3 Stretching energies a and singlet/triplet splitting b of silylene Molecule
r = r0
r = 1.5r 0
r = 2.0r 0
Sill 2 1A 1 SCF
- 289.997774
142.868
310.374
MP2 c
-290.084770
134.486
267.565
G V B C C (2PR)
- 290.037486
101.042
186.410
G V B - P P (3PR)
- 290.042869
105.631
191.731
G V B C C (3PR)
- 290.042875
105.573
191.337
CASSCF d
- 290.042911
103.129
189.993
FCI d
--290.110207
115.823
202.136
S i l l 2 3B 2 SCF
- 289.990049
150.102
272.429
MP2 c
- 290.063167
135.615
246.330
G V B - P P (2PR)
- 290.010613
137.233
263.269
G V B C C (2PR)
- 290.010613
137.221
263.270
CASSCF d
-- 290.016813
115.502
174.360
FCI d
--290.082313
115.302
181.023
- 30.220
ro SCF
7.725
14.959
MP2
21.603
22.732
0.369
G V B C C (2PR-2PR)
26.873
63.064
103.733
G V B C C (3PR-2PR)
32.262
63.910
104.195
CASSCF
26.098
38.471
10.465
FCI
27.894
27.373
6.781
a Absolute energies are reported for the equilibrium geometry along with the relative energy increases (in millihartrees) for stretched the configurations. b Singlet/triplet energy splittings are reported as the differences between the absolute energies of the two states for each geometry in millihartrees. c Computed using G A M E S S , Ref. [26]. r = 1.5r 0 and r = 2 . 0 r 0
d Ref. [24].
active space self consistent field (CASSCF) method as well as full CI results are also given. The double zeta + polarization basis sets used in these calculations were also employed in our work. However since our programs are implemented within GAMESS [26], six Cartesian d functions were used rather than the normal five. The additional resulting s function though should have only a minor effect (less than a millihartree) on the relative energy curves. Results for closed shell species show once again the SCF energies rapidly deteriorating in bond breaking regions. The M¢ller-Plesset second order (MP2) perturbation calculations partially correct these large errors but cannot overcome the inadequacy of a single reference RHF starting point. The GVB-PP and GVBCC results on the other hand closely parallel one another and are not markedly different from the more computationally intensive CASSCF values. Semiquantitative accuracy with errors of 10-30 mH is found for the relative energies which systematically underestimate the full CI results. Maximum differences between GVB-PP and GVBCC of 0.4-0.1 mH occur at twice the equilibrium bond length. Further calculations at more extended bond lengths show a decrease in this error. For example in the case of H 2 0 at r = 20.0r 0, this error is of the order of a microhartree (see Table 4). For the open shell 2B 1 NH 2 a n d 3 B 2 Sill 2 cases, similar trends as seen in the closed shell results are found for the SCF, MP2 and CASSCF methods. For the singlet/triplet splittings in Sill 2, there is some cancellation of errors with CASSCF results displaying semiquantitative accuracy and the MP2 values showing systematic errors between 5.3 and 6.4 mH. In contrast the GVB-PP and GVBCC calculations overestimate open shell relative
228
J. Cullen / Chemical Physics 202 (1996) 217-229
Table 4 Comparison of total energies and orbital coefficients between GVB-PP and GVBCC Molecule H20 (2.0r 0) (20.0r 0)
GVB-PP energy
GVBCC energy
Coefficient RMS error
# active pairs
- 75.831803 - 75.731796
- 75.832082 - 75.731797
2.90 X 10 4 6.89 x 10 5
4 4
- 290.042869 - 290.010613
- 290.042875 - 290.010613
9.75 X 10- 5 4.53 X 10- 7
3 2
C2H4 09 = 60°)
- 78.072197
- 78.072150
1.63 X 10- 4
6
benzene (-n) (~r + ~cc)
- 230.752709 -230.832772
- 230.752723 -230.832797
4.39 X 10- 5 1.44 X 10-4
3 9
glycine
-283.019065
- 283.019104
6.96 X 10 4
10
Sill 2 (1A1) (3B2)
energies with much larger errors than found in the singlet states where relative energies were underestimated. These errors add together in the resulting s i n g l e t / t r i p l e t splittings producing with the exception of the equilibrium geometry, nonphysical results as seen in Table 3. Carter and Goddard [27] have advocated using a balanced approach in which the singlet nonbonding electrons are allowed to occupy two orbitals just as is the case for the nonbonding electrons in the triplet state. For the perfect pairing model this leads to an energy difference between the GVB-PP (3PR) singlet energy and the G V B - P P (2PR) triplet energy. However our results for S i l l 2 show if the nonbonding pair is not correlated at all (i.e. using a GVB-PP (2PR) singlet energy) the resulting Te is slightly more accurate. In addition to examining energy differences between G V B - P P and G V B C C results we have also compared wave functions. In Table 4 the root mean square differences between the two methods for orbital coefficients of the correlated pairs are presented. Errors are typically of the order 1 0 - 4 - 1 0 -5 indicating the closeness of the wave functions between the two methods.
5. Conclusions The constrained coupled cluster algorithm developed here produces results in excellent agreement with GVB-PP calculations. This demonstrates that the optimization of the orbital space within a single reference coupled cluster procedure is quite capable of picking up the nondynamic electron correlation. This method is competitive with conventional GVB-PP treatment and should execute more rapidly for pseudospectral calculations involving a large number of GVB-PP pairs. The Head-Gordon, Pople algorithm used, also allows within the multireference calculation for the simultaneous optimization of molecular geometry with orbital coefficients. If analytical second derivatives are required then the G V B C C method provides a computationally attractive route since derivatives of CI coefficients are not required.
Acknowledgement W e thank Joe Paldus for providing references to the early work done in this area of research. J.C. gratefully acknowledges the financial support and the very stimulating environment provided by the Goddard group during
J. Cullen / Chemical Physics 202 (1996) 217-229
229
his sabbatical. Special thanks to Jean-Marc Langlois and Richard Muller for their enlightening e-mail discussions on GVB theory.
References [1] D.E. Woon and T.H. Dunning Jr., J. Chem. Phys. 99 (1993) 1914; K.A. Peterson, R.A. Kendall and T.H. Dunning Jr., J. Chem. Phys. 99 (1993) 1920, 9790. [2] F.W. Bobrowicz and W.A. Goddard III, in: Modem theoretical chemistry: methods of electronic structure theory, Vol. 3, ed. H.F. Schaefer III (Plenum Press, New York 1977) p. 79. [3] E.V. Ansyln, M.J. Brusich and W.A. Goddard III, Organometallics 7 (1988) 98. [4] W.A. Goddard III and L.B. Harding, Annu. Rev. Phys. Chem 29 (1978) 363. [5] J.M. Langlois, R.P. Muller, T.R. Coley, W.A. Goddard III, M.N. Ringnalda, R.A. Friesner and W.A. Goddard III, J. Chem. Phys. 92 (1990) 7488. [6] R.P. Muller, J.M. Langlois, M.N. Ringnalda, R.A. Friesner and W.A. Goddard III, J. Chem. Phys. 100 (1994) 1226. [7] I.I. Ukrainskii, Theoret. Math. Phys. 32 (1978) 816; Soviet Phys. JETP 49 (1979) 381; Phys. Stat. Sol. 106 (1981) 55. [8] J.D. Watts and R.J. Bartlett, Intern. J. Quantum Chem. $27 (1993) 51 [9] J. Paldus, in: Methods in computational molecular physics, eds. S. Wilson and G.H.F. Diercksen (Plenum Press, New York, 1992) p. 99. [10] M. Takahashi and J. Paldus, Phys. Rev. B 31 (1984) 5121. [11] Y. Yamaguchi, M.J. Frisch, T.J. Lee, H.F. Schaefer III and J.S. Binkley, Theoret. Chim. Acta 69 (1986) 337. [12] M. Duran, Y. Yamaguchi, R.B. Remington and H.F. Schaefer III, Chem. Phys. 122 (1988) 201. [13] M. Duran, Y. Yamaguchi, R.B. Remington, Y. Osamura and H.F. Schaefer III, J. Chem. Phys. 90 (1989) 334. [14] M. Head-Gordon and J.A. Pople, J. Phys. Chem. 92 (1988) 3063. [15] W.H. Press, B.P. Flannery, S.A. Teukolsky and W.T. Vetterling, in: Numerical recipes: the art of scientific computing (Cambridge Univ. Press, New York, 1986) p. 308. [16] T.H. Fischer and J. Alml6f, J. Phys. Chem. 96 (1992) 9768. [17] R.B. Murphy, R.A. Friesner, M.N. Ringnalda and W.A. Goddard Ili, J. Chem. Phys. 101 (1994) 2986 [18] R.B. Murphy, M.D. Beachy, R.A. Friesner and M.N. Ringnalda, J. Chem. Phys. 103 (1995) 1481. [19] J.M. Cullen, Intern. J. Quantum Chem. $25 (1991) 193. [20] J.M. Cullen, Intern. J. Quantum Chem. 46 (1995) 97. [21] M. Haser and R. Alhrichs, J. Comput. Chem. 10 (1989) 104. [22] C.W. Baushlicher Jr, S.R. Langhoff, P.R. Taylor, N.C. Handy and P.J. Knowles, J. Chem. Phys. 85 (1986) 1469. [23] C.W. Baushlicher Jr. and P.R. Taylor, J. Chem. Phys. 85 (1986) 2779. [24] C.W. Baushchlicher Jr. and P.R. Taylor, J. Chem. Phys. 86 (1987) 1420. [25] S.A. Kucharski, J. Noga and R.J. Bartlett, J. Chem. Phys. 90 (1989) 7282. [26] M.W. Schmidt, K.K. Baldridge, J.A. Boatz, J.H. Jensen, S. Koseki, M.S. Gordon, K.A. Nguyen, T.L. Windus and S.T. Elbert, QCPE Bull. 10 (1990) 52. [27] E.A. Carter and W.A. Goddard Ill, J. Chem. Phys. 86 (1987) 862.