Opt Keyword
Last Update: 6/25/2001

DESCRIPTION

This keyword requests that a geometry optimization be performed. The geometry will be adjusted until a stationary point on the potential surface is found. Gradients will be used if available. For the Hartree-Fock, CIS, MP2, MP3, MP4(SDQ), CID, CISD, CCD, QCISD, CASSCF, and all DFT and semi-empirical methods, the default algorithm for both minimizations (optimizations to a local minimum) and optimizations to transition states and higher-order saddle points is the Berny algorithm using redundant internal coordinates [101] (specified by the Redundant option). The default algorithm in Gaussian 92 was the Berny algorithm using internal coordinates (Opt=Z-matrix) [88,100,253]. The default algorithm for all methods lacking analytic gradients is the eigenvalue-following algorithm [96] (Opt=FP).

The remainder of this quite lengthy section discusses various aspects of geometry optimizations, and it includes these subsections:

  • Options to the Opt keyword.

  • Overview of geometry optimizations in Gaussian 98.

  • Ways of generating initial force constants.

  • Optimizing to transition states and higher-order saddle points.

  • Summary of the Berny optimization algorithm.

  • Notes on optimizing in redundant internal coordinates, including examples of Opt input and output and using the ModRedundant option.

  • Examples for Opt=Z-Matrix.

Users should consult those subsection(s) that apply to their interests and needs. Basic information as well as techniques and pitfalls related to geometry optimizations are discussed in detail in chapter 3 of Exploring Chemistry with Electronic Structure Methods[304]. See also Appendix B if you are interested in details about setting up Z-matrices for various types of molecules.


GENERAL PROCEDURAL OPTIONS

MaxCycle=N
Sets the maximum number of optimization steps to N. The default is the maximum of 20 and twice the number of redundant internal coordinates in use (for the default procedure) or twice the number of variables to be optimized (for other procedures).

StepSize=N
Sets the maximum size for an optimization step to 0.01N Bohr or radians. The default value for N is 30.

TS
Requests optimization to a transition state rather than a local minimum.

Saddle=N
Requests optimization to a saddle point of order N.

QST2
Search for a transition structure using the STQN method. This option requires the reactant and product structures as input, specified in two consecutive groups of title and molecule specification sections. Note that the atoms must be specified in the same order in the two structures. TS should not be specified with QST2.

QST3
Search for a transition structure using the STQN method. This option requires the reactant, product, and initial TS structures as input, specified in three consecutive groups of title and molecule specification sections. Note that the atoms must be specified in the same order within the three structures. TS should not be specified with QST3.

Path=M
In combination with either the QST2 or the QST3 option, requests the simultaneous optimization of a transition state and an M-point reaction path in redundant internal coordinates [307]. No coordinate may be frozen during this type of calculation.

If QST2 is specified, the title and molecule specification sections for both reactant and product structures are required as input as usual. The remaining M-2 points on the path are then generated by linear interpolation between the reactant and product input structures. The highest energy structure becomes the initial guess for the transition structure. At each step in the path relaxation, the highest point at each step is optimized toward the transition structure.

If QST3 is specified, a third set of title and molecule specification sections must be included in the input as a guess for the transition state as usual. The remaining M-3 points on the path are generated by two successive linear interpolations, first between the reactant and transition structure and then between the transition structure and product. By default, the central point is optimized to the transition structure, regardless of the ordering of the energies. In this case, M must be an odd number so that the points on the path may be distributed evenly between the two sides of the transition structure.

In the output for a simultaneous optimization calculation, the predicted geometry for the optimized transition structure is followed by a list of all M converged reaction path structures.

The treatment of the input reactant and product structures is controlled by other options: OptReactant, OptProduct, and BiMolecular.

Note that the SCF wavefunction for structures in the reactant valley may be quite different from that of structures in the product valley. Guess=Always can be used to prevent the wavefunction of a reactant-like structure from being used as a guess for the wavefunction of a product-like structure.

OptReactant
Specifies that the input structure for the reactant in a simultaneous optimization calculation should be optimized to a local minimum. The default is NoOptReactant, which retains the input structure as a point that is already on the reaction path (which generally means that it should have been previously optimized to a minimum). OptReactant may not be combined with BiMolecular.

OptProduct
Specifies that the input structure for the product in a simultaneous optimization calculation should be optimized to a local minimum. The default is NoOptProduct, which retains the input structure as a point that is already on the reaction path (which generally means that it should have been previously optimized to a minimum). Optproduct may not be combined with BiMolecular.

BiMolecular
Specifies that the reactants or products are bimolecular and that the input structures will be used as an anchor point. This anchor point will not appear as one of the M points on the path. Instead, it will be used instead to control how far the product side spreads out from the transition state. By default, this option is off.

Conical
Searches for a conical intersection or avoided crossing using the state-averaged CASSCF method. See the discussion of the CASSCF keyword for details and examples. Avoided is a synonym for Conical. Note that CASSCF=SlaterDet is needed in order to locate a conical intersection between a singlet state and a triplet state.

Restart
Restarts a geometry optimization from the checkpoint file. In this case, the entire route section will consist of the Opt keyword and the same options to it as specified for the original job (along with Restart). No other input is needed (see the examples).

NoFreeze
Activates (unfreezes) all variables (normally used with Geom=Check).

ModRedundant
Add, delete or modify redundant internal coordinate definitions (including scan and constraint information). This option requires a separate input section following the geometry specification. When used in conjunction with QST2 or QST3, a ModRedundant input section must follow each geometry specification. AddRedundant is synonymous with ModRedundant.

Lines in a ModRedundant input section use the following syntax: 

[Type] N1 [N2 [N3 [N4]]] [[+=]value] [A | F] [[min] max]]
[Type] N1 [N2 [N3 [N4]]] [[+=]value] B [[min] max]]
[Type] N1 [N2 [N3 [N4]]] K | R [[min] max]]
[Type] N1 [N2 [N3 [N4]]] [[+=]value] D [[min] max]]
[Type] N1 [N2 [N3 [N4]]] [[+=]value] H diag-elem [[min] max]]
[Type] N1 [N2 [N3 [N4]]] [[+=]value] S nsteps stepsize [[min] max]]

N1, N2, N3 and N4 are atom numbers or wildcards (discussed below). Atom numbering begins at 1, and any dummy atoms are not counted. Value specifies a new value for the specified coordinate, and +=value increments the coordinate by value.

The atom numbers and coordinate value are followed by a one-character code letter indicating the coordinate modification to be performed; the action code is sometimes followed by additional required parameters as indicated above. If no action code is included, the default action is to add the specified coordinate. These are the available action codes:

A
Activate the coordinate for optimization if it has been frozen.

F
Freeze the coordinate in the optimization.

B
Add the coordinate and build all related coordinates.

K
Remove the coordinate and kill all related coordinates containing this coordinate.

R
Remove the coordinate from the definition list (but not the related coordinates).

D
Calculate numerical second derivatives for the row and column of the initial Hessian for this coordinate.

H
Change the diagonal element for this coordinate in the initial Hessian to diag-elem.

S
Perform a relaxed potential energy surface scan. Set the initial value of this coordinate to value (or its current value), and increment the coordinate by stepsize a total of nsteps times, performing an optimization from each resulting starting geometry.

An asterisk (*) in the place of an atom number indicates a wildcard. Min and max then define a range (or maximum value if min is not given) for coordinate specifications containing wildcards. The action specified by the action code is taken only if the value of the coordinate is in the range.

Here are some examples of wildcard use:

*
All atoms specified by Cartesian coordinates

* *
All defined bonds

3 *
All defined bonds with atom 3

* * *
All defined valence angles

* 4 *
All defined valence angles around atom 4

* * * *
All defined dihedral angles

* 3 4 *
All defined dihedral angles around the bond connecting atoms 3 and 4

When the action codes K and B are used with one or two atoms, the meaning of a wildcard is extended to include all applicable atoms, not just those involving defined coordinates. By default, the coordinate type is determined from the number of atoms specified: Cartesian coordinates for 1 atom, bond stretch for 2 atoms, valence angle for 3 atoms and dihedral angle for 4 atoms. Optionally, Type can be used to designate these and additional coordinate types:

X
Cartesian coordinates. In this case, value, min and max are each triples of numbers, specifying the X,Y,Z coordinates.

B
Bond length

A
Valence angle

D
Dihedral angle

L
Linear bend specified by three atoms (or if N4 is -1) or by four atoms, where the fourth atom is used to determine the 2 orthogonal directions of the linear bend. In this case, value, min and max are each pairs of numbers, specifying the two orthogonal bending components.

O
Out-of-plane bending coordinate for a center (N1) and three connected atoms.

See the examples later in this section for illustrations of the use of this keyword.

HBond
Specifies how hydrogen bonds are generated as defined redundant internal coordinates. HBond says to generate hydrogen bonds but not to build extra bond angles and dihedral angles for them, and it is the default. NoHBond turns off the hydrogen bond generation, and AllHBond defines all hydrogen bonds and all related bond angles and dihedral angles involving those bonds. These options are valid only for optimizations in redundant internal coordinates.

Hindered
Requests the identification of internal rotation modes for the optimized structure during the vibrational frequency analysis at the end of an Opt=CalcAll optimization [342]. If any normal modes are identified as internal rotation, hindered or free, the thermodynamic functions are corrected. Hindered is valid only for optimizations performed in redundant internal coordinates and must be combined with the CalcAll option.


COORDINATE SYSTEM SELECTION OPTIONS

Redundant
Perform the optimization using the Berny algorithm in redundant internal coordinates. This is the default for methods for which analytic gradients are available.

Z-matrix
Perform the optimization in internal coordinates. In this case, the keyword FOpt rather than Opt requests that the program verify that a full optimization is being done (i.e., that the variables including inactive variables are linearly independent and span the degrees of freedom allowed by the molecular symmetry). The POpt form requests a partial optimization in internal coordinates. It also suppresses the frequency analysis at the end of optimizations which include second derivatives at every point (via the CalcAll option).

Cartesian
Requests that the optimization be performed in Cartesian coordinates, using the Berny algorithm. Note that the initial structure may input using any coordinate system. No partial optimization or freezing of variables can be done with purely Cartesian optimizations; the mixed optimization format with all atoms specified via Cartesian lines in the Z-matrix can be used along with Opt=Z-matrix if these features are needed (see Appendix B for details and examples).

When a Z-matrix without any variables is used for the molecule specification,and Opt=Z-matrix is specified, then the optimization will actually be performed in Cartesian coordinates.

Note that a variety of other coordinate systems, such as distance matrix coordinates, can be constructed using the ModRedundant option.


OPTIONS RELATED TO INITIAL FORCE CONSTANTS

EstmFC
Estimate the force constants using the old diagonal guesses. Only available for the Berny algorithm.

NewEstmFC
Estimate the force constants using a valence force field. This is the default.

ReadFC
Extract force constants from a checkpoint file. These will typically be the final approximate force constants from an optimization at a lower level, or the force constants computed correctly by a lower-level frequency calculation (the latter are greatly preferable to the former).

StarOnly
Specifies that the specified force constants are to be estimated numerically but that no optimization is to be done. This has nothing to do with computation of vibrational frequencies. In order to pass force constants estimated in this way to the Murtaugh-Sargent program, it is necessary to do one run with Opt=StarOnly to produce the force constants, and then run the actual optimization with Opt(MS,ReadFC).

FCCards
Reads the Cartesian forces and force constants from the input stream after the molecule specifications. This can be used to read force constants recovered from the Quantum Chemistry Archive using its internal FCList command. The format for this input is:

Energy (format D24.16)
Cartesian forces (lines of format 6F12.8)
Force constants (lines of format 6F12.8)

The force constants are in lower triangular form--((F(J,I),J=1,I),I=1,NAt3), where NAt3 is the number of Cartesian coordinates.

RCFC
Specifies that the computed force constants in Cartesian coordinates from a frequency calculation are to be read from the checkpoint file. This is used when the definitions of variables are changed, making previous internal coordinate force constants useless. ReadCartesianFC is a synonym for RCFC.

CalcHFFC
Specifies that the analytic HF force constants are to be computed at the first point. Note that this option is equivalent to CalcFC for DFT methods.

CalcFC
Specifies that the force constants be computed at the first point using the current method (available for the HF, MP2, CASSCF, DFT, and semi-empirical methods only).

CalcAll
Specifies that the force constants are to be computed at every point using the current method (available for the HF, MP2, CASSCF, DFT, and semi-empirical methods only). Note that vibrational frequency analysis is automatically done at the converged structure and the results of the calculation are archived as a frequency job.

VCD
Calculates VCD intensities at each point of an Opt=CalcAll optimization.

NoRaman
Specifies that Raman intensities are not to be calculated at each point of an Opt=CalcAll job (since it includes a frequency analysis using the results of the final point of the optimization). The Raman intensities add 10-20% to the cost of each intermediate second derivative point.


CONVERGENCE-RELATED OPTIONS

These options are available for the Berny algorithm only.

Tight
This option tightens the cutoffs on forces and step size that are used to determine convergence. An optimization with Opt=Tight will take several more steps than with the default cutoffs. For molecular systems with very small force constants (low frequency vibrational modes), this may be necessary to ensure adequate convergence and reliability of frequencies computed in a subsequent job step. This option can only be used with Berny optimizations.

VeryTight
Extremely tight optimization convergence criteria. VTight is a synonym for VeryTight.

EigenTest
Requests and NoEigenTest suppresses testing the curvature in Berny optimizations. The test is on by default only for transition states in internal (Z-matrix) or Cartesian coordinates, for which it is recommended. Occasionally, transition state optimizations converge even if the test is not passed, but NoEigenTest is only recommended for those with large computing budgets.

Expert
Relaxes various limits on maximum and minimum force constants and step sizes enforced by the Berny program. This option can lead to faster convergence but is quite dangerous. It is used by experts in cases where the forces and force constants are very different from typical molecules and Z-matrices, and sometimes in conjunction with Opt=CalcFC or Opt=CalcAll. NoExpert enforces the default limits and is the default.

Loose
Sets optimization conversion criteria to a maximum step size of 0.01 au and an RMS force of 0.0017 au. These values are consistent with Int(Grid=SG1) and may be appropriate for initial optimizations of large molecules using DFT methods which are intended to be followed by a full convergence optimization using the default (Fine) grid. It is not recommended for use by itself.


ALGORITHM-RELATED OPTIONS

Big
Requests the optimization to be done using the fast equation solving methods [308] for the coordinate transformations and the Newton-Raphson or RFO step. This option is default for semiempirical calculations. The use of Opt=Big is recommended for optimization of very large systems when the energy and gradient is calculated by O(N) scaled methods. This option can be turned off using Opt=Small. Large is a synonym for Big.

This method avoids the matrix diagonalizations. Consequently, the eigenvector following methods (Opt=TS) cannot be used in conjunction with it. QST2 and QST3 calculations are guided using an associated surface approximation, but this may not be as effective as the normal method involving eigenvector following.

Micro
The use of microiterations in geometry optimizations is the default for MM optimizations and ONIOM optimizations with an MM component. Use the Opt=NoMicro option to turn off microiterations.

CheckCoordinates
Rebuild the connectivity matrix before each optimization step. If there is any change in it, rebuild the redundant internal coordinate system. This option is off by default.

Linear
Requests and NoLinear suppresses the linear search in Berny optimizations. The default is to use the linear search whenever possible.

TrustUpdate
Requests and NoTrustUpdate suppresses dynamic update of the trust radius in Berny optimizations. The default is to update for minima.

RFO
Requests the Rational Function Optimization [261] step during Berny optimizations. It is the default.

GDIIS
Specifies the use of the modified GDIIS algorithm 309,310,311]. Recommended for use with large systems, tight optimizations and molecules with flat potential energy surfaces. It is the default for semiempirical calculations. This option is turned off by the RFO and Newton options.

Newton
Requests the Newton-Raphson step rather than the RFO step during Berny optimizations.

NRScale
Requests that if the step size in the Newton-Raphson step in Berny optimizations exceeds the maximum, then it is be scaled back. NoNRScale causes a minimization on the surface of the sphere of maximum step size [264]. Scaling is the default for transition state optimizations and minimizing on the sphere is the default for minimizations.

EF
Requests an eigenvalue-following algorithm [261,265,263 ]. Available for both minima and transition states, with second, first, or no analytic derivatives as indicated by CalcAll, CalcFC, the defaults, or EnOnly. EigFollow, EigenFollow, and EigenvalueFollow are all synonyms for EF. Note that when analytic gradients are available and the lowest eigenvector is being followed, then the default Berny algorithm has all of the features of the eigenvalue-following algorithm.

TVector=N
Follow the eigenvector with Nth eigenvalue of the actual Hessian. The use of Opt=CalcFC in conjunction with this option is highly recommended. TVector cannot be used in conjunction of Opt=Big.

Steep
Requests steepest descent instead of Newton-Raphson steps during Berny optimizations. This is only compatible with Berny local minimum optimizations. It may be useful when starting far from the minimum, but is unlikely to reach full convergence.


OPTIONS FOR INTERPRETING NUMERICAL ERROR

HFError
Assume that numerical errors in the energy and forces are those appropriate for HF and PSCF calculations (1.0D-07 and 1.0D-07, respectively). This is the default for optimizations using those methods.

FineGridError
Assume that numerical errors in the energy and forces are those appropriate for DFT calculations using the default grid (1.0D-07 and 1.0D-06, respectively). This is the default for optimizations using a DFT method and using the default grid (or specifying Int=FineGrid). SEError is a synonym for this option, as these values are also appropriate for semi-empirical calculations (for which it is also the default).

SG1Error
Assume that numerical errors in the energy and forces are those appropriate for DFT calculations using the SG-1 grid (1.0D-07 and 1.0D-05, respectively). This is the default for optimizations using a DFT method and Int(Grid=SG1Grid).

ReadError
Read in the accuracy to assume for the energy and forces, in format 2F10.6 (there is no terminating blank line for this input section since it is always a single line).


OVERVIEW OF GEOMETRY OPTIMIZATIONS IN GAUSSIAN 98

By default, Gaussian 98 performs the optimization in redundant internal coordinates. This is a change from previous versions of the program. There has been substantial controversy in recent years concerning the optimal coordinate system for optimizations. For example, Cartesian coordinates were shown to be preferable to internal coordinates (Z-matrices) for some cyclic molecules [254]. Similarly, mixed internal and Cartesian coordinates were shown to have some advantages for some cases [255] (among them, ease of use in specifying certain types of molecules).

Pulay has demonstrated [256; 257;258], however, that redundant internal coordinates are the best choice for optimizing polycyclic molecules, and Baker reached a similar conclusion when he compared redundant internal coordinates to Cartesian coordinates [259]. By default, Gaussian 98 performs optimizations via the Berny algorithm in redundant internal coordinates; these new procedures are also the work of H. B. Schlegel and coworkers [101].

In addition to employing a new coordinate system, this optimization procedure operates somewhat differently from those traditionally employed in the program:

  • The choice of coordinate system for the starting molecular structure is, quite literally, irrelevant, and it has no effect on the way the optimization proceeds. All of the efficiency factors in the various coordinate systems are of no consequence, since all structures are converted internally to redundant internal coordinates.

  • All optimizations in redundant internal coordinates are full optimizations unless variables are explicitly frozen using the ModRedundant option. Including a separate constants variable section in the molecule specification does not result in any frozen variables. Similarly, the requirement that all variables in the Z-matrix be linearly independent does not apply to these optimizations. Optimizations in redundant internal coordinates do make use of geometry constraint information and numerical differentiation specifications. See the examples subsection for details.

Optimizations in internal coordinates, which was the default procedure in Gaussian 92, is still available, via the Opt=Z-Matrix option.


WAYS OF GENERATING INITIAL FORCE CONSTANTS

Unless you specify otherwise, a Berny geometry optimization starts with an initial guess for the second derivative matrix--also known as the Hessian--which is determined using connectivity determined from atomic radii and a simple valence force field 260;101]. The approximate matrix is improved at each point using the computed first derivatives.

This scheme usually works fine, but for some cases, such as Z-matrices with unusual arrangements of dummy atoms, the initial guess may be so poor that the optimization fails to start off properly or spends many early steps improving the Hessian without nearing the optimized structure. In addition, for optimizations to transition states (see also below), some knowledge of the curvature around the saddle point is essential, and the default approximate Hessian must always be improved.

In these cases, there are several methods for providing improved force constants:

  • Use force constants from a lower-level calculation:
    The force constants can be read from the checkpoint file (Opt=ReadFC). These will typically be the final approximate force constants from an optimization at a lower level or (much better) the force constants computed correctly at a lower level during a frequency calculation.

  • Extract Cartesian force constants from a checkpoint file:
    The Cartesian (as opposed to internal) force constants can be read from the checkpoint file. Normally it is preferable to pick up the force constants already converted to internal coordinates as described above. However, a frequency calculation occasionally reveals that a molecule needs to distort to lower symmetry. Usually this means that a new Z-matrix with fewer symmetry constraints must be specified to optimize to the lower energy structure. In this case the computed force constants in terms of the old Z-matrix variables cannot be used, and instead the command Opt=RCFC is used to read the Cartesian force constants and transform them to the current Z-matrix variables.

    Note that Cartesian force constants are only available on the checkpoint file after a frequency calculation. You cannot use this option after an optimization dies because of a wrong number of negative eigenvalues in the approximate second derivative matrix. In that case, you may want to start from the most recent geometry and compute some derivatives numerically.

  • Calculate initial force constants at the HF level:
    You can also request that the analytic Hartree-Fock second derivatives be calculated at the first point of the optimization. This can be used with HF, DFT or post-SCF gradient optimizations. This is done by specifying Opt=CalcHFFC. Note that this option is equivalent to CalcFC for DFT methods.

  • Calculate initial force constants at the current level of theory:
    You can request that the second derivatives of the method being used in the optimization be computed at the first point by specifying Opt=CalcFC. This is only possible for HF, DFT, MP2, and semi-empirical methods.

  • Calculate new force constants at every point:
    Normally after the initial force constants have been decided upon, they are updated at each point using the gradient information available from the points done in the optimization. For a Hartree-Fock, MP2, or semi-empirical optimization, you can specify Opt=CalcAll, which requests that second derivatives be computed at every point in the optimization. Needless to say, this is very expensive.

  • Input new guesses:
    The default approximate matrix can be used, but with new guesses read in for some or all of the diagonal elements of the Hessian. This is specified in the ModRedundant input or on the variable definition lines in the Z-matrix. For example:

    Redundant Internals     Z-matrix
     1 2 3 104.5          A 104.5
     1 2 1.0 H 0.55       R 1.0 H 0.55

    The first line specifies that the angle formed by atoms 1, 2 and 3 (the variable A in the Z-matrix) is to start at the value 104.5, and the second line sets the initial value of the bond between atoms 1 and 2 (the variable R in the Z-matrix) to 0.55 Angstroms. The letter H on the second line indicates that a diagonal force constant is being specified for this coordinate and that its value is 0.55 hartree/au2. Note that the units here are Hartrees and Bohrs or radians.

    This option is valid only with the Berny algorithm.

  • Compute some or all of the Hessian numerically:
    You can ask the optimization program to compute part of the second derivative matrix numerically. In this case each specified variable will be stepped in only one direction, not both up and down as would be required for an accurate determination of force constants. The resulting second-derivatives are not as good as those determined by a frequency calculation but are fine for starting an optimization. Of course, this requires that the program do an extra gradient calculation for each specified variable. This procedure is requested by a flag (D) on the variable definition lines: 

            Redundant Internals            Z-matrix
        1 2 1.0 D                      R1 1.0 D
        2 3 1.5                        R2 1.5
        1 2 3 104.5 D                  A1 104.5 D
        2 3 4 110.0                    A2 110.0

    This input tells the program to do three points before taking the first optimization step: the usual first point, a geometry with the bond between atoms 1 and 2 (R1) incremented slightly, and a geometry with the angle between atoms 1, 2 and 3 (A1) incremented slightly. The program will use the default diagonal force constants for the other two coordinates and will estimate all force constants (on and off diagonal) for bond(1,2)/R1 and angle(1,2,3)/A1 from the three points. This option is only available with the Berny and EF algorithms.


OPTIMIZING TO A TRANSITION STATE OR HIGHER-ORDER SADDLE POINT

Transition State Optimizations Using Synchronous Transit-Guided Quasi-Newton (STQN) Methods. Gaussian 98 includes a new method for locating transition structures. This STQN method, implemented by H. B. Schlegel and coworkers [102, 101], uses a quadratic synchronous transit approach to get closer to the quadratic region of the transition state and then uses a quasi-Newton or eigenvector-following algorithm to complete the optimization. Like the default algorithm for minimizations, it performs optimizations by default in redundant internal coordinates. This method will converge efficiently when provided with an empirical estimate of the Hessian and suitable starting structures.

This method is requested with the QST2 and QST3 options. QST2 requires two molecule specifications, for the reactants and products, as its input, while QST3 requires three molecule specifications: the reactants, the products, and an initial structure for the transition state, in that order. The order of the atoms must be identical within all molecule specifications. See the examples for sample input for and output from this method.

Despite the superficial similarity, this method is very different from the Linear Synchronous Transit method for locating transition structures requested with the now-deprecated LST keyword. Opt=QST2 generates a guess for the transition structure that is midway between the reactants and products in terms of redundant internal coordinates, and it then goes on to optimize that starting structure to a first-order saddle point automatically. The Linear Synchronous Transit method merely locates a maximum along a path connecting two structures which may be used as a starting structure for a subsequent manually-initiated transition state optimization; LST does not locate a proper stationary point. In contrast, QST2 and QST3 do locate proper transition states.

Traditional Transition State Optimizations Using the Berny Algorithm. The Berny optimization program can also optimize to a saddle point using internal coordinates, if it is coaxed along properly. The options to request this procedure are Opt=TS for a transition state (saddle point of order 1) or Opt(Saddle=N) for a saddle point which is a maximum in N directions.

When searching for a local minimum, the Berny algorithm uses a combination of rational function optimization (RFO) and linear search steps to achieve speed and reliability (as described below). This linear search step cannot be applied when searching for a transition state. Consequently, transition state optimizations are much more sensitive to the curvature of the surface. A transition state optimization should always be started using one of the options described above for specifying curvature information. Without a full second derivative matrix the initial step is dependent on the choice of coordinate system, so it is best to try to make the reaction coordinate (direction of negative curvature) correspond to one or two redundant internal coordinates or Z-matrix variables (see the examples below).

In the extreme case in which the optimization begins in a region known to have the correct curvature (e.g., starting with Opt=CalcFC) and steps into a region of undesirable curvature, the Opt=CalcAll option may be useful. This is quite expensive, but the full optimization procedure with correct second derivatives at every point will usually reach a stationary point of correct curvature if started in the desired region. For suggestions on locating transition structures, refer to the literature [100].

An eigenvalue-following (mode walking) optimization method [98, 99] can be requested by Opt=EF. This was sometimes superior to the Berny method in Gaussian 88, but since the RFO step 261{Simons, 1983 #118} has now been incorporated into the Berny algorithm, EF is seldom preferable unless its ability to follow a particular mode is needed, or gradients are not available (in which case Berny can't be used anyway). This algorithm has a dimensioning limit of 50 active variables. By default, the lowest mode is followed. This is correct when already in a region of correct curvature and when the softest mode is to be followed uphill. This default can be overridden in two ways:

  • The mode having largest magnitude component for a specific Z-matrix variable can be requested by placing a 4 on the variable definition line: 

    Ang1 104.5 4

  • The Nth mode in order of increasing Hessian eigenvalue can be requested by placing a 10 after the Nth variable definition line, as in this input file: 

    	# Opt=(EF,TS)

	HCN -> HNC transition state search
	This job deliberately follows the wrong (second) mode!

	0,1
	N
	C,1,CN
	H,1,CH,2,HCN

	CN 1.3
	CH 1.20 10	       Requests the second mode.
	HCN 60.0

By default, the Berny optimization program checks the curvature (number of negative eigenvalues) of its approximate second derivative matrix at each step of a transition state optimization. If the number is not correct (1 for a transition state), the job is aborted. A search for a minimum will often succeed in spite of bad real or approximate curvature, because the steepest descent and RFO parts of the algorithm will keep the optimization moving downward, although it may also indicate that the optimization has moved away from the desired minimum and is headed through a transition state and on to a different minimum. On the other hand, a transition state optimization has less chance of success if the curvature is wrong at the current point. However, the test can be suppressed with the NoEigenTest option. If NoEigenTest is used, it is best to MaxCycle to a small value (e.g. 5) and check the structure after a few iterations.


THE BERNY OPTIMIZATION ALGORITHM

The Berny geometry optimization algorithm in Gaussian 98 is based on an earlier program written by H. B. Schlegel which implemented his published algorithm [88]. The program has been considerably enhanced since this earlier version using techniques either taken from other algorithms or never published, and consequently it is appropriate to summarize the current status of the Berny algorithm here. At each step of a Berny optimization the following actions are taken:

  • The Hessian is updated unless an analytic Hessian has been computed or it is the first step, in which case an estimate of the Hessian is made. By default, the update is done using an iterated BFGS for minima and an iterated Bofill for transition states in redundant internal coordinates, and using a modification of the original Schlegel update procedure for optimizations in internal coordinates.Normally, this is derived from a valence force field [260], but upon request either a unit matrix or a diagonal Hessian can also be generated as estimates.

  • The trust radius (maximum allowed Newton-Raphson step) is updated if a minimum is sought, using the method of Fletcher 262,554,555].

  • Any components of the gradient vector corresponding to frozen variables are set to zero or projected out, thereby eliminating their direct contribution to the next optimization step.

    If a minimum is sought, perform a linear search between the latest point and the best previous point (the previous point having lowest energy). If second derivatives are available at both points and a minimum is sought, a quintic polynomial fit is attempted first; if it does not have a minimum in the acceptable range (see below) or if second derivatives are not available, a constrained quartic fit is attempted. This fits a quartic polynomial to the energy and first derivative (along the connecting line) at the two points with the constraint that the second derivative of the polynomial just reach zero at its minimum, thereby ensuring that the polynomial itself has exactly one minimum. If this fit fails or if the resulting step is unacceptable, a simple cubic is fit is done

    Any quintic or quartic step is considered acceptable if the latest point is the best so far but if the newest point is not the best, the linear search must return a point in between the most recent and the best step to be acceptable. Cubic steps are never accepted unless they are in between the two points or no larger than the previous step. Finally, if all fits fail and the most recent step is the best so far, no linear step is taken. If all fits fail and the most recent step is not the best, the linear step is taken to the midpoint of the line connecting the most recent and the best previous points.

  • If the latest point is the best so far or if a transition state is sought, a quadratic step is determined using the current (possibly approximate) second derivatives. If a linear search was done, the quadratic step is taken from the point extrapolated using the linear search and uses forces at that point estimated by interpolating between the forces at the two points used in the linear search. By default, this step uses the Rational Function Optimization (RFO) approach [98,99,261,263]. The RFO step behaves better than the Newton-Raphson method used in earlier versions of Gaussian when the curvature at the current point is not that desired. The old Newton-Raphson step is available as an option.

  • Any components of the step vector resulting from the quadratic step corresponding to frozen variables are set to zero or projected out.

  • If the quadratic step exceeds the trust radius and a minimum is sought, the step is reduced in length to the trust radius by searching for a minimum of the quadratic function on the sphere having the trust radius, as discussed by Jorgensen [264]. If a transition state is sought or if NRScale was requested, the quadratic step is simply scaled down to the trust radius.

  • Finally, convergence is tested against criteria for the maximum force component, root-mean-square force, maximum step component, and root-mean-square step. The step is the change between the most recent point and the next to be computed (the sum of the linear and quadratic steps).


CHANGE IN TRADITIONAL CONVERGENCE CRITERIA

There has been one small but significant change in the criteria for determining when a geometry has converged. When the forces are two orders of magnitude smaller than the cutoff value (i.e., 1/100th of the limiting value), then the geometry is considered converged even if the displacement is larger than the cutoff value. This test was introduced to facilitate optimizations of large molecules which may have a very flat potential energy surface around the minimum.


AVAILABILITY

Analytic gradients are available for the HF, all DFT methods, CIS, MP2, MP3, MP4(SDQ), CID, CISD, CCD, QCISD, CASSCF, and all semi-empirical methods.

The Tight, VeryTight, Expert, Eigentest and EstmFC options are available for the Berny algorithm only.


RELATED KEYWORDS

Scan, IRC, NoSymm


EXAMPLES ILLUSTRATING OPTIMIZATIONS IN REDUNDANT INTERNAL COORDINATES

The examples in the subsection will focus on optimization procedures new to Gaussian 98. However, at the end of the subsection, examples illustrating traditional, Z-matrix-based optimizations using the Berny algorithm will also be given.

Basic Optimization Input. Traditionally, geometry optimizations required a Z-matrix specifying both the starting geometry and the variables to be optimized. For example, the input file in the left column below could be used for a such an optimization on water: 

# HF/6-31G(d) Opt Test     # HF/6-31G(d) Opt Test

Water opt                  Water opt

0  1                       0  1
O1                         O  0.00  0.00  0.00
H1 O1 R                    H  0.00  0.00  1.00
H2 O1 R H1 A               H  0.97  0.00 -0.25
   Variables:
R=1.0
A=104.5

This Z-matrix specifies the starting configuration of the nuclei in the water molecule. It also specifies that the optimization should determine the values of R and A which minimize the energy. Since the OH bond distance is specified using the same variable for both hydrogen atoms, this Z-matrix also imposes (appropriate) symmetry constraints on the molecule.

The Cartesian coordinate input in the right column is equivalent to the Z-matrix in the left column. In earlier versions of Gaussian, such input would lead to an optimization performed in Cartesian coordinates; in Gaussian 92, Z-matrix input could be used for optimizations in either coordinate system.

By contrast, in Gaussian 98 these two input files are exactly equivalent. They both will result in a Berny optimization in redundant internal coordinates, giving identical final output.

Output from Optimization Jobs. The string GradGradGrad... delimits the output from the Berny optimization procedures. On the first, initialization pass, the program prints a table giving the initial values of the variables to be optimized. For optimizations in redundant internal coordinates, all coordinates in use are displayed in the table (not merely those present in the molecule specification section): 

 GradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGrad
 Berny optimization. 
 Initialization pass.
                     ----------------------------
                     !    Initial Parameters    !
                     ! (Angstroms and Degrees)  !
 --------------------                            --------------------
 ! Name  Definition              Value          Derivative Info.    !
 --------------------------------------------------------------------
 ! R1    R(2,1)                  1.             estimate D2E/DX2    !
 ! R2    R(3,1)                  1.             estimate D2E/DX2    !
 ! A1    A(2,1,3)              104.5            estimate D2E/DX2    !
 --------------------------------------------------------------------

The manner in which the initial second derivative are provided is indicated under the heading Derivative Info. In this case the second derivatives will be estimated.

Each subsequent step of the optimization is delimited by lines like these: 

 GradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGrad
 Berny optimization. THE OPT. ALGORITHM IS IDENTIFIED BY THE HEADER FORMAT & THIS LINE
 Search for a local minimum.
 Step number   4 out of a maximum of  20

Once the optimization completes, the final structure is displayed: 

 Optimization completed.
    -- Stationary point found.
                     ----------------------------
                     !   Optimized Parameters   !
                     ! (Angstroms and Degrees)  !
 --------------------                            --------------------
 ! Name  Definition              Value          Derivative Info.    !
 --------------------------------------------------------------------
 ! R1    R(2,1)                  0.9892         -DE/DX =    0.0002  !
 ! R2    R(3,1)                  0.9892         -DE/DX =    0.0002  !
 ! A1    A(2,1,3)              100.004          -DE/DX =    0.0001  !
 --------------------------------------------------------------------

The redundant internal coordinate definitions are given in the second column of the table. The numbers in parentheses refer to the atoms within the molecule specification. For example, the variable R1, defined as R(2,1), specifies the bond length between atoms 1 and 2.

When a Z-matrix was used for the initial molecule specification, this output will be followed by an expression of the optimized structure in that format, whenever possible.

The energy for the optimized structure will be found in the output from the final optimization step, which precedes this table in the output file.

More detailed information about the output from geometry optimizations is provided in chapter 3 of Exploring Chemistry with Electronic Structure Methods.

Compound Jobs. Optimizations are commonly followed by frequency calculations at the optimized structure. To facilitate this procedure, the Opt keyword may be combined with Freq in the route section of an input file, and this combination will automatically generate a two-step job.

It is also common to follow an optimization with a single point energy calculation at a higher level of theory. The following soroute section automatically performs a HF/6-31G(d,p) optimization followed by an MP4/6-31G(d,p) single point energy calculation :

# MP4/6-31G(d,p)//HF/6-31G(d,p) Test

Note that the Opt keyword is required in this case. However, it may be included if setting any of its options is desired.

Specifying Redundant Internal Coordinates. The following input file illustrates the method for specifying redundant internal coordinates within an input file: 

# HF/6-31G(d) Opt=ModRedun Test

Opt job

0,1
C1  0.000   0.000   0.000
C2  0.000   0.000   1.505
O3  1.047   0.000  -0.651
H4 -1.000  -0.006  -0.484
H5 -0.735   0.755   1.898
H6 -0.295  -1.024   1.866
O7  1.242   0.364   2.065
H8  1.938  -0.001   1.499

3  8 
2  1  3

This structure is acetaldehyde with an OH substituted for one of the hydrogens in the methyl group; the first input line for ModRedundant creates a hydrogen bond between that hydrogen atom and the oxygen atom in the carbonyl group. Note that this line adds only the bond between these two atoms. The associated angles and dihedral angles would need to be added as well if they were desired.

Displaying the Value of a Desired Coordinate. The second input line for ModRedundant specifies the C-C=O bond angle, ensuring that its value will be displayed in the summary structure table for each optimization step.

Using Wildcards in Redundant Internal Coordinates. A distance matrix coordinate system can be activated using the following input: 

* * B   Define all bonds between pairs of atoms
* * * K           Remove all other redundant internal coordinates

The following input defines partial distance matrix coordinates to connect only the closest layers of atoms:

* * B 1.1      Define all bonds between atoms within 1.1 À
* * * K      Remove all other redundant internal coordinates

The following input sets up an optimization in redundant internal coordinates in which atoms N1 through Nn are frozen. Note that the lines containing the B action code will generate Cartesian coordinates for all of the coordinates involving the specified atom since only one atom number is specified: 

N1 B          Generate Cartesian coordinates involving atom N1
...
Nn B          Generate Cartesian coordinates involving atom Nn
* F           Freeze all Cartesian coordinates

The following input defines special "spherical" internal coordinate appropriate for molecules like C60 [312] by removing all dihedral angles from the redundant internal coordinates: 

* * * *  R     Remove all dihedral angles

The following input rotates the group about the N2-N3 bond by 10 degrees: 

* N2 N3 *  +=10.0     Add 10.0 to the values to dihedrals of the N2-N3 bond

Additional examples are found in the section on relaxed PES scans below.

Performing Partial Optimizations in Gaussian 98. The following job illustrates the method for freezing variables during a redundant internal coordinate optimization: 

# HF/6-31G* Opt=ModRedundant Test

Partial optimization

1 1
C 
H 1 R1
H 1 R1 2 A1
O 1 R2 2 A2 3 120.0
H 4 R3 3 A3 2 180.0

A1=120.0
...
R3=1.1

4 5        1.3 F
5 4 3 2        F

The structure is specified as a traditional Z-matrix, with its variables defined in a separate section. The final input section gives the values for the ModRedundant option. This input fixes the O-H bond and the dihedral angle for the final hydrogen atom. Note that any value specified in this manner need not be the same as the one listed in the preceding Z-matrix (as is the case for the O-H bond length); the structure is adjusted to enforce this constraint. The constrained value is optional. For example, in this case the value of second modified redundant internal coordinate defaults to the value from the Z-matrix (180.0).

Modifying Optimized Structures (Why You DonÀt Need a Z-matrix). Use the Cartesian coordinates version of the optimized structure as your starting point. It can be generated by a route like this one: 

# Guess=Only Geom=Check

(It can also be extracted from an archive entry.) Once you have the structure in Cartesian coordinates, you can use it in a variety of ways:

  • Add and/or remove atoms from it. Additional atoms may be specified in either Cartesian or internal coordinates.

  • Modify it by substituting atoms or groups: For example, you could change a hydrogen to a methyl group by editing the structure, replacing the desired hydrogen with a carbon atom, and then adding three additional hydrogen atoms bonded to that carbon. The latter could be given in internal coordinates: 

    H6 1.2 2.3 1.1            H6 1.2 2.3 1.1
H7 1.2 0.0 -.9            C7 1.2 0.0 -.9
H8 0.0 -.9 0.0            H8 0.0 -.9 0.0
                          H9  C7 R H5 A C2 180.0
                          H10 C7 R H6 A C2 180.0
                          H11 C7 R H8 A C2 -180.0

                          R=1.0
                          A=120.0
                          7 2 1.5

    The new structure on the right also uses an additional redundant internal coordinate (specifying Opt=ModRedundant on the final job) to alter the bond distance for the new carbon atom which is replacing the hydrogen (bonded to atom 2).

    If all you want to do is change the value or actiivate/frozen status of one or more variables, then you can use Geom=ModRedundant rather than this approach.

Restarting an Optimization. A failed optimization may be restarted from its checkpoint file by simply repeating the route section of the original job, adding the Restart option to the Opt keyword. For example, this route section restarts a Berny optimization to a second-order saddle point:

# RHF/6-31G(d) Opt=(Saddle=2,Restart,MaxCyc=50) Test

Reading a Structure from the Checkpoint File. Redundant internal coordinate structures may be retrieved from the checkpoint file with Geom=Checkpoint as usual. The read-in structure may be altered by specifying Geom=ModRedundant as well; modifications have a form identical to the input for Opt=ModRedundant:

[Type] N1 [N2 [N3 [N4]]] [[+=]Value] [Action [Params]] [[Min] Max]]

Locating a Transition Structure with the STQN Method. The QST2 option initiates a search for a transition structure connecting specific reactants and products. The input for this option has this general structure:

#HF/6-31G(d,p)Opt=QST2Test #Hf/6-31G(d)(Opt=QST2,ModRedun)Test 
 First title section
Moleculespecificationforthereactants
Second title section
Molecule specificationfortheproducts		 First title section
Molecule specification for the reactants ModRedundant input for the reactants (optional) Second title section
Molecule specification for the products ModRedundant input for the products (optional) 

Note that each molecule specification is preceded by its own title section (and separating blank line). If the ModRedundant option is specified, then each molecule specification is followed by any desired modifications to the redundant internal coordinates.

Gaussian 98 will automatically generate a starting structure for the transition structure midway between the reactant and product structures, and then perform an optimization to a first-order saddle point.

The QST3 option allows you to specify a better initial structure for the transition state. It requires the two title and molecule specification sections for the reactants and products as for QST2 and also additional, third title and molecule specification sections for the initial transition state geometry (along with the usual blank line separators), as well as three corresponding modifications to the redundant internal corrdinates if the ModRedundant option is specified. The program will then locate the transition structure connecting the reactants and products closest to the specified initial geometry.

The optimized structure found by QST2 or QST3 appears in the output in a format similar to that for other types of geometry optimizations: 

                    ----------------------------
                    !   Optimized Parameters   !
                    ! (Angstroms and Degrees)  !
---------------------                          ---------------------
! Name  Definition    Value    Reactant  Product  Derivative Info. !
--------------------------------------------------------------------
! R1    R(2,1)        1.0836    1.083     1.084   -DE/DX =    0.   !
! R2    R(3,1)        1.4233    1.4047    1.4426  -DE/DX =   -0.   !
! R3    R(4,1)        1.4154    1.4347    1.3952  -DE/DX =   -0.   !
! R4    R(5,3)        1.3989    1.3989    1.3984  -DE/DX =    0.   !
! R5    R(6,3)        1.1009    1.0985    1.0995  -DE/DX =    0.   !
! ...                                                              !
--------------------------------------------------------------------

In addition to listing the optimized values, the table includes those for the reactants and products.

Performing a Relaxed Potential Energy Surface Scan. The Opt=Z-matrix and Opt=ModRedundant keywords may also be used to perform a relaxed potential energy surface (PES) scan. Like the scan facility provided by previous versions of Gaussian, a relaxed PES scan steps over a rectangular grid on the PES involving selected internal coordinates. It differs from the operation of the Scan keyword in that a constrained geometry optimization is performed at each point.

Relaxed PES scans are available only for the Berny algorithm. If any scanning variable breaks symmetry during the calculation, then you must include NoSymm in the route section of the job, or it will fail with an error.

Redundant internal coordinates specified with the Opt=ModRedundant option may be scanned using the S code letter: N1 N2 [N3 [N4]] [[+=]value] S steps step-size. For example, this input adds a bond between atoms 2 and 3, setting its initial value to 1.0 Å, and specifying three scan steps of 0.05 Å each: 

2 3 1.0 S 3 0.05

Wildcards in the ModRedundant input may also be useful in setting up relaxed PES scans. For example, the following input is appropriate for a potential energy surface scan involving the N1-N2-N3-N4 dihedral angle. Note that all other dihedrals around the bond should be removed: 

* N2 N3 * R	       Remove all dihedrals involving the N2-N3 bond
N1 N2 N3 N4 S 20 2.0	   Relaxed PES scan of 20 steps in 2° increments


EXAMPLES FOR OPTIMIZATIONS IN INTERNAL (Z-MATRIX) COORDINATES

Full vs. Partial Optimizations. When it is performed in internal (Z-matrix) coordinates, the Berny optimization algorithm makes a distinction between full and partial optimizations. Full optimizations optimize all specified variables in order to find the lowest energy structure, while partial optimizations optimize only a specified subset of the variables. Note that the Fopt keyword is used to request that the optimization variables be tested for linear independance prior to beginning the optimization.

Those variables whose values should be held fixed are specified in a separate input section, separated by the usual variables section by a blank line or a line containing a space in the first column and the string Constants:. For example, the following input file will optimize only the bond distance R, but not the angle A, which will be held fixed at 105.4 degrees throughout the optimization: 

# HF/6-31G(d) Opt Test

Partial optimization for water

0 1
O
H1 O R
H2 O R H1 A
   Variables:
R 1.0
   Constants:
A 105.4

Freezing and Unfreezing Optimization Variables. When a geometry is retrieved from a checkpoint file, the active/frozen status as well as the values of variables may be changed for a subsequent optimization by using the Geom=Modify keyword. Consider the following input file, which performs a partial optimization on methanol, optimizing only the dihedral angles: 

%Chk=methanol
# HF/6-31G(d) Opt=Z-matrix Test

Methanol--optimize just dihedral angles

0,1
C
C 1 CO
H 2 CH 1 COH
H 1 CHA 2 OCHA 3  180.0
H 1 CHB 2 OCHB 3  HOCH
H 1 CHB 2 OCHB 3 -HCCH
   Variables:
HOCH=65.0
   Constants:
CC=1.4
CH=1.0
remaining variables ...

The geometry from the checkpoint file created by the previous job may be used for another optimization, using the following input file, which optimizes the bond angles and dihedral angles within the molecule: 

%Chk=methanol
# HF/6-31G(d) Opt(ReadFC,Z-matrix) Geom=Modify Guess=Read Test

Methanol--optimize just bond and dihedral angles, 
changing the values of CC, CH, and COH

0,1

CC=1.45
CH=0.95
COH=105.0 A
OCHA A
OCHB A

The final input section contains modifications to the geometry retrieved from the checkpoint file. The A code letter indicates that a variable is to be made active (i.e., optimized); the code letter F freezes a variable's value. If no code letter is specified, the variable's status remains the same as it was in the original job. Variables of either type may also be assigned new values (as illustrated).

The following input file illustrates the method for activating all variables in a subsequent optimization, using the NoFreeze option to the Opt keyword: 

%Chk=methanol
# HF/6-31G(d) Opt(Z-matrix,NoFreeze,ReadFC) 
# Geom=Check Guess=Read Test

Methanol--optimize everything

0,1

Note that NoFreeze is not needed for a subsequent Berny optimization using redundant internal coordinates, as this procedure ignores the distinction between active and frozen variables.

Breaking Symmetry During an Optimization in Internal Coordinates.Below are two geometry specifications for water. The one on the left has been constrained to C2v symmetry; since the same variable is used for both bond lengths, their values will always be the same: 

O                 O
H 1 R1            H 1 R1
H 1 R1 2 A        H 2 R2 2 A

R=0.9             R1=0.9
A=105.4           R2=1.1
                  A=105.4

By contrast, the Z-matrix on the right is unconstrained since the two bond lengths are specified by different variables having different initial values. Note that an optimization in redundant internal coordinates which begins from a C2v structure will retain that symmetry throughout the optimization.

Relaxed PES Scans. For Opt=Z-matrix, a relaxed PES scan is requested simply by tagging the Z-matrix variables whose values are to be incremented with the S code letter and the number of steps and the increment size. For example, the following input file requests a relaxed PES scan for the given molecule: 

# HF/6-31G(d) Opt=Z-matrix Test

Relaxed PES scan

0 1
O
H 1 R1
C 1 R2 2 A2
...
   Variables:
R1 0.9 S 5 0.05
R1 1.1
A2 115.4 S 2 1.0
...

This causes the variable R1 to be incremented five times, by 0.05 Å each time, and the variable A2 to be incremented twice, by 1 degree each time, resulting in a total of 18 geometry optimizations (the initial values for each variable also constitute a point within the scan).

Transition State Optimizations Using Internal Coordinates. We'll conclude the examples of geometry optimizations by considering traditional ways of approaching transition state optimizations (using the Opt=TS). Here is a simple example input file for determining the transition state for the 1,2-hydrogen shift which converts HCN into HNC: 

# HF/6-31G(d) Opt(TS) Test

HCN --> HNC TS Opt

0 1
C
N 1 NC
X 2 QN 1 90.0
H 3 HQ 2 90.0 1 0.0

NC=1.18
HQ=0.8 D
QN=1.15

This Z-matrix uses a dummy atom so that the shift of the hydrogen is largely given by a single coordinate (HQ). The force constants for this variable are estimated numerically giving (hopefully) one negative eigenvalue for the force constant matrix. This is noted in the output in the table labeled Initial Parameters: 

 !   HQ        0.8      calc D2E/DXDY, stepsize =    0.0026       !

The first point of the optimization proceeds as usual. At the start of the second point, Link 103 notes that the second derivatives are being estimated numerically : 

 GradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGrad
 Berny optimization.
 Numerically estimating second derivatives.

The first step is taken once this estimation is complete. The output for that point includes the eigenvector corresponding to the negative eigenvector of the Hessian matrix (i.e., the approximate reaction coordinate): 

 The second derivative matrix:
                          NC        HQ        QN
           NC           1.15446
           HQ            .07430   -.05078
           QN            .00000   -.02113    .28679
     Eigenvalues ---    -.05664    .28808   1.15902
 Eigenvectors required to have negative eigenvalues:
                          NC        HQ        QN
         1              -.06112    .99625    .06129

This output can be helpful in determining whether an optimization is proceeding toward the desired transition state.

The following input file is a two-step job illustrating a useful technique for locating transition state structures. It locates the 1,2-hydrogen shift transition structure for hydrogen cyanide: 

%Chk=hcn_ts
# HF/6-31G(d) Opt=Z-matrix Test

HCN --> CNH Initial partial opt

0 1
C
N 1 CN
H 1 CH 2 A

CN 1.3
CH 1.2

A 70.0

In the first step, the bond angle has been guessed to be about 70 degrees. The corresponding variable is held fixed while the energy is minimized with respect to the two bond lengths. The transition state is located in a second job step, starting from this partially-optimized structure: 

--Link1--
%Chk=hcn_ts
%NoSave
# HF/6-31G(d) Opt(TS,Z-matrix,ReadFC) Geom=Modify Guess=Read Test

CN --> CNH Final TS opt

0 1

A D

The single line in the geometry modification input section specifies numerical differentiation with respect to the variable A be performed before proceeding with the optimization. Note that requesting numerical differentiation with the D code letter automatically activates the variable. Since no new value is given for A, it initially has the same value as was retrieved from the checkpoint file (70 degrees). Since the angle contributes significantly to the reaction coordinate, differentiating with respect to it ensures that the initial Hessian in the second job step will have the correct curvature (one negative eigenvalue).

Return to TOC