J. Phys. Chem. 1991, 95, 5104-5108
5104
Constrained Optimization of ab Initio and Semiempirlcal Hartree-Fock Wave Functions Using Mrect Minimization or Simulated Annealing Martin J. Field Division of Computer Research and Technologies, Building 12A. National Institutes of Health, Bethesda, Maryland 20892 (Received: October 25, 1990)
A method recently developed for the simultaneous minimization of the energy of a quantum mechanical system with respect to its nuclear coordinates and the variables that determine its wave function is extended to include ab initio Hartree-Fock wave functions. The orthonormality constraints on the molecular orbitals are enforced by using an iterative technique, SHAKE, and it is shown how extra constraints that permit the calculation of excited-state wave functions can be introduced. The results of a small number of calculations are presented as examples to illustrate the direct minimization and classical molecular dynamics/simulated annealing optimization procedures.
1. Introduction
Quantum mechanical (QM) techniques have traditionally been employed to study relatively small regions of the potential surfaces of molecular systems, using, for example, singlapoint calculations, local geometry optimizations, or reaction path search procedures.’ It is increasingly apparent that such approaches have their limitations and that, in many cases, it is important to perform classical or quantum dynamical simulations to obtain information about processes of chemical significance.* Classical dynamics simulations using QM potentials, in which the wave function of the system is computed along the trajectory, have been prohibitively expensive to perform. However, large increases in the processing speeds of computers and important algorithmic improvements are now making these calculations feasible. Recently, Car and Parrinello introduced a combined classical molecular dynamics (MD)/simulated annealing3 technique for use with density functional QM potentials4 in which the variables describing the wave function as well as the Cartesian coordinates defining the positions of the atoms are treated as dynamic variables. The method was used to simultaneously optimize the geometry and wave function of a crystal and it was claimed that, if the temperature of the electronic variables was held close to zero, trajectories could be generated that mimicked the true atomic dynamics of the system. Although less accurate, the approach has the advantage that it is much quicker than the exact method which solves for the wave function at each point along the atomic trajectory. In a previous paper (paper I)5 the author outlined a method, similar in spirit to that of Car and Parrinello, but applicable to neglect of differential diatomic overlap (NDDO) Hartree-Fock (HF) wave functionsa6 The first part of the present work extends the theory behind the method to a b initio H F wave functions and discusses the implementation of simulated annealing and direct minimization optimization algorithms. Much of the theory necessary was derived by Head-Gordon and Pople’ for their simultaneous optimization method for ab initio wave functions but the method in this work is more efficient because the molecular orbital (MO) coefficients themselves are taken to be the electronic (1) Hehre, W.
J.; Radom, L.; Schleyer, P. v. R.; Pople, J. A. Ab-Initio
Molecular Orbital Theory; Wiley: New York, 1986. (2) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids, 1st ed.; Oxford University Press: London, 1987. Brooks, C. L.; Karplus, M.; Pettitt, B. M. Adti. Chem. Phys. 1!M, 71. (3) Kirkpatrick, S.; Gelatt, Jr., C. D.; Vecchi, M. P. Science 1983, 220, 671. (4) Car, R.; Parrinello, M. Phys. Reti. Lett. 1985,55, 2473. For a recent review, see: Remler, D.K.; Madden, P. A. Mol. Phys. 1990, 70, 921. (5) Field, M. J. Chem. Phys. Lett. 1990, 172, 83. (6) Pople, J. A.; Beveridge. D. L. Approximate Molecular Orbital Theory; McGraw-Hill: New York, 1970. (7) Head-Gordon, M.; Pople, J . A. J . Phys. Chem. 1988, 92, 3063.
variables, rather than the Jacobi rotation parameters. To enforce the orthonormality conditions on the MO’s, the Car-Parrinello method and the procedures given here use the versatile iterative SHAKE* scheme. The second part of this paper shows how, by the introduction of extra SHAKE constraints in the minimization procedure, the HF wave functions of excited states can be computed with little extra cost compared to the equivalent ground-state calculations. The theory is applicable to both a b initio and NDDO wave functions. Section 2 of the current paper outlines the background theory required for the constrained optimization of the energy functional. It parallels closely the theory section of paper I and uses the same notation. Section 3 discusses the implementation of the direct minimization and simulated annealing MD optimization algorithms whilst section 4 provides two examples of the use of the method. The theory for the calculation of excited-state wave functions follows in section 5 and representative excited-state calculations are given in section 6. Section 7 summarizes the paper and concludes with some areas for future work. Appendix A provides a description of the implementation of the SHAKE procedure for excited state constraints. 2. Background Theory The closed-shell restricted HF wave function, q ,for a molecule which has 2M electrons is written as a single antisymmetric determinant of M doubly occupied orthonormal molecular orbitals, ~9
* = I4la(l)dl8(2)42a(3)...4MB(21M))
(1)
where a and p are electron spin eigenfunctions and the MO’s are linear combinations of N basis functions, {q,,]: N
4i =
(2)
Ecpiqp p= 1
obey the orthonormality constraint
The M O coefficients conditions, (Aij]
Aij = ~ C , , ~ S ,-, ~6,C=~ 0 ,
(3)
CY
where S,, is an overlap integral between the basis functions q,, and qua The energy, E, of the system is the expectatio? value of the wave function over the electronic Hamiltonian, H: E =
(*lfil*)
=
+ FCvI + Vnuc
CCpCvilH~v ,vi
(4)
Here HCY is a one-electron matrix element, V,,, is the nuclear(8) Ryckaert, J. P. Mol. Phys. 1985, 55, 549. (9) Szabo, A.; Ostlund, N . S.Modern Quantum Chemistry: Introduction to Adtianced Electronic Structure Theory, 2nd 4.;McGraw-Hill: New York, 1989.
This article not subject to US.Copyright. Published 1991 by the American Chemical Society
The Journal of Physical Chemistry, Vol. 95, No. 13, 1991 5105
Optimization of HartreeFock Wave Functions nuclear repulsion energy, and Fp is a Fock matrix element which is itself a function of the MO coefficients. It is written as
F,, = H,,+ E:cA~G~I~(PvIAQ) - (PAIUY)I Ad
(5)
where (~ulAa)and ( ~ A ~ I Jare u ) two-electron repulsion integrals. Normally the optimum wave function for a system is obtained by optimizing the energy expression (eq 4 ) with respect to the MO coefficients and solving the Roothaan equations that result?JO If the minimum energy conformation of the molecule is to be found at the same time as the optimum wave function then it is necessary to minimize eq 4 with respect to the atomic coordinates, {Ra},as well. For NDDO wave functions the MO coefficients, {c,,,},and the geometrical parameten, {Ra},can be chosen to be independent variables for the minimization because the NDDO constraint conditions (eq 3) are independent of the overlap matrix and so, also, of the atomic coordinates. The same choice is possible for ab initio wave functions but it leads to more complex constraint equations than in the NDDO case. Instead, the constraint equations (eq 3) can be made independent of the overlap matrix by transforming the MOs to an orthogonal basis. The new coefficients, la,}, are taken to be independent of {R,] and are defined as awl
= E("J"1 Y
where Ma is the mass of the nucleus of atom a and ma is the "mass" of the MO coefficients. The equations of motion for the MO coefficients, derived from the Lagrangian, are
(6)
urn
The exact form of the matrix, X,depends the type of orthogonalizing transformation that is employed but it must satisfy the relationship
xx,$~Ju@ =
(7)
PU
The constraints, in terms of the new variables, become
All =
Equivalent expressions for canonical orthogonalization may also be derived but, unlike eqs 12 and 13, they are odd with respect to the eigenvectors of the overlap matrix and so it is necessary to take account of their phase throughout the optimization process. Symmetric orthogonalization is preferred for the present work. To perform a local minimization of the energy expression (eq 4) it is only necessary to have defined the independent variables, their constraint conditions, and their energy derivatives, but for a MD/simulated annealing calculation a Lagrangian, &, is needed to describe the dynamics of the variables as a function of the time, t . A suitable Lagrangian is4
EU,,,~,~ - 6,] = 0 Ir
(8)
Both the simulated annealing and direct minimization optimization techniques implemented in section 3 require the first derivatives of the energy expression (eq 4) with respect to the MO and coordinate variables. The MO coefficient derivatives are easily derived and they are written as
ac, aE -aa,aE1- - E-aa,, ac,
(9)
PI
where
The coordinate derivatives of the energy are more complicated but the expression is given by Head-Gordon and Pople.' It is
Equation 11 differs in its last term from the usual expressions for ab initio wave function coordinate derivatives because the MO coefficients, {c,,,l, do not obey the Roothaan equations1° except at the end of an optimization. Expressions for the derivatives of the orthogonalizing transformation, X,depend upon its form. For symmetric orthogonalization? the matrix, X, is
X,, = (rl/2),,v = ET,t;112T,, a
(12)
where the It,] and the ITm]are the eigenvalues and eigenvectors, respectively, of the overlap matrix, S. The derivatives, derived in ref 7, are (IO) Roothaan, C. C. J. Rev. Mol. Phys. 1951, 23, 69.
and, for the atomic positions
aE Madt2= --aRa The (tklIare a set of Lagrange multipliers introduced to satisfy the constraints. The derivatives of the constraints (the second term on the right-hand side of eq 15) are d2Ra
aAkl/aarl = a , k h + a p k k
(17)
3. Implementation of the Constrained Optimization Procedure Two optimization schemes based upon the theory given in section 2 have been implemented in the DYNAMO program." One is a standard first-derivative conjugate gradient algorithm.12 It makes use of the SHAKE algorithm in the way suggested by van Gunsteren and Karplus in their procedure for the constrained geometry optimization of macrom~lecules.~~ The second is a simulated annealing3 MD technique based on the algorithms developed by Car and Parinello4 and the author in paper Iss In broad outline the schemes are as follows: Initialization: Choose starting values for the geometry and MO coefficients. Transform the M O s to the orthogonal basis using the overlap matrix at the current position. In the case of simulated annealing choose a mass for the MO coefficients and assign velocities to all the variables. I. Conjugate Gradient Minimization a. Calculate the energy at the current position and the derivatives of the energy with respect to the transformed MO coefficients and the atomic coordinates (eqs 4, 9, and 11, respectively). b. Constrain the MO derivatives by using SHAKE to remove any components that violate the constraints. c. Check for convergence of the optimization procedure. If convergence has been achieved proceed to the end. d. Perform a step of conjugate gradient optimization.I2 e. Constrain the new MO coefficients to retain orthonormality using SHAKE. f. Return to step a. (1 1) Field, M. J. DYNAMO: A Program for Performing Molecular Dynumics Simuluriom with Quuntum Mechunicul Potenriuls. Unpublished work. (12) Press, W. H.; Flannery, B. P.; Teukolsky, S.A.; Vetterling, W. T. Numeric41 Recipes: rhe Arr of Scientific Computing. Cambridge University Press: Cambridge, U.K., 1986. (13) van Gunsteren, W. F.; Karplus, M. J . Comput. Chem. 1980, I , 266.
5106 The Journal of Physical Chemistry, Vol. 95, No. 13, 1991
Field 300
TABLE I: Geometries of the Water MoleculeO
) I I I
,
,
-----
1
,
___
Total Energy __Potential Energy _ _ - Atom finetic Energy MO i-hnetic Energy
optimized geometry
starting geometry 0-H bond length 0.957 0.967 HQH bond angle 104.5 107.7 a Included are the starting model-built structure and the final optimized structure for the MD/simulated annealing calculation on water. Bond lengths are in angstroms and angles are in degrees.
i
TABLE 11: Geometries of the Formmide Moleculea
starting geometry
optimized symmetrized geometry geometry
bond lengths
C-O C-H C-N N-H (cis to oxygen) N-H (trans to oxygen) bond angles H-C-O N-C-O C-N-H (cis) C-N-H (trans)
1.200 1.100 1.350 1.Ooo 1.om
1.212 1.085 1.353 0.997 0.995
1.212 1.085 1.353 0.995 0.995
120.0 120.0 120.0 120.0
121.9 125.0 118.7 122.9
121.9 125.0 120.0 120.0
0.00
0.02
0.04
0.06
0.08
0.10
Time (ps)
Figure 1. A plot of the evolution of the total, potential, and kinetic energies along the MD trajectory for the simulated annealing simulation on water. All energies are in kJ mol-' and the time is in picoseconds. The total energy and potential energy are plotted as differences from the final optimized value of the potential energy (-75.585 960 hartrees). The values of the kinetic energies are absolute.
Included are the starting model-built structure and the final optimized structure obtained in the MD/simulated annealing calculation and the 'symmetrized" optimized structure used for the excited-state calculations. Bond lengths are in angstroms and angles are in degrees.
,
1
1400 - m 1200
I I . M D Simulated Annealing a. Calculate the energy a t the current position and the derviatives of the energy with respect to the transformed M O coefficients and the atomic coordinates (eqs 4, 9, and 11, respectively). b. Compute the new positions and atomic velocities. DYNAMO uses a Verlet integration method.I4 c. Compute the new MO coefficients. d. Constrain the MO coefficients using SHAKE. e. Calculate the MO coefficient velocities. f. Perform any analysis or velocity modification. g. Return to step a for as many steps as are required. Termination: Back-transform the MOs to the nonorthogonal basis and canonicalize them if desired. For each scheme, only one transformation between the sets of MO coefficients, {a,,,)and {c,,,),is necessary per iteration (in steps Ia and Ha). Otherwise, all manipulations are performed on the set of transformed M O coefficients, la,,). The iterative procedure, SHAKE, is used to constrain the MO coefficients and their derivatives (steps Ib, Ie, and IId). The formulas for applying the constraints are given in ref 8. No problems with their implementation or the convergence of the method have been encountered. For M orbitals there are M(M 1)/2constraints and each constraint requires O(N) manipulations, giving an overall scaling dependence of o(n,NW) for the algorithm, if ni is the number of iterations required to reach convergence. It has h e n found that n,,even for the largest cases, is usually less than 20, and so the scaling properties of the procedure compare favorably with the scaling dependence for the method of Head-Gordon and Pople' which is O(N(N-M)M). No steps in the procedure above, other than the two-electron integral calculation and Fock matrix construction (which for a standard a b initio method are proportional to O(N4)),have a scaling dependence of greater than O ( p ) .
1000
I
___ ~
- - -
1;,
I Total Energy Potential Energy Atom finetic Energy MO finetic Energy
-1
\ \ , - _ - - - - _ - - _ _ _ . _ _ _ _
+
4. Wave Function and Geometry Optimizations The combined geometry/MO coefficient optimization algorithm for a b initio wave functions is illustrated by MD/simulated annealing calculations on water and formamide with the 3-21Gbasis set.I5 Both simulations were started using model built geometries (14) Verlet, L. Phys. Reu. 1967, 159, 98. (15) Binkley, J. S.; Pople, J. A,; Hehre, W.J. J . Am. Chem. Soc. 1980, 102, 939.
~~~~
~
(16) Pople, J. A.; Beveridge, D. L.; Dobosh, P. A. J . Chem. Phys. 1967,
47, 2026.
Optimization of Hartree-Fock Wave Functions fommide, mpectively. During cooling the total energy decreased uniformly for both molecules coming within 1 kJ mol-' of the exact optimized energy after 550 cooling step for water and 850 cooling steps for formamide. In both cases the MO coefficients provide the dominant contribution to the total kinetic energy because they are more numerous, but only slightly less massive, than the coordinate degrees of freedom. The optimizations converged satisfactorily to the optimum geometry and wave functions for the molecules using the particular heating and cooling strategy outlined above. Experience with the method shows that the convergence of the calculation is relatively insensitive to the procedure adopted as long as the cooling rate is not too rapid. Undoubtedly the scheme employed here can be improved upon. By contrast, the simulations are sensitive to the values of the time step and the MO coefficient mass. If the ratio of the square of the time step to the mass is too small then convergence is very slow, but if it is too large the changes in the MO coefficients are too large and the SHAKE procedure fails to converge. In comparison to the NDDO simulations (paper I) the time step used needs to be smaller and the M O mass larger for comparable performance of the method.
5. Calculation of Excited-State Wave Functions To calculate the wave function, */,of the Ith excited state of a system within the framework of the algorithms discussed above, (I - 1) additional constraints, {AolJ],are introduced which force orthogonality between the state and the (I- 1) lower ones. They take the form Ao/j=(*/l*j)-6/j=0 VJ=l,(I-I) (18) In principle, to give an exact upper bound to the excited-state energy, additional (I - 1) constraints would need to be imposed, {Ah/,! V J = l , ( I - l ) (19) Ah/,= ( \ k / ( Q \ k j ) - E P / j = O but, because the matrix elements of states over the Hamiltonian are often very small and the extra constraints add considerably to the complexity of the algorithm, they are neglected in this work. A fuller discussion of the point is to be found in the paper by Colle, Fortunelli, and Salvetti." The M O s of different states will, in general, not be orthogonal as only orthogonality between entire wave functions is required. Formulas for the overlap matrix elements of wave functions with nonorthogonal orbitals may be found in refs 17 and 18 but two examples are given here. For two closed-shell wave functions, and W,, the overlap is given by
The Journal of Physical Chemistry, Vol. 95, No. 13, 1991 5107 The scheme currently implemented in DYNAMO^^ for the calculation of a series of excited states of a system at a single geometry is as follows. III. Excited-State Calculations a. Choose the form of the restricted HF wave function for the ground state of the system (e.g., closed-shell or open-shell singlet wave functions). b. Perform the ground-state calculation. c. Choose the form of the excited-state wave function. d. Generate a starting set of orthonormal MO's for the excited state which satisfy the state orthogonality constraints (eq 18). e. Perform an excited-state calculation using either of the schemes in section 3. The extra SHAKE constraints are introduced into steps Ib, Ie, and IId. f. If convergence has been achieved, return to step c for the next excited-state calculation. The excited-state constraints are straightforwardly implemented with the SHAKE technique and they do not appear to adversely affect either its convergence properties or those of the direct minimization or simulated annealing algorithms. Each excitedstate calculation is only slightly more expensive than a ground-state calculation as the calculation of the SHAKE constraints scales as only O(MNz). More details are given in the Appendix. It should be emphasized that the current method produces a wave function for the excited state that is optimum for that state subject only to the orthogonality restrictions to the lower states and that, in contrast to the method of Colle et al.," an expensive four-index transformation of the two-electron integrals is not required. In scheme I11 the ground- and excited-state calculations are performed separately. In many cases this is advantageous as it is possible to study the excited states of a system after the ground-state reaction path or dynamics trajectory has been determined. In some cases, though, such as when the ordering of the states is uncertain, a scheme in which the wave functions of the states are optimized simultaneously would be preferable. The method proposed by Dutta and Bhattacharyya is an example of that a p p r o a ~ h . ' ~A simultaneous optimization of several states is also the most convenient framework for optimizing the geometry of one of the excited states at the same time as the wave functions. No examples of such calculations are given in this paper.
6. Excited-State Calculations As an example of the ab initio excited-state method, direct minimization calculations were performed with the 3-21G basis seti5to determine how the excitation energy between the ground state and the first singlet excited state of the formamide molecule changed as a function of the torsion angle about the carbon(*/l*J) = (D,/,)2 (20) nitrogen bond. The geometry chosen for the molecule is given where Data is an M X M determinant of orbital overlaps in Table 11. It is similar to the optimized planar geometry for the molecule but the amide group has been "symmetrizedn so that Data = det 1 ( 4 ~ 1 4 ~ ) ( 4 ~ I 4 ~ ) . . . ( 4 ~ 1 4 (21) ~)l the molecule has C, symmetry when the C-N torsion angle is at with 90° as well as a t Oo. Calculations were performed for values of the torsion angle from 0' to 90' at approximately 5 O increments. (22) (dI4;) = C
