Gaussian 03 Online Manual
Last update: 21 November 2006
Opt
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, CCSD, 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 [149,15] (specified by the Redundant option). The default algorithm for all methods lacking analytic gradients is the eigenvalue-following algorithm (Opt=EF). The Berny algorithm using internal coordinates (Opt=Z-matrix) is also available [136,148,529].
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 03.
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 [308]. 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).
MaxStep=N
Sets the maximum size for an optimization
step (the initial trust radius) 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 [164]. 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, 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. This is the default. NoOptReactant
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.
BiMolecular
Specifies
that the reactants or products are bimolecular and that the input structure 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 reactant
side spreads out from the transition state. By default, this option is off.
OptProduct
Specifies that the input structure for the product in a simultaneous optimization
calculation should be optimized to a local minimum. This is the default. NoOptProduct
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.
Conical
Search 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 interpreted as the X, Y and Z coordinates (respectively).
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.
InitialHarmonic=N
Add harmonic constraints to the initial structure with force constant N/1000
Hartree/Bohr2. IHarmonic is a synonym for this option.
ChkHarmonic=N
Add harmonic constraints to the initial structure saved on the chk file with force
constant N/1000 Hartree/Bohr2. CHarmonic is a synonym
for this option.
ReadHarmonic=N
Add harmonic constraints
to a structure read in the input stream (in the input orientation), with force
constant N/1000 Hartree/Bohr2. RHarmonic is a synonym
for this 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 be 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.
OldRedundant
Use the Gaussian 94 redundant internal coordinate generator.
Note that a variety of other coordinate systems, such as distance matrix coordinates, can be constructed using the ModRedundant option.
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. CalcHFFC is used with MP2 optimizations, and it 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
Calculate VCD intensities at each point of a Hartree-Fock Opt=CalcAll optimization.
NoRaman
Specifies that Raman intensities are not to be calculated at each point of a Hartree-Fock
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. For DFT calculations, Int=UltraFine
should be specified as well.
VeryTight
Extremely tight optimization
convergence criteria. VTight is a synonym for VeryTight. For DFT
calculations, Int=UltraFine should be specified
as well.
EigenTest
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 the optimization convergence criteria to
a maximum step size of 0.01 au and an RMS force of 0.0017 au. These values are
consistent with the Int(Grid=SG1) keyword, 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
Micro
Use microiterations in ONIOM(MO:MM) optimizations. The default, with selection
of L120 or L103 for the microiterations depending on whether electronic embedding
is on or off. NoMicro forbids microiterations during ONIOM(MO:MM) optimizations.
Mic120 says to use microiterations in L120 for ONIOM(MO:MM), even for mechanical embedding. This is the default for electronic embedding. Mic103 says to perform microiterations in L103 for ONIOM(MO:MM). It is the default for mechanical embedding, and it does not work for electronic embedding.
QuadMacro
Controls whether
the coupled, quadratic macro step is used during ONIOM(MO:MM) geometry optimizations.
This is possible with mechanical embedding but not with electronic embedding.
NoQuadMacro is the default.
Linear
Linear requests and NoLinear suppresses
the linear search in Berny optimizations. The default is to use the linear search
whenever possible.
TrustUpdate
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 [530] step during Berny optimizations.
It is the default.
GDIIS
Specifies the use of the modified GDIIS
algorithm [531,532,533].
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
Use the Newton-Raphson step rather than the RFO step during Berny optimizations.
NRScale
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 [534].
Scaling is the default for transition state optimizations and minimizing on the
sphere is the default for minimizations.
EF
Requests an eigenvalue-following
algorithm [530,535,536].
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.
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.
UpdateMethod=keyword
Specifies the Hessian update method. Keyword is one of: Powell,
BFGS, PDBFGS, ND2Corr, OD2Corr, D2CorrBFGS,
Bofill, D2CMix and None.
Big
Requests
the optimization to be done using the fast equation solving methods [537]
for the coordinate transformations and the Newton-Raphson or RFO step. This option
is default for semiempirical calculations. 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.
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
By default, Gaussian 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 [538]. Similarly, mixed internal and Cartesian coordinates were shown to have some advantages for some cases [539] (among them, ease of use in specifying certain types of molecules).
Pulay has demonstrated [540,541,542], 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 [543]. By default, Gaussian performs optimizations via the Berny algorithm in redundant internal coordinates; these procedures are also the work of H. B. Schlegel and coworkers [149].
This optimization procedure operates somewhat differently from those traditionally employed in electronic structure programs (including Gaussian 94 and earlier versions):
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 constant 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 [149,544]. 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 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 includes the STQN method for locating transition structures. This method, implemented by H. B. Schlegel and coworkers [149,150], 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 [148].
An eigenvalue-following (mode walking) optimization method [146,147] can be requested by Opt=EF. This was sometimes superior to the Berny method in Gaussian 88, but since the RFO step [530] 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 the 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 is based on an earlier program written by H. B. Schlegel which implemented his published algorithm [136]. 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. Normally 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.By default, this is derived from a valence force field [544], 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 [545,546,547].
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 [146,147,530,536]. 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 [534]. 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 BEGINNING WITH GAUSSIAN 98
Gaussian 98 introduced 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.
The generation of redundant internal coordinates for weakly bound complexes was also updated with Gaussian 98. We include Hydrogen bonds automatically. In addition, in connecting different fragments which are only weakly bound (hydrogen-bonded and otherwise), all pairs of atoms with one atom in each fragment having distance within a factor of 1.3 of the closest pair have their distances added to the internal coordinates. If at least 3 such pairs are found, then no angles or dihedrals involving both fragments are added. However, if only 1 or two pairs of atoms are close, then the related angles and dihedrals are added in order to ensure a complete coordinate system. As usual, the ModRedundant option can be used to add or remove any coordinates manually.
Analytic gradients are available for the HF, all DFT methods, CIS, MP2, MP3, MP4(SDQ), CID, CISD, CCD, CCSD, QCISD, CASSCF, and all semi-empirical methods.
The Tight, VeryTight, Expert, Eigentest and EstmFC options are available for the Berny algorithm only.
The examples in the subsection will focus on normal optimization procedures in Gaussian 03. 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 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 early versions of Gaussian, such input would lead to an optimization performed in Cartesian coordinates; however, by Gaussian 92, Z-matrix input could be used for optimizations in either coordinate system.
By contrast, beginning with Gaussian 98 these two input files are exactly equivalent, and this holds for Gaussian 03 as well. 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. The opt. algorithm is identified by the header format & this line.
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. 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 out put from geometry optimizations is provided in Chap. 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 route section automatically performs an 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 not 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 (such jobs may require the NoSymm keyword). 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 [548] 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 involving N2-N3 bond
Additional examples are found in the section on relaxed PES scans below.
Performing Partial Optimizations. 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 atoms, 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 activate/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) Opt=QST2 # HF/6-31G(d) (Opt=QST2,ModRedun)
First title section First title section
Molecule specification for the reactants Molecule specification for the reactants
Second title section ModRedundant input for the reactants
Molecule specification for the products 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 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 coordinates 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 Specify a relaxed PES scan of 20 steps in 2º increments
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 form is used to request that the optimization variables be tested for linear independence 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=Z-matrix Test Partial optimization for water 0 1 O H1 O R H2 O R H1 A Variables: R 1.0 Constants: A 105.4
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
R1=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 R2 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).