it it it

(1)

WALKING ON MCSCF POTENTIAL ENERGY SURFACES: APPLICATION TO H202 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 restricted- step 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 H202 (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 sta-

tes 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.

(2)

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-phy- sical" 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. Conse- quently, 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, re- duced 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 ⁽²⁾

(3)

WALKING ON MCSCF POTENTIAL ENERGY SURFACES

F. 1

H ••

1J

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

(3)

(4)

(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

=

^{kx _ kx}

- =0 (6)

r k E(k!) - E(k!o)

R(k!. kX ) 1 +

E(2) (k!) E(k!o) ⁼⁰ ⁽⁷⁾

where

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

=0 E (2) (k!) E(k!o)

(8) 231

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:

(1) if r min <rk<r goo d or 2-r goo d <r<2-r. , then m1n = ; (2) if r good <r<2-r good' then hk+l ⁼ ;

(3) if rk<rmin or 2-r min <rk ' then the step is rejected and

(4)

, 1 d ' h d , - l

a new step 1S eva uate W1t trust ra 1US a

Simons et a1. [2] suggest the parameter values rgood = 0.85, rmin =0.7, and a =-r.S-. We have found through numerical experience on small systems that a more optimal choice may be rgood = 0.90, rmin = 0.35, a = 1.3, and hl = 0.01. Near convergence or when one or more Hessian eigenvalues are small even more stringent criteria are used:

r goo d = 0.95, r, = 0.75, and m1n a = 1.3 • 2.2. The Optimal Step Length and Direction

Let us assume that the local quadratic energy function and the current trust radius h are known. The optimal point is then a stationary point inside or on the boundary of the hypersphere. The Newton-Raphson step is optimal if it is inside the trust region and if the Hessian has the desired index (number of negative eigenvalues). If the Newton-Raphson step is longer than the trust radius or if the Hessian does not have the desired index then the optimal step is on the boundary of the trust region. The optimal point on the boundary may be determined by looking for stationary points of the Lagrangian function

(9) where V is a Lagrange multiplier introduced to keep the step of lengh h • We now consider how to determine V • Setting = 0 determines

the optimal step on the boundary

= -(M -

V!) -1

!

(10)

Examination of = h2 reveals there are several possible choi- ces for V which givea step length of h • Using a representation in which the Hessian is diagonal

H U

= = U

= -

E: ⁽¹¹⁾

E:, ,

1J <5" E:,

1J 1 (12)

we may express the step-length constraint equation as '" 2 -2

= L: F, (E:, - v)

1 1 (13)

where! and;;X are the molecular gradient and+the step vector, trans- formed to the diagonal representation (using Q) • In order to determine the appropriate value of V it is also useful to examine the predicted first- and second-order changes in energy [2]

(5)

WALKING ON MCSCF POTENTIAL ENERGY SURFACES 233

t,E(1)

Q9

E(l) - E(!o) L

.

...., 2 F. (v - E.) -1 ⁽¹⁴⁾ t,E(2)

Q9

^E(2)(!) _{- E(X )}₌₀

=

L _L

.

^{...., 2}F. ^{(V -} E.) ^-2^(V- IE.) ^L ⁽¹⁵⁾ Assume now that the Hessian eigenvalues at X are ordered

El E2 •••• In walking toward a minimum first-order and the second-order energy change for each mode should be minimized. This can only be fulfilled if

(16) (note that E may be negative). This gives a unique choice for v such that the step length is h.

In walking from a minimum to a transition state along a streambed we would like to maximize the energy in one direction (leading eventua- lly to a negative Ej at the stationary point) and to minimize the energy in all other directions (leading to positive E.'S, i j at the stationary point). We therefore require both the preaicted first- and second-order energy change for mode j to be positive and the predicted first- and second-order energy change to be negative for all other modes. From Eqs. (14) and (15) we see that V should therefore be chosen such that

V > !E., and V > E.

J - J ⁽¹⁷⁾

and that

V < E. and V < IE. ; i j

- L - L (18)

Eq. (18) cannot always be fulfilled within the coordinate system described thus far, particularly if we wish to walk along modes corresponding to higher Hessian eigenvalues. In order to ensure that Eq. (18) is fulfilled, a scaling (i.e., stretching or compression) of the coordinates can be introduced

[;X.

=

n. [;X!

- I . - L (19)

Straightforward reexpression of the molecular gradient and Hessian in the scaled coordinate system gives

[;X. _n.AX....

'

L L L

F.

_L

2 ...

n. _LF • _L V - n.

z

E.

L L

F.

_L

I

V. - E. _L _L

(20)

(21)

(6)

That is, by scaling the individual coordinates by n. one effectively scales the Hessian eigenvalues by n? • Thus individGal modes can be scaled so that Eq. (18) is fulfillea (i.e. different level shifts v!

may be chosen for each of the modes). Alternatively, if there are

ral near-zero modes or a walk is along a mode near a bifurcation (nearly equal eigenvalues) Eq. (21) can be used to give more desired characte- ristics (e.g., to scale low Hessian eigenvalues upward). In practice, usually only one mode needs to be scaled. If streambed walking along the El rnodelis desired and Eq. (18) is not fulfilled then we have chosen

ni =

4

E2/El ' and all other n. = 1 • If streambedlwalking along the Ej mode is desired (j F 1) then have chosen =

4

El/E. and all

other ni = 1 . J J

2.3. Imposition of Geometrical Constraints on Potential Energy Surface Walking Procedures

Recently Taylor and Simons [13] proposed a procedure to apply geometrical constraints to multidimensional surface walking algorithms. In par- ticular they discussed center-of-mass displacements, infinitesimal rotations, planes of symmetry, bond lengths, bond angles, substituent group internal rotations, and dihedral angels. These constraints are useful in eliminating center-of-mass translations and rotations, as well as other

"uninteresting" degrees of freedom. They proposed the use of several linearized constraints

(22) and several corresponding Lagrange multipliers

n ..

This leads to a set of coupled inhomogeneous equations to determine fhe Lagrange multipliers for the constraints which must be solved for each value of v in a line search for the step-length constrained to h.

In our opinion it is easier to simply use projectors to eliminate unwanted motions, i.e.

t,x +

!:

t,X (23)

where

P 1 - m E iv.><v.i i=l 1 1

(24)

The sum in Eq. (24) goes over all constraints that are imposed on the walk, i.e. all constraints given in Eq. (22). Undesired motions are eleminated by using the projected second-order energy

(25) where

!:

is the projected gradient, and

!: !:

is the projected Hessian.

(7)

WALKING ON MCSCF POTENTIAL ENERGY SURFACES 235

2.4. Walks With a MCSCF Wave Function

When evaluating the MCSCFmolecular Hessian, the linear response of the MCSCF wave function to a nuclear displacement must be evaluated [14] • In a preceeding paper, Jensen showed how this linear response may also be used to obtain a first-order estimate of the MCSCF wave function at the displaced geometry + • Denoting the linear response of the MCSCF wavefunction to a nuclear distortion a by t a , Jensen showed that to first order the MCSCF wave function at + is determined by the parameter set

L: ta

a - -a ⁽²⁶⁾

where the vector A contains two sets of parameters, one describing an orbital rotation and the other the change in the configuration state function amplitudes. The total energy of this state may be used in the walking algorithm to give a first-order estimate of the exact MCSCF energy at + • If an unacceptable ratio, see Eq. (7), is found with this fLrst-order energy replacing the exact energy the step should be rejected and the trust radius updated according to the walking algorithm. When an unacceptable ratio is determined in this way, it simply reflects the fact that the size of the nuclear step is so large that non-quadratic terms dominate in the total energy expansion.

The exact MCSCF energy must also fulfill the conditions on r k • For instance,near an avoided crossing, the first-order estimate may give a good ratio while the exact MCSCF wave function changes a lot and gives a poor ratio r k thus in this case the step should also be rejected. Use of the first-order estimate for the MCSCF wave function is a necessary but not a sufficient condition; however, the use of this necessary condition can save a substantial amount of computer time.

3. APPLICATIONS TO H202 AND NH3

We applied the technique described in the previous sections to geometrical surface walking on the electronic ground state of hydrogen peroxide, H202 ' and of ammonia, •

In NH3 we used a douBle-zeta contracted Gaussian basis set [15]

with a set of six d polarization functions (exponent 0.80) on N and a set of p polarization functions (exponent 1.00) on each H. The hydrogen s exponents were scaled to the equivalent of an effective Slater exponent of 1.2. A complete active space of all valence orbitals (2s, 2p on oxygen and Is on hydrogen) was chosen (8 electrons in 7 orbitals) re- sulting in 490 configuration state functions. The optimized MCSCF ground state geometry is given in Table 1 together with the experimental equilibrium geometry.

(8)

236

TABLE 1

Stationary Point Obtained in NH3•

Ground statea Transition state

Total

Energy(a.u.) N-Hl

-56.28179 1.028 (1.019)

Bond Distances (A) N-H2

1.028 (1.019) -56.11949 1.042

(1.024=0.005 ) 1.042

(1.024=0.005) N-H3

1.028 (1.019) 3.457

D. L. YEAGER ET AL.

Bond Angle(o) Hl-N-H2

105.0 (109.0 ) 101.9 (103.3:0.5)

a) Numbers in parenthesis are the experimental values reported by O.

Bastiansen and B. Beagley, Acta Chem. Scan. 18 2077(1964).

b) Numbers in parenthesis are the experimental values for the B state 2 of NH2 reported by K. Dressler and D.A. Ramsay, Phil. Trans. loy.

Soc., London Ser. A 251 553 (1959).

When a transition state search is initiated from a minimum, the number of streambed directions is twice the number of internal Hessian modes as for each mode we may take the first step in opposite directions.

We arbitrarily designate one of these two directions by the positive "+"

and the other by the native "-".

The walk which was initiated by taking a step of 0.1 a.u. in the negative eigen-direction of the lowest Hessian eigenvalue at the equilibrium geometry resulted in a sequence of iterations which approached the geometry of the inversion barrier. When very close to the inversion barrier, the step sizes had to be constrained to be 0.015 a.u., thus leading to very slow convergence toward the inversion barrier. This very slow convergence reflects the fact that the chosen active space is inad- equate to describe the inversion barrier because there are no active orbitals to correlate the lone pair TI orbital on nitrogen at the planar geometry. For this reason, we have not reported the rotation barrier in NH)"

In the walk which was initiated in the positive eigendirection of the lowest Hessian eigenvalue (initially pushing the hydrogens together) the surface is fairly quadratic and, consequently, the trust radius in- creases to 0.38 a.u. (starting from 0.10 a.u.) prior to the appearance of a negative Hessian eigenvalue. The walk resulted in the dissociation products H + NH2 with no barrier, that is the transition state appear at infinite separation. In Table 1 we report the geometry of the dissociation product as it occurs when the leaving hydr2gen atom is 3.457 a.u.

from the nitrogen atom. The NH2 radical is in the Bl ground state with

(9)

TABLE 2 D:tai1s of the surface walks at selected points in NH3 walking up the gentlest mode away from inver- S10n. Ratioa) Step Trust Lowest Total number radius r Hessian energy Bond distances (A) Bond Angle (0) (a. u.) (see text) eigenvalues (a. u.) (a.u.) N-H1 N-H2 N-H3 H1-N-H2 0 1.27E-2 -56.28179 1.028 1.028 1.028 105.0 4.22E-2 8.90E-2

15 0.38 1.18 1.85E-2 -56.22025 1.051 1.099 1.227 106.5 (0.99)

4.75E-2 8.10E-2

21 0.38 1.26 -2.08E-2 -56.13109 1.039 1.071 1.834 102.7 (0.94 )

2.14E-2 4.01E-2

23 0.19 0.90 -2.63E-2 -56.12756 1.031 1.066 1.987 103.5 (1.01)

2.23E-2 3.13E-2

32 0.71 1.81 -2.27E-3 -56.11949 1.042 1.042 3.457 101.9 (1.16) 2.42E-3 4.22E-3

a) The ratio in parenthesis is obtained using the exact MCSCF total energy and the other ratio is obtained using the first order estimate of the MCSCF total energy.

::IE > r :>:: Z Cl 0 z s:: () en () '1l , '"" 0 >-l ttl z >-l ;; r ttl z ttl ::0 Cl ...:: en c ::0 '1l > () ttl en N w -.l

(10)

a bond length of 1.042 a.u. and a bond angle of 101.90 • 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 conse- qutive iterations. As the transition state is at infinity the Newton- Raphson 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 H202 we used a sub-minimum basis set. The exponents and contrac- tion coefficients are given in Table 3.

TABLE 3. Subminimum basis set used for H202 • Exponent Contrac don

coefficient Atom Function

o

s 420.0 0.032265

91.9805 0.156410 24.4515 0.447813 7.2230 0.481602

s 1.0631 0.504708

0.3227 0.616743

p 1.0000 0.520000

0.36503 0.604484

H s 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-

(11)

TABLE 4 H1 Transition States Obtained in H202 by Walking up the Gentlest Modes "0 1---{) 2" H2 Mode Total Energy(a.u.) Bond Distances (!) Bond Angles (0) Dihedral Angle (0) 01-H1 02-H2 01-02 Hl-01-02 H2-02-01 H1-01-02-H2 Ground state a) -149.20379 1.194 1.194 1.557 100.16 100.16 89.89 (0.950) (0.950) (1.475) (94.8) (94.8) (119.8) -1 (trans) -149.20121b) 1.187 1.187 1.575 96.42 96.42 180.00 +1 (cis) -149.19458b) 1.184 1.184 1.580 102.53 102.54 0.00 +3 (H2O+O) -149.12211 1.520 1.171 1.986 37.54 95.22 89.51 a) Numbers in parenthesis are the experimental values reported byRL Redington, Olson and P.C Cross J. Chem. Phys. 1311 (1962). 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).

:;:: ;.. l"""' Z C'l 0 z ;;:: n <Zl n 'Tl "I:i 0 ..., t"l Z ..., ;; l"""' t"l Z t"l ;<l C'l ><: <Zl c::: ;<l 'Tl ;.. n t"l <Zl IV W '-0

(12)

ting from the MCSCF optimized ground state geometry will be discussed.

The walks initiated along the two directions of the lowest Hessian mode led to the cis and trans rotational barriers, respectively,and the transition state geometries are reported in Table 4. The walk along the positive second eigendirection (+2) was facing a wall and after two steps the lowest Hessian eigenvalue had increased by a factor 50 as a result of which this walk was abandoned. The walk along the negative second eigendirection (-2) resulted in the dissociation products HOO + H ,

(transition state at infinite separation) with no barrier. When walking up the positive third eigendirection (+3) a very interesting transition state was obtained. This transition state consisted of a pseudo water molecule loosely bound to the other oxygen atom. The 0-0 bond length was 1.986 a.u. and the 0-0 bond was orthogonal to the plane of the pseudo water molecule. The transition state geometry is given in Table 4.

Although the above H20 calculations are too small to reliably re- produce the physics of the molecule, they do demonstrate the effec- tiveness of the surface walklng algorithm for finding transition states.

4. SUMMARY AND CONCLUSION

We have presented a surface walking algorithm and have used it for walking up the lowest eigenvalue valleys on MCSCF potential energy surfaces in H202 and NH3 • These calculations were ferformed entirely on a VAX 11/780 computer using the SIRIUS¹⁷/ABACUS 8 MCSCF program system. Each point (MCSCF + MC Hessian) required ^tV3 hours for NH3 • The MCSCF optimization and the evaluation of the molecular Hessian each require approxi- mately the same amount of computer time. Thus, although this surface walking to locate transition states is time consuming on a VAX 11/780

for somewhat larger systems, this kind of "real" chemistry for small mo- lecules is very feasible on faster computers, i.e. supercomputers.

ACKNOWLEDGEMENTS

Danny L. Yeager would like to acknowledge research support from the National Science Foundation (Grant CHE-84l3442) and the Robert A. Welch Foundation (Grant A-770).

REFERENCES

1. C.J. Cerjan and W. Miller, J. Chern. Phys. 2800 (1981).

2. J. Simons, P. J¢rgensen, H. Taylor, and J. Ozment, J. Phys. Chern.

2745 (1983).

3. A. Banerjee, N. Adams, J. Simons, and R. Shepard, J. Phys. Chern.

89, 52 (1985).

(13)

WALKING ON MCSCF POTENTIAL ENERGY SURF ACES

4. S.M. Crippen and H.A. Scheraga, Arch. Biochern. Biophys., 144, 462 (1971) •

5. S. Bell, J.S. Cirghton, and R. Fletcher, Chern. Phys. Lett., 82,

241

122 (1981); S. Bell and J. Cirghton, J. Chern. Phys. 80, 2464--(1984).

6. See, for example: J.W. McIver and A. Komornicki, J. Amer. Chern.

Soc., 94, 2625 (1972); A. Komornicki, K. Oshida, K. Morokuma, R.

Ditchfie1d, and M. Conrad, Chern. Phys. Lett., 595 (1977).

7. R. Fletcher, 'Practical Methods of Optimization', John Wiley, New York, 1980.

8. P. J¢rgensen, P. Swanstr¢m, and D.L. Yeager, J. Chern. Phys., 1§., 347 (1983).

9. J. Golab, D.L. P. J¢rgensen, Chern. Phys., 1§., 175 (1983).

10. J. Golab, D.L. Yeager, and P. J¢rgensen, Chern. Phys. 83 (1985).

11. H.J. Aa. Jensen and P. J¢rgensen, J. Chern. Phys. 80, 1204 (1984).

12. H.-J. Werner and P. Knowles, J. Chern. Phys. 5055 (1985).

13. H. Taylor and J. Simons, J. Phys. Chern. 89, 684 (1985).

14. P. J¢rgensen and J. Simons, Second Quantization Based Methods in Quantum Chemistry (Academic, New York, 1981).

15. T.H. Dunning, J. Chern. Phys., 53 2823 (1970).

16. K. Dressler and D.A. Ramsay, Phil. Trans. Ray. Soc., London Ser.

A 251 553 (1959).

17. H.J. Aa. Jensen and H. Agren, Chern. Phys. Lett., 110, 2 (1984).

18. T.U. He1gaker, J. A1m1of, H.J. Aa. Jensen and P. J¢rgensen, "Mole- cular Hessians for Large-Scale MCSCF Wavefunctions", submitted J.

Chern. Phys.