34.2.2 Automatic geometry optimization (OPTG)

The OPTG command is used to perform automatic geometry optimizations for all kinds of wavefunctions. It must directly follow the input for the wavefunction used in the geometry optimization. OPTG will call FORCE, OPT, INT, and, as needed, HF, RHF, MCSCF, CI, CCSD etc. For each of these programs, the input file is automatically repositioned to the last corresponding input before the OPTG card; so any input for RHF, MCSCF, CI, CCSD etc. can be used and will be correctly processed. It is essential, however, that the most recently optimized orbitals are used in the wavefunction for which the geometry is optimized. Any input needed for OPTG must directly follow the OPTG card. The gradients are computed analytically if the geometry is optimized for SCF or MCSCF wavefunctions; otherwise the gradients are computed by finite differences (see OPTG, NUMERICAL). Davidson corrected energies or excited state energies can be optimized using the VARIABLE and STATE subdirective.

The standard MOLPRO convergency criterion requires the maximum component of the gradient to be less then $3 \cdot 10^{-4}$ [a.u.] and the maximum energy change to be less than $1 \cdot 10^{-6}$ [H] or the maximum component of the gradient to be less then $3 \cdot 10^{-4}$ [a.u.] and the maximum component of the step to be less then $3 \cdot 10^{-4}$ [a.u.].

You may also use the convergency criterion of the Gaussian program package. It is somewhat weaker than the MOLPRO criterion and requires the the maximum component of the gradient to be less then $4.5 \cdot 10^{-4}$ [a.u.] and the root mean square (RMS) of the gradient to be less then $3 \cdot 10^{-4}$ [a.u.] as well as the maximum component of the optimization step to be less then $0.0018$ [a.u.] and the RMS of the optimization step to be less then $0.0012$ [a.u.].

The automatic geometry optimization program is invoked with

OPTG,key1=value, key2=value,......

where key can be

MAXIT

to set the maximum number of optimization cycles. The default is 50.

GRAD

sets the required accuracy of the optimized gradient. The default is $3 \cdot 10^{-4}$.

STEP

to set the convergence threshold for the geometry optimization step; if $value \ge 1$, the threshold is set to $10^{-value}$. The default is $3 \cdot 10^{-4}$.

ENERGY

sets the required accuracy of the optimized energy. The default is $1 \cdot 10^{-6}$.

GAUSSIAN

Use Gaussian convergency criteria.

SRMS

sets (for Gaussian convergency criterion) the required accuracy of the RMS of the optimization step. The default is $0.0012$.

GRMS

sets (for Gaussian convergency criterion) the required accuracy of the RMS of the gradient. The default is $3 \cdot 10^{-4}$.

BAKER

Use Baker's convergency criteria (see J. Baker, J. Comp. Chem. 14,1085 (1993)).

The defaults for the convergence parameters can also be changed by using a global GTHRESH directive, i.e.

GTHRESH, OPTSTEP=step, OPTGRAD=grad, ENERGY=energy;

By default, the geometry is optimized with respect to all variables on which the geometry depends, unless GEOMTYP=XYZ is used or the COORD,3N option is set. In the latter two cases all 3N coordinates are optimized. In other cases, subsets of these can be defined using ACTIVE or INACTIVE subdirectives.

The MOLPRO geometry optimization now utilizes a force field approximation to the hessian (``Model Hessian'', see R. Lindh, A. Bernhardsson, G. Karlström and P. Malmqvist Chem. Phys. Lett. 241, 423 (1995).) which speeds up convergence significantly. The Model Hessian is parameterized for the elements up to the third row and is used by default unless the molecule contains atoms from higher rows. Use of the Model Hessian can be switched off by the NOMODEL option. If the geometry is defined as a Z-matrix, the optimization will normally be performed in all active Z-matrix coordinates. Use of the COORD,3N option allows to optimize all 3N cartesian coordinates, even if the input is given as a Z-matrix with less parameters. Similarly, if a cartesian XYZ-input is given (GEOMTYP=XYZ), all 3N coordinates are optimized, irrespective of whether they are specified using variables or numbers. In the latter cases optimization will take place in local normal coordinates. If all 3N coordinates are optimized (GEOMTYP=XYZ or COORD,3N) it is recommended to use the BMAT option in conjunction with the Model Hessian. For other options see also COORD and UPDATE sections.

The BMAT program constructs it's own internal coordinates (natural internal coordinates, see G. Fogarasi, X. Zhou, P. W. Taylor and P. Pulay J. Am. Chem. Soc. 114, 8191 (1992). For a detailed description of the coordinates see P. Pulay, G. Fogarasi, F. Pang, J. E. Boggs J. Am. Chem. Soc. 101, 2550 (1979).) from the cartesian input. These coordinates resemble in part the valence coordinates used by vibrational spectroscopists, and have the advantage of decreasing coupling between different modes. This often increases the speed of convergence (see also COORD, section 30.2.16). The use of this option is highly recommended, especially in minimization of large organic molecules with rings. Nevertheless you should keep in mind that these coordinates are constructed automatically, and there exist exotic bond structures which might not be treated properly (e.g. weakly bonded species as in transition state optimizations). In such a case, if the BMAT optimization converges slowly or leads to symmetry-breaking errors, you should try another optimization method and/or cartesian or Z-Matrix coordinates.

The geometry, energies, and gradients are stored on a geometry record (default 700). Normally, the information of each optimization step is lost after the optimization is finished. However, using the SAVE directive it is possible to define separate geometry records for each optimization, which can be used in later optimizations using the START or GETHES commands. You may use these options to improve convergency, by using the information from this geometry records in constructing the hessian. Note that use of the Model Hessian makes the START and GETHES options obsolete. Normally convergence is reached faster by using the Model Hessian (which is recalculated in every optimization step) than by using information from old geometry records.

MOLPRO allows minimization (i.e. search for equilibrium geometries), transition state optimization (i.e. search for saddle points on energy surfaces) and reaction path following. The standard algorithms are based on the rational function approach and the geometry DIIS approach (see F. Eckert, P. Pulay and H.-J. Werner, J. Comp. Chem 18, 1473 (1997)). Also available is the quadratic steepest descent following method of Sun and Ruedenberg (see J. Sun and K. Ruedenberg, J. Chem. Phys. 99, 5257 (1993)). This method is often advantageous in Transition State searches.

In transition state searches and reaction path following Z-Matrix coordinates should be used. Although it is also possible to use cartesian or BMAT coordinates, the computational effort is usually much larger, since the hessian matrix has to be calculated numerically in all $3*N$ possible degrees of freedom.

For a detailed discussion of the various minimization algorithms see (see F. Eckert, P. Pulay and H.-J. Werner, J. Comp. Chem 18, 1473 (1997)).