WALKING ON MCSCF POTENTIAL ENERGY SURFACES ...

4 downloads 0 Views 3MB Size Report
Hans. Aa. Jensen and Poul. Department of. Chemistry, Aarhus University, DK-8000 .... d' h d,-l a new step 1S eva uate W1t trust ra 1US a. Simons et a1. [2] suggest the ... cond-order energy change for mode j to be positive and the predicted.
WALKING ON MCSCF POTENTIAL ENERGY SURFACES: APPLICATION TO H20 2 AND NH3

Danny L. Yeager, Chemistry Department, Texas A & M University College Station, TX 77843, USA Hans Aa. Jensen and Poul Department of Chemistry, Aarhus University, DK-8000 Arhus C, Denmark Trygve U. Helgaker, Department of Chemistry, University of Oslo, Blindern, N-03l5, Oslo 3, Norway ABSTRACT. We develop a surface walking procedure based on a restrictedstep algorithm combined with the Newton-Raphson method for MCSCF potential 'energy surfaces. We discuss some of the details of the algorithm including how to optimally use it to obtain a first-order estimate of the MCSCF wave function at a new point. We apply the algorithm to the problems of the rotation barrier and dissociation in H20 2 (minimum basis set) and the dissociation in NH3 (double zeta plus polar1zation basis). 1.

INTRODUCTION

An important task in ab-initio quantum chemistry is to reliably determine equilibrium and transition-state geometries [1-6]. Such geometries are stationary points on a Born-Oppenheimer potential energy surface and are therefore characterized by having a zero molecular gradient. Equilibrium geometries are usually efficiently determined with algorithms which use only molecular gradient information. However, transition states are more difficult to determine as a result of which algorithms for determining transition states usually require both the molecular gradient and the molecular Hessian to be evaluated. Recent methods for finding transition states [1-3] have emphasized the importance of staying close to a "stream bed". These walks proceed by maximizing the total energy in one direction and minimizing it in all others, firmly maintaining some step-length control. The step length control is to ensure a controlled movement along a reaction coordinate until the Hessian has the desired structure (number of negative eigenvalues) and until the Newton-Raphson step can be taken with confidence. The restricted-step algorithm [7] combined with the Newton-Raphson algorithm can easily be adapted for walking on energy hypersurfaces. For walks to energy minima this procedure guarantees convergence [7,8]. To our knowledge, it is the only second-order procedure in use in quantum chemistry which has been proven mathematically to have this feature. (Demonstrating simple energy decrease on each iteration locally is usually not sufficient to show guaranteed convergence). For walks to other 229

P. J¢rgensen and J. Simons (eds.), Geometrical Derivatives of Energy Surfaces and Molecular Properties, 229-241. © 1986 by D. Reidel Publishing Company.

230

D. L. YEAGER ET AL.

stationary points it has repeatedly been shown to be rapid and reliable [2, 9-10].

The restricted step algorithm was first used in quantum chemistry for MCSCF wavefunction optimization by et al. [8]. Jensen and [11] and Werner and Knowles [12] subsequently used the procedure in large-scale second-order MCSCF optimization methods. Golab et al.[lO] used the method to demonstrate the existence of many "non-physical" stationary points on a MCSCF energy hypersurface. Simons et al. [2] adapted and modified the procedure for walking on a geometrical-surface and applied it to model surfaces and to two potential energy surfaces of HCN. The restricted-step algorithm is based on defining a hypersphere "trust region" with a trust radius h inside which a quadratic approximation to the total energy is a good representation of the exact energy function. Starting from some initial value, hI' which determines the maximum allowable step-length norm for the in1tial iteration, subsequent maximum radii, hk ' increase, decrease, or stay the same based on how well the second-order Taylor series expansion at a point with a step agrees with the exact value of the function at !. If the step-length is smaller than the trust radius and the Hessian has the correct structure (desired number of negative eigenvalues) the Newton-Raphson step is taken. This occurs in the local region of a stationary point. Consequently, the procedure has second-order convergence. Thus the restricted-step algorithm allows the maximum extraction of information from a second-order procedure. If a second-order expansion for a given step predicts a total energy which agrees fairly well with the exact value subsequently found at that point, then the second-order expansion is a good approximation. If the second-order expansion value does not agree very well with the exact value then the step is rejected and a new, reduced step is taken. For large-scale MCSCF optimization or for MCSCF geometry walks, step-rejection may not be inexpensive, however, it is not nearly as costly as taking a step that is too long or in the wrong direction. Such steps may lead to regions on the energy hypersurface which are far from the desired stationary point and thus result in poor convergence. 2.

THEORY

2.1.

Second-Order Expansions and the Restricted-Step Algorithm

At a point on a multidimensional energy surface the energy can be expanded in a Taylor series to obtain the energy at X (1)

where /:::'X = X - X

--

-

.::0

(2)

WALKING ON MCSCF POTENTIAL ENERGY SURFACES

231

F.

(3)

H ••

(4)

1

1J

The infinite series in Eq. (1) can be truncated at second-order: (5)

As a measure of how well Eq. (5) approximates Eq. (1) the ratio of Eq. (1) to Eq. (5) in the k'th iteration can be examined b.kX rk

= kx-

_ kx

(6)

=0

E(k!) - E(k!o) E(2) (k!)

E(k!o)

1 + R(k!. kX )

(7)

=0

where R(k!. kX ) =0

E(kX) _ E(2)(kX) E (2) (k!)

(8)

E(k!o)

Thus if thi second-oiderkexpansion is a good approximation to the exact energy at ! then R(!, !o) is small • . The more the ratio r k deviates from one the more important are the rema1nder terms. Fletcher defines a hypersphere("trust region") with radius ("trust radius") in which the quadratic approximation to the energy resembles the exact energy function. If the Hessian has the correct structure (i.e. the correct number of negative eigenvalues) and the Newton-Raphson step-length norm lies within the trust-region then the Newton-Raphson step is taken. Otherwise, depending on r k a modified is taken on the boundary of the hypersphere or the step is rejected and a new, smaller trust radius is defined. Fletcher [7] has suggested a set of rules to follow for determining for minimization. For walks to transition states where positive and negative energy contributions may cancel, Simons et al. [2] suggest: or 2-r min ()

'1l

::0

c

en

...::

Cl

::0

zttl

ttl

r

>-l

z ;;

>-l ttl

0

'""

'1l

()

()

s::

z

0

Cl

Z

:>::

> r

238

D. L. YEAGER ET AL.

a bond length of 1.042 a.u. and a bond angle of 101.9 0 • The experimental bond length and bond angle is l6 1.024±0.005 and 103.3±0.5 . The details of the walk are given in Table 2 and show it is almost a valley walk. For example, in step 15 the step in the lowest mode is 0.376 a.u. and the trust radius is 0.38 a.u. At step 29 it was clear that the walk resulted in the dissociation products NH +H and the algorithm was forced to take the Newton-Raphson step in the consequtive iterations. As the transition state is at infinity the NewtonRaphson steps as expected did gradually increase. The increase of the Newton-Raphson step reflects that the molecular Hessian goes toward zero with one higher inverse power in the internuclear distance than the molecular gradient. In H20 2 we used a sub-minimum basis set. The exponents and contraction coefficients are given in Table 3. TABLE 3. Subminimum basis set used for H20 2 • Atom

o

Function s

s

p

H

s

Exponent

Contrac don coefficient

420.0

0.032265

91.9805

0.156410

24.4515

0.447813

7.2230

0.481602

1.0631

0.504708

0.3227

0.616743

1.0000

0.520000

0.36503

0.604484

4.0000

0.100000

0.6535

0.334495

0.1776

0.674700

There were seven inactive and four active orbitals g1v1ng twenty configuration state functions. The stationary points obtained from surface walking along the "gentlest" (smallest eigenvalue) three valleys, star-

-149.12211

+3 (H 2O+O)

1.171

1.184

1.187

1.194 (0.950)

02-H2

1.986

1.580

1.575

1.557 (1.475)

01-02

37.54

102.53

96.42

100.16 (94.8)

Hl-01-02

95.22

102.54

96.42

100.16 (94.8)

H2-02-01

Bond Angles (0)

"0 1---{) 2" H2

H1

Olson and P.C Cross

89.51

0.00

180.00

89.89 (119.8)

H1-01-02-H2

(0)

Dihedral Angle

b) The trans barrier is calculated to 1.62 kca1/mo1 and the cis barrier to 5.78 kca1/mo1. The experimental values are 1.10 kca1/mo1 and 7.57 kca1/mo1 respectively. (C.S. Ewig and D.O. Harris, J. Chem. Phys. 52, 6268 (1970).

a) Numbers in parenthesis are the experimental values reported byRL Redington, J. Chem. Phys. 1311 (1962).

1.520

1.184

-149.19458b )

+1 (cis)

1.194 (0.950) 1.187

-149.20379

01-H1

Bond Distances (!)

-149.20121 b )

a)

Total Energy(a.u.)

-1 (trans)

Ground state

Mode

Transition States Obtained in H20 2 by Walking up the Gentlest Modes

TABLE 4

:;::

IV W

'-0