A numerically stable orbital connection for the calculation of analytical Hessians using perturbation-dependent basis sets

PHYSICS LETTERS ELSEVIER

Chemical Physics Letters 235 (1995) 4 7 - 5 2

A numerically stable orbital connection for the calculation of analytical Hessians using perturbation-dependent basis sets

Kenneth Ruud a, Trygve Helgaker a, Jeppe Olsen b, Poul Jcrgensen c, Keld L. Bak d,1

a Department of Chemistry, University ofOslo, Blindern, N-0315 Oslo, Norway

b Theoretical Chemistry, Chemical Centre, University of Lund, P.O. Box 124, S-221 O0 Lund, Sweden c Department of Chemistry, Aarhus University, DK-8OOO,Arhus C, Denmark

d Department of Physical Chemistry, H.C. Orsteds Institute, Universitetsparken 5, DK-2100 Copenhagen ¢), Denmark Received 9 December 1994

Abstract

It is demonstrated that calculations of second-order molecular properties using perturbation-dependent basis sets are most efficiently carried out using the natural connection.

1. Introduction and only the total magnetizability remains unaffected

by the choice of orbital connection. The relaxation When second-order molecular properties are cal- term is obtained by solving a set of linear response culated from perturbation-dependent basis sets it is equations to some prescribed numerical accuracy.

necessary to define an orbital connection, establish- The most efficient way to solve these equations is ing a one-to-one correspondence between the sets of thus to choose an orbital connection where the relax- orthonormal orbitals at different values of the pertur- ation term is as small as possible. This is done by bation strength x. Molecular properties are indepen- choosing an orbital connection where the relaxation dent of the choice of orbital connection. However, term in the limit of a complete basis set does not different choices shift contributions between the var- contain reorthonormalization contributions. We have ious terms contributing to a given molecular prop- previously described an orbital connection, the natu- erty. This may create numerical problems as well as ral connection, which fulfils this requirement [1].

clutter the physical interpretation of the individual Here we compare the magnitudes of the relaxation contributions. For instance, in calculations of the and static contributions obtained with the natural molecular magnetizability using London atomic or- connection and with the commonly used symmetric bitals the magnitudes of the relaxation and static connection.

terms change when different connections are used, The relaxation term of the magnetizability in- volves the angular momentum operator and if we do not use perturbation-dependent basis sets this term

1 Present address: Uni-C, Olof Palmes All6 38, DK-8200 .Arhus becomes identical to the paramagnetic contribution

N, Denmark. tO the magnetizability. This term also represents the

0 0 0 9 - 2 6 1 4 / 9 5 / $ 0 9 . 5 0 © 1995 Elsevier Science B.V. All rights reserved SSD! 0 0 0 9 - 2 6 1 4 ( 9 5 ) 0 0 0 9 2 - 5

4 8 K. Ruud et al. / Chemical Physics Letters 235 (1995) 47-52

electronic contribution to the molecular g-factor if ble to ~OMO(x0). This can be obtained by minimiz- the gauge origin is chosen as the center of mass. In ing

contrast, for perturbation-dependent basis sets the DnC(x) = ~ II ~/,OMO(X) -- ~/'mOMO(x0) II (6)

relaxation contribution will only correspond to the m '

paramagnetic term if the natural connection is em-

which implies ployed. Any other choice of connection will destroy

the physical interpretation of the relaxation term, as Br/" =

TtW*,

(7) it will then in addition contain spurious reorthonor-

where malization contributions. This is clearly demon-

strated below where we carry out calculations using

Win,,(x)

= (~OMO [~OMO(X0))

the natural and symmetric connections. = (~OMO [ ~ M O ( x 0 ) ) " (8) In order to obtain an explicit expression for the

2. Theoretical background connection matrix T we write

T= W-1M.

(9)

The orbital connection defines a one-to-one corre-

spondence between the so-called orthonormalized Inserting this into Eq. (7) we find that the matrix M must be Hermitian, M = Mr. Inserting Eq. (9) into molecular orbitals (OMOs)~OMO(x)at each strength Eq. (4) and using the Hermiticity of M we obtain x of the perturbation. The OMOs are related to the

the following expression for the connection matrix:

unmodified molecular orbitals (UMOs) through

~OMO(x) =

~_,¢nUMO(x)T,,,,(X).

⁽¹⁾

T=W-I( WS-1W*)I/2

⁽¹⁰⁾

n in the natural connection. In the symmetric connec-

Here tion, we minimize instead the difference function

~OMO(x) =

~,X~(x)C~,m(X),

(2)

DS¢(x)

= E I I ~ ° M ° ( x ) -- t~2MO(X)

II,

(11)

/x m

which gives the symmetric connection matrix

~mUMO(X)

= E X~,( x)C~()m,

(3)

t* T = S - 1 / 2 . (12)

where the

C~,,(x)

are the orbital expansion coeffi- Thus, the natural connection minimizes the overlap cients of the perturbation-dependent atomic basis between the perturbed OMOs and the unperturbed functions

Xt,(x). c(O)-j,,n

are the coefficients at x 0, the UMOs, whereas the symmetric connection mini- reference perturbation strength at which the molecu- mizes the overlap between the UMOs and OMOs at lar property is calculated.

T,m(X)

is the orbital con- a given perturbation strength. This means that for a nection matrix, which satisfies complete basis the OMOs at x 0 will correspond to T * ( x ) S ( x ) T ( x ) = 1, (4) the OMOs at any other x for the natural connection.

Any other connection will introduce a rotation at an where

S(x)

is the UMO overlap matrix. According arbitrary x among the OMOs at x 0 and therefore to Eqs. (1) and (3) the OMOs should reduce to the introduce spurious reorthonormalization contribu- UMOs at x 0, and we therefore require that

T(x o)

⁼ tions in the calculation even in the limit of a com- 1. By rewriting Eq. (4) as plete basis. For a detailed discussion of connections,

see Ref. [1].

(S'/2T)t(S1/2T)

= 1, (5)

we see that S 1 / 2 T is unitary. Thus the number of

undetermined parameters in T equals the number of 3. Numerical examples independent parameters in a unitary matrix. In the

natural connection these parameters are determined In order to demonstrate the numerical superiority by requiring the ~OMO(x) to be as similar as possi- of the natural connection, we have performed a

n u m b e r o f s e l f - c o n s i s t e n t f i e l d ( S C F ) c a l c u l a t i o n s o f c r e a s e , b o t h b y g e n e r a l q u a l i t y ( p V D Z , p V T Z , p V Q Z s o m e s e c o n d - o r d e r p r o p e r t i e s o f the n i t r o g e n a n d p V 5 Z ) a n d b y a u g m e n t a t i o n , the t w o c o n t r i b u - m o l e c u l e u s i n g p e r t u r b a t i o n - d e p e n d e n t b a s i s sets. t i o n s in the s y m m e t r i c c o n n e c t i o n g o to i n f i n i t y T h e m o l e c u l a r H e s s i a n h a s b e e n c a l c u l a t e d w i t h the w h e r e a s n o s u c h e f f e c t is o b s e r v e d for the n a t u r a l a t o m i c o r b i t a l s c l a m p e d to the n u c l e i [2]. W e h a v e c o n n e c t i o n . T h e e f f e c t o f u s i n g the n a t u r a l c o n n e c - t a k e n a d v a n t a g e o f t r a n s l a t i o n a l a n d r o t a t i o n a l in- t i o n is m o r e p r o n o u n c e d f o r m o l e c u l a r H e s s i a n s a n d v a r i a n c e a n d o n l y r e p o r t t h e o n e i n d e p e n d e n t c o m p o - m a g n e t i z a b i l i t i e s t h a n for n u c l e a r s h i e l d i n g s , d u e to n e n t . N u c l e a r s h i e l d i n g s a n d m a g n e t i z a b i l i t i e s are the fact t h a t t h e s e t w o p r o p e r t i e s d e p e n d to s e c o n d c a l c u l a t e d u s i n g L o n d o n a t o m i c o r b i t a l s as d e s c r i b e d o r d e r o n t h e p e r t u r b a t i o n - d e p e n d e n t b a s i s . F o r n u - in Refs. [ 3 - 5 ] . T h e g e o m e t r y c o r r e s p o n d s to the r 0 c l e a r s h i e l d i n g s t h i s d e p e n d e n c y is o n l y linear.

v a l u e f r o m Ref. [6], a n d the b a s i s sets u s e d are the It is o f s o m e i n t e r e s t to c o m p a r e the r e l a x a t i o n c o r r e l a t i o n c o n s i s t e n t s e t s o f D u n n i n g a n d c o - w o r k e r s t e r m s o f t h e t w o c o n n e c t i o n s w i t h the r e s u l t o b t a i n e d [ 7 - 1 0 ] . A l l c a l c u l a t i o n s h a v e b e e n p e r f o r m e d w i t h w i t h o u t the u s e o f L o n d o n o r b i t a l s . U s i n g t h e d - a u g - t h e H E R M I T / S I R I U S / A B A C U S - p r o g r a m s y s t e m c c - p V Q Z b a s i s set, w h i c h s h o u l d b e c l o s e to t h e [ 1 1 - 1 3 ] , in w h i c h w e h a v e i m p l e m e n t e d the n a t u r a l H a r t r e e - F o c k l i m i t f o r b o t h s h i e l d i n g s a n d m a g n e t i - a n d s y m m e t r i c c o n n e c t i o n s , z a b i l i t i e s , w e o b t a i n - 4 6 3 . 3 p p m f o r the s h i e l d i n g F r o m T a b l e s 1 - 3 w e see a s t r i k i n g d i f f e r e n c e r e l a x a t i o n a n d 3 . 8 2 7 2 a u for the m a g n e t i z a b i l i t y b e t w e e n t h e t w o c o n n e c t i o n s . A s t h e b a s i s sets in- r e l a x a t i o n . T h e s e p a r a m a g n e t i c t e r m s are d i r e c t l y

Table 1

Different contributions to the molecular Hessian as well as the total Hessian and vibrational frequency obtained with the symmetric (upper row) and natural connection (lower row). Only the component along the molecular axis is reported. Hessian elements in atomic units, and vibrational frequencies in cm-1

Basis Static Relaxation Total Vibrational/freq.

cc-pVDZ - 84.249396 82.499993 - 1.749403 2569.52

- 81.588101 79.838698 - 1.749403 2569.52

aug-ce-pVDZ - 2051.861708 2050.129022 - 1.732686 2557.22

- 87.563116 85.830430 - 1.732686 2557.22

d-aug-cc-pVDZ - 4305.559664 4303.830477 - 1.729187 2554.63

- 89.020980 87.291785 - 1.729195 2554.64

t-aug-cc-pVDZ - 5462.467546 5460.738922 - 1.728624 2554.22

- 89.068214 87.339563 - 1.728651 2554.24

q-aug-cc-pVDZ - 5788.424342 5786.695793 - 1.728544 2554.16

- 89.080681 87.352108 - 1.728573 2554.18

cc-pVTZ - 2116.640552 2115.038126 - 1.602426 2459.22

- 181.309533 179.707109 - 1.602424 2459.22

aug-cc-pVTZ - 11358.140928 11356.546498 - 1.594430 2453.07

- 194.845845 193.251424 - 1.594421 2453.07

d-aug-cc-pVTZ - 20942.56049 20940.967004 - 1.593486 2452.35

- 196.800884 195.207457 - 1.593427 2452.30

t-aug-cc-pVTZ - 25560.155023 25558.561852 - 1.593171 2452.10

- 197.128018 195.534933 - 1.593085 2452.04

q-aug-cc-pVTZ - 27419.500416 27417.907306 - 1.593110 2452.06

- 197.216327 195.623318 - 1.593009 2451.98

cc-pVQZ - 10821.364790 10819.782077 - 1.582713 2444.04

- 368.075082 366.492382 - 1.582700 2444.03

aug-cc-pVQZ - 39821.782959 39820.201061 - 1.581898 2443.4l

- 378.808857 377.226985 - 1.581872 2443.29

d-aug-cc-pVQZ - 69905.750188 69904.168252 - 1.581936 2443.44

- 383.283770 381.701937 - 1.581833 2443.36

cc-pV5Z - 134786.232291 134784.651125 - 1.581166 2442.85

- 500.538789 498.957717 - 1.581072 2442.78

50 K. Ruud et al. / Chemical Physics Letters 235 (1995) 47-52

related to the spin-rotation and molecular g-factor Table 3

constants, respectively. The results in Tables 2 and 3 Diamagnetic and paramagnetic contributions as well as the total

magnetizability of the nitrogen molecule, obtained using the sym-

clearly demonstrate that only the natural connection metric (upper row) and the natural connection (lower row). All

can be used in order to obtain physically well-de- numbers in atomic units fined paramagnetic contributions.

Basis Static Relaxation Total

In all calculations in Tables 1 - 3 the response (diamag.) (paramag.)

equations have been converged to an accuracy of cc-pVDZ -83.7582 81.3003 -2.4579

l 0 - 4 relative to the norm of the right-hand side. The -5.8493 3.3914 - 2 . 4 5 7 9

fact that we converge to a relative accuracy m e a n s aug-cc-pVDZ - 6 5 5 8 . 9 6 3 3 6556.3638 - 2 . 5 9 9 5

that the numbers of iterations for the solution of the -5.8849 3.2854 - 2 . 5 9 9 5

response equations are about the same for the natural d-aug-cc-pVDZ - 3 4 7 3 1 . 6 1 4 0 34729.0218 - 2 . 5 9 2 3 - 5 . 8 8 1 3 6 3.2894 - 2 . 5 9 4 2

and symmetric connections. For example, when cal-

t-aug-cc-pVDZ - 69566.6695 69564.0914 - 2.5781

culating the magnetic properties with the q-aug-cc- - 5 . 8 8 0 0 3.3058 - 2 . 5 7 4 3

pVTZ basis five iterations were required for the q-aug-cc-pVDZ -91463.8780 91461.3228 -2.5552 symmetric connection and six for the natural connec- -5.8821 3.3338 - 2 . 5 4 8 4

tion. cc-pVTZ - 1 9 7 4 . 4 4 9 6 1971.9118 - 2 . 5 3 7 8

- 6 . 1 2 3 8 3.5861 - 2 . 5 3 7 8 aug-cc-pVTZ -27962.6021 27960.0152 - 2 . 5 8 4 5 - 6 . 1 7 7 8 3.5936 - 2 . 5 8 4 2

Table 2 d-aug-cc-pVTZ - 9 3 6 2 3 . 0 9 8 4 93620.5154 - 2 . 5 8 3 0

Diamagnetic and paramagnetic contributions to the nuclear mag- - 6 . 1 8 3 8 3.6026 - 2 . 5 8 1 3 netic shielding constant as well as the total isotropic chielding

constant of the nitrogen atom, obtained with the symmetric (upper t-aug-cc-pVTZ - 191095.8816 191093.3255 - 2.5561 - 6 . 1 8 2 6 3.6306 - 2 . 5 5 2 0 row) and natural connection (lower row). All numbers in ppm q-aug-cc-pVTZ a -260611.1421 260608.6647 - 2 . 4 7 7 5

Basis Static Relaxation Total - 6.1590 3.6821 - 2.4769

(diam.) (param.) cc-pVQZ - 8862.4173 8859.8590 - 2.5583

cc-pVDZ 359.18 - 4 4 6 . 2 5 - 8 4 . 0 8 - 6 . 3 4 6 1 3.7878 - 2 . 5 5 8 3

326.32 - 4 1 3 . 3 9 - 8 7 . 0 7 aug-cc-pVQZ - 7 5 9 9 8 . 2 1 2 8 75995.6315 - 2 . 5 8 1 3

aug-cc-pVDZ 1368.45 - 1447.62 - 7 9 . 1 7 - 6.3658 3.7851 - 2.5807

359.36 - 4 3 8 . 5 0 - 79.14 d-aug-cc-pVQZ -215674.2268 215671.6464 - 2.5804

d-aug-cc-pVDZ 2721.56 - 2 7 9 9 . 8 7 - 7 8 . 3 2 - 6 . 3 6 6 9 3.7892 - 2 . 5 7 7 7

360.08 - 4 3 8 . 2 6 - 7 8 . 3 7 cc-pV5Z - 1 1 5 1 1 9 . 9 7 9 9 115117.4118 - 2 . 5 6 8 1

t-aug-cc-pVDZ 3419.93 - 3 4 9 7 . 8 7 - 7 7 . 9 3 - 6 . 4 1 0 8 3.8425 - 2 . 5 6 8 3

360.26 - 4 3 8 . 4 4 - 7 8 . 1 8 a Increasing the convergence of the response vector to 10 -5 as q-aug-cc-pVDZ 3703.45 - 3 7 8 1 . 1 5 - 7 7 . 7 0 well as the accuracy of the wave function to 10 -1° in order to 360.23 - 438.50 - 78.27 obtain a magnetizability accurate to four decimal points, we obtain cc-pVTZ 1342.79 - 1448.69 - 105.90 a diamagnetic term of - 260611.1431 and a paramagnetic term of 353.80 - 4 5 9 . 7 1 - 1 0 5 . 9 1 260608.6662. This gives a total isotropic magnetizability of aug-cc-pVTZ 6328.87 - 6 4 3 2 . 9 1 - 1 0 4 . 0 4 - 2 . 4 7 6 9 , identical to the one obtained with the natural connec-

355.37 - 459.45 - 104.08 tion.

d-aug-cc-pVTZ 11618.39 - 11722.80 - 104.41 355.21 - 459.70 - 104.50 t-aug-cc-pVTZ 13659.28 - 13763.48 - 104.20

355.32 - 4 5 9 . 6 8 - 1 0 4 . 3 6 Although the CPU times are much the same, the

q-aug-cc-pVTZ 14456.58 - 1 4 5 6 0 . 4 7 -103.88 absolute accuracies differ significantly in the two 355.43 -459.56 -104.13 calculations. This can be seen from the results ob- cc-pVQZ 2888.41 -2998.47 -110.06 tained for the magnetizability with the q-aug-cc-

354.81 - 464.87 - 110.06

aug-cc-pVQZ 8023.26 - 8132.54 - 109.29 p V T Z basis set and the above convergence threshold.

354.04 - 4 6 3 . 2 3 - 1 0 9 . 1 8 The chosen relative accuracy threshold gives an ab-

d-aug-cc-pVQZ 6620.70 - 6 7 3 0 . 0 4 - 1 0 9 . 3 4 solute accuracy approximately equal to this threshold 354.25 -463.49 - 1 0 9 . 2 4 multiplied by the square root of the paramagnetic

cc-pV5Z 2231.11 - 2 3 4 3 . 9 3 - 1 1 2 . 8 2 term. Thus, in the symmetric connection the solution

353.28 - 466.09 - 112.82

vector is accurate to about 0.05, whereas in the

natural connection it is accurate to 0.0002. As we dicular components of the magnetizability in ben- have taken advantage of the presence of quadratic zene using the symmetric connection, a modest basis errors in the relaxation term [14], the error in the set and a modest active space.

paramagnetic term for the symmetric connection is Some formulations of energy derivatives for per- of the order 0.003, while for the natural connection it turbation-dependent basis sets do not distinguish be- is 4 × 10 -8. Our wave function has been converged tween reorthonormalization and relaxation contribu- to 1 × 10 -8, so the results obtained using the natural tions, see for example Gerratt and Mills [16] and connection are almost as accurate as the wave func- Pulay [17]. In such cases the physical interpretation tion. For the symmetric connection the results cannot of the individual contributions to the properties is be regarded as accurate to more than two decimal cluttered. The numerical problems described above points as also seen from differences in the third may also arise depending on the implicit choice decimal point. For the same cost we have thus made for the orbital connection. For example Gerratt obtained an accuracy of • 10 -8 using the natural and Mills [16] use a symmetric orbital connection for connection and -- 3 × 10 -3 using the symmetric the redundant (occupied-occupied and unoccupied- connection. It is thus more efficient to use the natural unoccupied) orbital parameters of the orbital connec- connection than the symmetric connection or any tion matrix in deriving the coupled Hartree-Fock other connections, as all connections except the natu- equations, and Hoffmann and Schaefer [18] use a ral have spurious reorthonormalization terms enter- Gram-Schmidt orbital connection for these parame-

ing the relaxation term. ters.

The above analysis can also be applied to the To summarize, all perturbation-dependent basis results for the shieldings and the molecular Hessian, set calculations benefit from the natural connection except that the shieldings do not display quadratic both for physical interpretation and for obtaining errors in the relaxation term. Therefore all results for computational efficiency and numerical stability.

the natural connection in Tables 1 - 3 are accurate to the numbers quoted, while this is not the case for the

symmetric connection. We also notice that although Acknowledgement vibrational frequencies are proportional to the square

root of the eigenvalues of the Hessian, differences in This work has received support from the Norwe- vibrational frequencies are observed due to numeri- gian Supercomputing Committee (TRU) through a cal inaccuracies for the symmetric connection when grant of computer time. Partial support by the Danish the basis set becomes bigger (see for instance aug- Natural Science Research Council (Grant No. 11- cc-pVQZ and d-aug-cc-pVQZ). 0924) and Nordisk Forskeruddannelsesakademi is

also gratefully acknowledged.

4. Conclusion

References We conclude that current ab initio methods for

calculating second-order molecular properties using [1] J. Olsen, K.L. Bak, K. Ruud, T. Helgaker and P. J0rgensen, perturbation-dependent basis sets are best carried out Theoret. Chim. Acta., in press.

in the natural connection. If this connection is not [e] T.U. Helgaker, J. AlmliSf, H.J.Aa. Jensen and P. J0rgensen, J. Chem. Phys. 84 (1986) 6266.

used, the results may in some cases become numeri- [3] K. Wolinski, J.F. Hinton and P. Pulay, J. Am. Chem. Soc.

cally unstable. For instance, we have observed that 112 (1990)8251.

the magnetizability in B e l l - is in error by a factor of [4] K. Ruud, T. Helgaker, K.L. Bak, P. JOrgensen and H.J.Aa.

twO for the symmetric connection with a highly Jensen, J. Chem. Phys. 99 (1993)3847.

correlated wave function and large basis sets [15]. [5] K. Ruud, T. Helgaker, R. Kobayashi, P. J~argensen, K.L Bak and H.J.Aa. Jensen, J. Chem. Phys. 100 (1994) 8178.

Instabilities may also occur for large molecules with [6] M.D. Harmony, V.W. Laurie, R.L. Kuczkowski, R.H.

moderate basis sets. We have for instance been Schwendeman, D.A. Ramsey, F.J. Lovas, W.J. Lafferty and unable to obtain identical results for the two perpen- A.G. Maki, J. Chem. Phys. Ref. Data 8 (1979) 619.

52 K. Ruud et aL / Chemical Physics Letters 235 (1995) 47-52

[7] T.H. Dunning Jr., J. Chem. Phys. 90 (1989) 1007. Kobayashi, H. Koch, K. Mikkelsen, J. Olsen, K. Ruud, P.R.

[8] D.E. Woon and T.H. Dunning Jr., J. Chem. Phys. 98 (1993) Taylor and O. Vahtras, ABACUS, a second-order MCSCF

1358. molecular property program.

[9] R.A. Kendall, T.H. Dunning Jr. and R.J. Harrison, J. Chem. [14] H. Sellers, Intern. J. Quantum Chem. 30 (1986) 433.

Phys. 96 (1992) 6796. [15] K. Ruud, T. Helgaker, K.L. Bak, P. J~argensen and J. Olsen,

[10] D.E. Woon and T.H. Dunning Jr., J. Chem. Phys. 100 (1994) submitted for publication.

2975. [16] J. Gerratt and I.M. Mills, J. Chem. Phys. 49 (1968) 1719.

[11] T. Helgaker, P.R. Taylor, K. Ruud, O. Vahtras and H. Koch, [17] P. Pulay, J. Chem. Phys. 78 (1983) 5043.

HERMIT, a molecular integral program. [18] M.R. Hoffmann and H.F. Schaefer III, in: Geometrical [12] H.J.Aa. Jensen and H..~gren, SIRIUS, a program for calcula- derivatives of energy surfaces and molecular properties,

tion of MCSCF wave functions. NATO ASI, eds. P. JCrgensen and J. Simons (Reidel, Dor-

[13] T. Helgaker, K.L. Bak, H.J.Aa. Jensen, P. J~rgensen, R. drecht, 1986)p. 63.