Hans. Aa. Jensen and Poul. Department of. Chemistry, Aarhus University, DK8000 .... d' h d,l a new step 1S eva uate W1t trust ra 1US a. Simons et a1. [2] suggest the ... condorder 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, DK8000 Arhus C, Denmark Trygve U. Helgaker, Department of Chemistry, University of Oslo, Blindern, N03l5, Oslo 3, Norway ABSTRACT. We develop a surface walking procedure based on a restrictedstep algorithm combined with the NewtonRaphson method for MCSCF potential 'energy surfaces. We discuss some of the details of the algorithm including how to optimally use it to obtain a firstorder 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 abinitio quantum chemistry is to reliably determine equilibrium and transitionstate geometries [16]. Such geometries are stationary points on a BornOppenheimer 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 [13] 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 steplength 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 NewtonRaphson step can be taken with confidence. The restrictedstep algorithm [7] combined with the NewtonRaphson 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 secondorder 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, 229241. © 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, 910].
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 largescale secondorder MCSCF optimization methods. Golab et al.[lO] used the method to demonstrate the existence of many "nonphysical" stationary points on a MCSCF energy hypersurface. Simons et al. [2] adapted and modified the procedure for walking on a geometricalsurface and applied it to model surfaces and to two potential energy surfaces of HCN. The restrictedstep 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 steplength norm for the in1tial iteration, subsequent maximum radii, hk ' increase, decrease, or stay the same based on how well the secondorder Taylor series expansion at a point with a step agrees with the exact value of the function at !. If the steplength is smaller than the trust radius and the Hessian has the correct structure (desired number of negative eigenvalues) the NewtonRaphson step is taken. This occurs in the local region of a stationary point. Consequently, the procedure has secondorder convergence. Thus the restrictedstep algorithm allows the maximum extraction of information from a secondorder procedure. If a secondorder expansion for a given step predicts a total energy which agrees fairly well with the exact value subsequently found at that point, then the secondorder expansion is a good approximation. If the secondorder expansion value does not agree very well with the exact value then the step is rejected and a new, reduced step is taken. For largescale MCSCF optimization or for MCSCF geometry walks, steprejection 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.
SecondOrder Expansions and the RestrictedStep 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 secondorder: (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 secondoiderkexpansion 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 NewtonRaphson steplength norm lies within the trustregion then the NewtonRaphson 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 2r 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 NewtonRaphson step in the consequtive iterations. As the transition state is at infinity the NewtonRaphson steps as expected did gradually increase. The increase of the NewtonRaphson 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 subminimum 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)
02H2
1.986
1.580
1.575
1.557 (1.475)
0102
37.54
102.53
96.42
100.16 (94.8)
Hl0102
95.22
102.54
96.42
100.16 (94.8)
H20201
Bond Angles (0)
"0 1{) 2" H2
H1
Olson and P.C Cross
89.51
0.00
180.00
89.89 (119.8)
H10102H2
(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
01H1
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