# Geometry optimization (OPTG)

Automatic geometry optimization is invoked using the `OPTG`

command. The `OPT`

command available in previous MOLPRO versions is no longer needed and not available any more.

`OPTG`

[, *key1=value, key2=value,……*]

The `OPTG`

command can be used to perform automatic geometry optimizations for all kinds of wavefunctions. For minimum searches, it is usually sufficient to give just the `OPTG`

command without further options or directives, but many options are available which are described in the following sections.

The `OPTG`

command must be given after the energy calculation to which it refers. If the commands for the energy calculation (e.g. `HF`

, `KS`

, `MP2`

, etc.) are in a procedure, the `OPTG`

must also be in the procedure. Furthermore, no procedures without an energy calculation must directly precede `OPTG`

.

Various optimization methods can be selected as described in section selecting the optimization method (METHOD).
Molpro allows minimization (i.e. search for equilibrium geometries), transition state optimization (i.e. search for saddle points on energy surfaces), and reaction path following. Since Molpro 2023.1 the optimizer developed by J. Baker and P. Pulay (PQSOPT) is available in Molpro and used by default. The method is described in
J. Baker et al., J. Chem. Phys. 105, 192 (1996);
J. Baker and P. Pulay, J. Chem. Phys. 105, 11100 (1996);
J. Baker, J. Comput. Chem. 18, 1079 (1997);
J. Baker and P. Pulay, J. Comput. Chem. 21, 69 (2000).
It uses delocalized internal coordinates and converges often better than previous optimisers in Molpro. Furthermore, it allows for constraint optimisations. Other available algorithms are based on the *rational function* approach and the *geometry DIIS* approach. Also available is the *quadratic steepest descent following* method of Sun and Ruedenberg (see J. Sun and K. Ruedenberg, J. Chem. Phys. **99**, 5257 (1993)). This method is often advantageous in Transition State searches. For a detailed discussion of the various older minimization algorithms see F. Eckert, P. Pulay and H.-J. Werner, *J. Comp. Chem* **18**, 1473 (1997). Reaction path following is described in F. Eckert and H.-J. Werner, *Theor. Chem. Acc.* **100**, 21, (1998). Please refer to the references section for citations of the analytic gradient methods.

When analytical gradients are available for the optimized energy these will be used. See section analytical energy gradients for a list of methods with analytical gradients. Otherwise the gradient will be computed numerically from finite energy differences. Normally, the last computed ground-state energy is used. But the `VARIABLE`

directive or option can be used to optimize, e.g., Davidson corrected energies, excited states, or counterpoise corrected energies.

By default the program repeats in each geometry optimization step those commands in the input that are needed to compute the last energy *and* any `variable=value`

assignments in between the last energy calculation and the `OPTG`

command. For example, for MP2 gradients the commands HF and MP2 are needed. The MP2 gradients will then be computed automatically. Variables set by these programs as, e.g., `ENERGY`

, are set as usual in each optimization step. In addition, there is a variable `OPTSTEP`

, which is incremented in each optimization step. Its initial value is zero. To save the results of the initial energy calculation, one can use, e.g.,

`hf`

if(optstep.eq.0) e1=energy

optg

It is also possible to define procedures for the energy calculation, or to specify the first command from which the input should be repeated in each step (see section options to select the wavefunction and energy to be optimized). The section of the input which is needed for the geometry optimization must not modify variables that are used in the geometry definition (changes of such variables are ignored, and a warning message is printed).

## Options

Most parameters can be given as options on the `OPTG`

command line, as described in this section. Alternatively, directives can be used, which will be described in section directives for OPTG.

### Options to select the wavefunction and energy to be optimized

By default, the last computed energy is optimized, and all commands on which the last energy calculation depends are automatically executed. For certain purposes, e.g., optimization of counter-poise corrected energies or Davidson corrected energies, the following options can be used to alter the default behaviour.

Specifies a start command. In each geometry optimization step all input beginning with`STARTCMD`

=*command**command*to the current`OPTG`

is processed. This input must not include numerical gradient or Hessian calculations. If numerical gradients are needed, these will be computed for the final energy (or specified variable) by`OPTG`

. It is assumed that these commands have been executed before entering the`OPTG`

program.specifies a procedure to be executed in each geometry optimization step. This must define a complete energy calculation (orbital optimization and correlation treatment), and must not include numerical gradient of Hessian calculations (numerical gradients will be computed automatically for the optimized energy or variable). However, the procedure can include the calculation of analytical gradients, for instance for counter-poise corrected optimizations in which a linear combination of several gradient calculations is needed.`PROC`

=*procname*Optimize the value of variable`VARIABLE`

=*varname**varname*. This implies numerical gradients.

### Options for optimization methods

`METHOD`

=`PQS|RF|AH|DIIS|QSD|QSDPATH|SRMIN|SRTRANS|STSTEEP`

Optimization method to be used. See section selecting the optimization method (METHOD) for details.

Minimum search (1, default) or transition state search (2).`ROOT`

=1$|$2Determines initial step length and direction in reaction path following, see section reaction path following options (OPTION).`DIRECTION|IDIR`

=*idir*For reaction path following, see section reaction path following options (OPTION).`STPTOL`

=*value*For reaction path following, see section reaction path following options (OPTION).`SLMAX`

=*value*Max step length in one optimization step. This also affects the step length in reaction path following. For more detailed specifications see section setting a maximum step size (STEP).`STEPMAX`

=*value*(logical). If .true., the Cartesian coordinates are transformed to minimize rotations (default=.true.)`ROTATE`

The following options are not available in PQSOPT:

Trust ratio for Augmented Hessian method (default 0.5).`TRUST`

=*value*Maximum step size allowed in the Augmented Hessian procedure. This refers to the scaled parameter space (default 0.5).`AHMAX`

=*value*Threshold for orthonormalization used in conjugate gradient update of Hessian (default 1.d-3).`CUT`

=*value*

Specific options for PQSOPT (for details see options for the PQS optimizer):

This determines the coordinate type (default -1)`OPTC`

=*value*This determines whether the internal coordinates are regenerated in each step (default 1)`IGEN`

=*value*Special flag for surface/cluster optimizations (default 0)`TYPE`

=*value*Special flag for constraint optimisations (default 0)`ICONS`

=*value*

### Options for the PQS optimizer

Convergence criteria can be specified as for other methods using the `ENERGY`

, `GRADIENT`

, and `STEP`

options (see options to modify convergence criteria).

The maximum stepsize can be specified using the `STEPMAX`

option (see options for optimization methods)

Transition state optimization can be activated by setting option `ROOT=2`

.

Options for numerical gradients and Hessian options work as for other optimization programs (see options for numerical gradients and options for computing Hessians).

There are many further options, but only the following may be sometimes relevant:

`OPTC:`

`OPTC= 1`

- generate and optimize in internal coordinates; if this fails abort`OPTC=-1`

- generate and optimize in internal coordinates; if this fails switch to Cartesians and continue (default)`OPTC= 2`

- optimize in Z-Matrix internal coordinates (not supported in Molpro)`OPTC=-2`

- optimize in Z-Matrix internal coordinates (not supported in Molpro)`OPTC= 3`

- large-molecule optimization in Cartesian coordinates (never tested so far)`OPTC= 4`

- large-molecule optimization in delocalized internals (never tested so far)

`IGEN:`

`IGEN=0`

- generate a set of non-redundant internals and use throughout the optimization`IGEN=1`

- generate a new set of non-redundant internals from the same primitive space every cycle (default)`IGEN=2`

- generate a new set of underlying primitives and a new set of non-redundant internals every cycle (not recommended)

`ITYPE:`

`ITYPE=0`

- standard molecular optimization (default)`ITYPE=1`

- optimization of molecule reacting/adsorbed on model surface`ITYPE=2`

- optimization of molecular clusters (see below)

`ICONS:`

`ICONS=0`

- full optimization (no constraints, default)`ICONS=1`

- constrained optimization using Lagrange multipliers (default if constraints are given see constraint optimisations using the PQS optimizer)`ICONS=2`

- constrained optimization using penalty functions

`Cluster optimizations:`

Cluster optimizations use inverse distance coordinates for intermolecular degrees of freedom.

`MOLECULES=[`

- For cluster optimisations. The molecules must come one after the other in the geometry input.*n1,n2,n3,…*]are the number of atoms in each molecule. The sum of these must be equal to the total number of atoms. This option sets automatically`n1`

,`n2`

,…`ITYPE=2`

`CUTOFF=`

- Distance cutoff for intermolecular coordinates in Angstrom. The program will generate inverse distance coordinates for between all atoms of different molecules, unless the distance is larger than this value. Default: 5 Angstrom.*value*

`SCALE=`

- Scaling factor for inverse distance coordinates. Typical values are between 1 and 10. Default is 1.0*value*

Unfortunately, according to our current limited experience, the cluster optimization in pqsopt is neither well working in the original PQS program nor in Molpro and should therefore be used with care.

`Optimizations with z-matrix input:`

By default, optimizations with z-matrix geometry input are done with the old RF optimizer, and only the active variables in the z-matrix are optimised (cf. Defining active geometry parameters). It is possible, however, to force pqs optimization by giving the optg option `PQS`

or `METHOD=PQS`

. In this case, all cartesian coordinates will be optimised, and in each iteration the new structure is back-transformed to the z-matrix. If the z-matrix contains constants, these may change and will be overwritten. If the same variable is used in several places of the z-matrix and the corresponding new parameters are not the same, the variable is only kept in the first place, and the further ones are overwritten by the new fixed parameters. The updated z-matrix is printed in the log-file in each iteration.

### Constraint optimisations using the PQS optimizer

Constraints can be be defined in an input block (with the OPTG command block) beginning with `CONSTRAINT`

and ending with `END`

. The possible coordinate types are:

* `STRE `

(Angstrom) *I J value*

* `BEND `

(degrees) *I J K value*

* `TORS `

(degrees) *I J K L value*

* `OUTP `

(degrees) *I J K L value*

* `LINC `

(degrees) *I J K L value*

* `LINP `

(degrees) *I J K L value*

where *I, J, K, L* refer to row numbers in the XYZ input. Alternatively, atom symbols like, e.g. H1, O, H2, can be given, but these must be unique in the input geometry. The coordinate types are:

`STRE`

distance constraint between any two (different) atoms.

`BEND`

planar bend constraint between any three (different) atoms. J is the middle atom of the bend.

`TORS`

dihedral angle (proper torsion) constraint between any four (different) atoms.
The connectivity is I-J-K-L, with J-K being the central “bond”.
The angle is between the planes I-J-K and J-K-L.

`OUTP`

out-of-plane-bend constraint between any four (different) atoms. The angle is made by bond I-L with the plane J-K-L (L is the central atom)

`LINC`

colinear bending constraint between any four (different) atoms. Bending of I-J-K in the plane J-K-L.

`LINP`

perpendicular bending constraint between any four (different) atoms. Bending of I-J-K perpendicular to the plane J-K-L.

`LINC/LINP`

are special angles used when four atoms are near-linear.

The constraint values have not to be satisfied in the initial geometry. If *value* is not given, the initial value will be kept.
The number of constraints is in principle arbitrary, but currently restricted to 20.

Example: Keep the distance between atoms 2 and 3 to 1.5 ANG:

`{OPTG`

CONSTRAINT

STRE,2,3,1.5

END}

** Composite constraints:**

It is possible to define linear combinations of primitive internal coordinates as constraints. These must be defined after the simple constraints (if present). Each composite constraint begins with the keyword `COMPOSITE`

. The input of each primitive is as described above, but the `weight`

(the factor in the linear combination) replaces `value`

. The value is given at the end of the first line. Angles in the linear combinations are in radians. The general structure is:

`{OPTG`

,*options*

`CONSTRAINT`

*simple constrains* (if present)

`COMPOSITE`

*type atom_list weight value*

*type atom_list weight*

…

`END}`

*value* is the total value of the combined constraint. Note that the weights are normalised. Example:

`COMPOSITE`

STRE I J 1.0

STRE J K -1.0

BEND I J K 2.0

END

will constrain the value of $6^{-1/2}(R_{IJ} - R_{JK} + 2 A_{IJK})$ to whatever it is in the starting geometry.

Currently, the total number of composite constraints is limited to 20. Note that there is no checking for incompatible or impossible constraints, and it is the user's responsibility to make correct definitions.

### Options to modify convergence criteria

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

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

maximum number of optimization cycles. The default is 100. It cannot be increased above 100. If an optimization did not work, restart from e.g. the structure with the lowest energy using the information in the correspoing log-file (suffix .log).`MAXIT`

=*maxit*required accuracy of the optimized gradient. The default is $3 \cdot 10^{-4}$.`GRADIENT`

=*thrgrad*required accuracy of the optimized energy. The default is $1 \cdot 10^{-6}$.`ENERGY`

=*threnerg*convergence threshold for the geometry optimization step. The default is $3 \cdot 10^{-4}$.`STEP`

=*thrstep*(logical). Use Baker’s convergency criteria (see J. Baker,`BAKER`

*J. Comp. Chem.***14**,1085 (1993)).(logical). Use Gaussian convergency criteria.`GAUSSIAN`

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

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

=*thrgrms*Freeze DFT grid and domains in local calculations if the step length is smaller than`FREEZE`

=*thrfreez**thrfreez*(default 0.01).

Note: The defaults for the convergence parameters can also be changed by using a global `GTHRESH`

directive, i.e.

`GTHRESH`

, `OPTSTEP`

=*step*, `OPTGRAD`

=*grad*, `ENERGY`

=*energy*;

### Options to specify the optimization space

If the geometry is given as Z-matrix, the default is to optimize the variables on which the Z-matrix depends. In case of `xyz`

input, always all 3N coordinates are optimized, even if the `xyz`

input depends on fewer variables. If Cartesian z-matrix input is used, optimization in full space is only enforced if automatic orientation is requested using the `MASS`

, or `CHARGE`

options on the `ORIENT`

directive. See *opt_space* in section optimization coordinates (COORD) for details.

Specifies the coordinates to be used in the optimization. Z-matrix optimization is only possible if the geometry is given as Z-matrix.`SPACE=ZMAT|3N`

(logical). Same as`OPT3N|3N`

`SPACE=3N`

(logical). Same as`ZMAT`

`SPACE=ZMAT`

### Options to specify the optimization coordinates

These options specify the coordinates in which the optimization takes place. The default is to use local normal coordinates. See *opt_coord* in section optimization coordinates (COORD) for details.

`COORD=DELO[CALIZED]|NORMAL|NONORMAL|BMAT`

(logical). Same as`NORMAL`

`COORD=NORMAL`

.(logical). Same as`NONORMAL`

`COORD=NONORMAL`

.(logical). Same as`BMAT`

`COORD=BMAT`

.

`DELO[CALIZED]`

enforces full 3N optimization in delocalized internal coordinates using the PQSOPT program (available in Molpro 2023.1 or later)

`BMAT`

uses an older version of optimisations using internal coordinates.

### Options for numerical gradients

Numerical gradients can be computed with respect to variables on which the Z-matrix depends or with respect to Cartesian coordinates. In the latter case, it is most efficient to use symmetrical displacement coordinates. These do not change the symmetry of the molecule and the number of displacements is minimal. Alternatively (mainly for testing purpose) the gradients can be computed using symmetry unique Cartesian displacements or all 3N Cartesian displacements. In these cases the symmetry of the molecule can be reduced by the displacements and using such displacements is normally not recommended.

`DISPLACE=ZMAT|SYMM|UNIQUE|CART`

Displacement coordinates to be used for numerical gradient. The default is `ZMAT`

if the geometry is given as a zmatrix which depends on variables, and `SYMM`

(symmetrical displacement coordinates) otherwise. The use of `UNIQUE`

or `CART`

is not recommended.

Symmetry to be used in wavefunction calculations of numerical gradients. This option is only relevant if`SYMMETRY=AUTO|NOSYM`

`DISPLACE=UNIQUE|CART`

. If`AUTO`

is given, the maximum possible symmetry is used for each displacement. This implies that the energy is independent of the symmetry used. Note that this often not the case in MRCI or CASPT2 calculations. The option can also not be used in local correlation calculations.(logical). Same as`AUTO`

`SYMMETRY=AUTO`

(logical). Same as`NOSYM`

`SYMMETRY=NOSYM`

Step length for distances in numerical gradient calculations (in bohr). The default is 0.01.`RSTEP`

=*rstep*Step length for symmetrical displacements (in bohr). The default is 0.01.`DSTEP`

=*dstep*Step length for angles in numerical gradient calculations (in degree). The default is 1.`ASTEP`

=*astep*(logical). Use 4-point formula for accurate numerical gradient.`FOURPOINT`

(logical). Force the use of numerical gradients, even if gradients are available.`NUMERICAL`

In special cases the energy calculations may require several steps with different occupations and/or symmetries. In this case the `FORCEINP`

directive must be given (see section miscellaneous options:).

### Options for computing Hessians

By default, an approximate Hessian (model Hessian) is used. Optionally, a Hessian can be computed in the optimization or read from a previous Hessian or frequency calculation.

If given, a numerical Hessian is computed in each`NUMHESS`

=*hstep**hstep’th*iteration. If*hstep=0*or not given, only an initial Hessian is computed.Read initial Hessian from the given record. If`HESSREC`

=*record**record*is not given or zero, the last computed Hessian is used.(logical). Same as`READHESS`

`HESSREC=0`

.specifies a procedure to be used for computing the Hessian. This procedure must be define a complete energy calculation (orbital optimization and correlation treatment). A different method can be used than for the optimized energy, but the basis must not be redefined in this procedure. For instance, an MP2 Hessian can be used for CCSD(T) optimizations, or a CASPT2 Hessian for MRCI optimizations. By default, the same procedure is used for the Hessian as for the optimized energy. Note: If a hessian procedure is used and two or more optimizations are done after each other in the same input, the complete energy calculation (including orbital optimization) must be defined before each optg command or in the optg procedure.`HESSPROC`

=*procname*Compute Hessian for variable`HESSVAR`

=*varname**varname*. This implies numerical calculation of the Hessian from energies. The default is to use the same variable as for the energy and gradient.Use central gradient differences for computing Hessian (only effective if gradients are available)`HESSCENT`

Use forward gradient differences for computing Hessian (only effective if gradients are available). This effectively computes the Hessian at a slightly displaced geometry, but needs only half the number of displacements. This is the default.`HESSFORW`

`UPDATE`

=`BFGS|IBFGS|CGRD|PMS|POWELL|MS|NONE`

Hessian update method to be used. See section hessian update (UPDATE) for details.

Max number of Hessian updates. The count is reset to zero each time a Hessian is computed.`MAXUPD`

=*maxupd*If true, replace diagonal elements of model hessian by diagonal numerical hessian (if available). This sometimes improves convergence, but since it may lead to symmetry breaking it is no the default.`NUMDIAG`

Note that there are restrictions for computing Hessians for multireference methods (MCSCF, MRCI, ACPF,AQCC,RS2). For these methods the symmetry must not change by any displacements, since this could change the occupations and states and may lead to non-contiguous potential energy surfaces. One of the following three options can be used in these cases:

- Use no symmetry from the beginning (
`NOSYM`

). - Use symmetric displacement coordinates. This is the default if the optimization is done in 3N cartesian coordinates. One can use
`OPTG,DISPLACE=SYMM`

to force the use of symmetrical displacements (this creates 3N cartesian coordinates if a Z-matrix is used in the geometry input). - Use a Z-matrix with the restriction that no variable in the Z-matrix may change the symmetry. For example,
`geometry={O;H1,O,r;H2,O,r,H1,theta}`

would work, but`geometry={O;H1,O,r1;H2,O,r2,H1,thetai}`

would not work. In this case the program prints a warning message. If an incorrect Z-matrix is used and the symmetry changes, the program will crash.

### Miscellaneous options:

Save Cartesian gradients in variables`VARSAVE`

`GRADX, GRADY, GRADZ`

.Do not compute gradients at lattice points.`NONUC`

Set debug print options.`DEBUG`

Print option for optimization.`PRINT`

=*iprint*Save the optimized coordinates in an xyz-file. One file is written for each step. The filename is`SAVEXYZ`

[=*file*]*file*_*nn*`.xyz`

, where*nn*is the iteration number. For the final geometry, _*nn*is omitted. If filename is not given,*file*is taken to be the root name of the input, i.e.`test.inp`

creates`test_1.xyz`

in the first iteration and`test.xyz`

for the converged geometry. By default, the xyz information is written to the log file in each step.Save optimized variables in each step. The file name is`SAVEACT`

[=*file*]*file*`.act`

. If*file*is not given the root name of the input is used. The file can be read later using`the READVAR`

command or copied into new input.Write in each step the Cartesian coordinates and gradients. The file name is`SAVEGRD`

[=*file*]*file*`.grd`

. If file is not given, the root name of the input appended by`.grd`

is used.(logical). If given, existing`APPEND`

`SAVEACT`

and/or`SAVEGRD`

files are appended.(logical). If given, the`REWIND`

`SAVEACT`

and/or`SAVEGRD`

files are rewound at each step, i.e. only the last geometry or gradient is saved, previous values are overwritten.(logical) Disables the mechanism that freezes the occupations during geometry optimizations. If this option is given, the occupations and wave function definitions given in the input for the energy calculations are strictly obeyed.`FORCEINP`

## Directives for OPTG

An alternative way to specify options is to use directives, as described in this section. In some cases this allows more detailed specifications than with the options on the `OPTG`

command. In particular, directives `ACTIVE`

or `INACTICE`

can be used to define the optimization space in more detail.

### Selecting the optimization method (METHOD)

`METHOD`

,*key*;

*key* defines the optimization method.

For *minimization* the following options are valid for *key*:

PQS optimizer (default for full (3N) optimisations).`PQS`

Rational Function method (default for Z-matrix optimisations).`RF`

Augmented Hessian method. This is similar to`AH`

`RF`

algorithm but uses a more sophisticated step restriction algorithm.Pulay’s Geometry DIIS method. As an an additional option you may add the number of geometries to be used in GDIIS interpolation (default 5) and the interpolation type (i.e. the subspace in which the GDIIS interpolation is made.`DIIS`

`METHOD,DIIS,`

*number*, *type*

*type* may be `GRAD`

interpolation using the gradients (default), working good for rigid molecules, `STEP`

interpolation using Quasi-Newton steps which could be advantageous in dealing with very floppy molecules, `ENER`

interpolation using energies, which is an intermediate between the above two.

Quadratic steepest descent method of Sun and Ruedenberg.`QSD`

Old version of`SRMIN`

`QSD`

.

For *transition state* searches (invoked with the `ROOT`

option, see section transition state (saddle point) optimization (ROOT)) *key* can be

PQS optimizer using internal coordinates (default for full (3N optimisations).`PQS`

Rational Function method (default for Z-matrix optimisations).`RF`

Pulay’s Geometry DIIS method (see above).`DIIS`

Quadratic Steepest Descent Transition State search using the image Hessian method (see J. Sun and K. Ruedenberg, J. Chem. Phys.`QSD`

**101**, 2157 (1994)) The use of this option is recommended for transition state searches – especially in complicated cases. The optimization step is checked and the Hessian is recalculated when approaching a troublesome region of the PES. Thus**this method is somewhat safer (and often faster) in reaching convergence than the RF or DIIS method**. The Hessian recalculation safeguard may be turned off using the`METHOD,QSD,NOHESS`

input card.Old version of`SRTRANS`

`QSD`

.

For *reaction path following* the input *key* is

Quadratic Steepest Descent reaction path following. This methods determines reaction paths (intrinsic reaction coordinates, IRCs) by following the exact steepest descent lines of subsequent quadratic approximations to the potential energy surface. The Hessian matrix is calculated numerically at the first optimization step and subsequently updated by Powell or BFGS update. If a given arc length of the steepest descent lines is exceeded, the Hessian is recalculated numerically (see`QSDPATH`

`OPTION`

section reaction path following options (OPTION)). For details see J. Sun and K. Ruedenberg, J. Chem. Phys.**99**, 5269 (1993) It is also possible to recalculate the Hessian after each*m*steps using the`NUMHES`

,*m*command (see section numerical Hessian (NUMHESS)). If the Hessian matrix is recalculated in every optimization step (`NUMHES`

,$1$) a algorithm different to the one with updated Hessians is used, which is very accurate. Using the`PRINT,OPT`

card, this algorithm prints in every optimization step a*reaction path point r*which is different from the point where the energy and the gradient is calculated but closer to the real reaction path (for further details of the algorithm see J. Sun and K. Ruedenberg, J. Chem. Phys.**99**, 5257 (1993)). For further input options of the QSD reaction path following see`OPTION`

section reaction path following options (OPTION).Old Version of`SRSTEEP`

`QSDPATH`

.

### Optimization coordinates (COORD)

It is possible to use various coordinate types and algorithms for the optimization. This can be controlled by additional subcommands as described in this and the following subsections.

`COORD`

,[*opt_space*],[*opt_coord*],[`NOROT`

]

These options choose the optimization space and the coordinate system in which the optimization takes place.

*opt_space* defines the parameters to be optimized. By default, if the geometry input is given in Z-matrix format, all variables on which the Z-matrix depends are optimized. Subsets of the variables on which the Z-matrix depends can be chosen using the `ACTIVE`

or `INACTIVE`

subdirectives. If the Z-matrix depends on no variables or `xyz`

input is used, all $3N$ cartesian coordinates are optimized.

*opt_space* can be one of the following:

Optimize all variables on which the Z-matrix depends (default if the geometry is given as Z-matrix).`ZMAT`

Optimize all $3N$ cartesian coordinates (default if the Z-matrix depends on no variables, or if xyz-input is used). Z-Matrix input coordinates will be destroyed if 3N is used.`3N`

*opt_coord* determines the coordinates in which the optimization takes place. Optionally, delocalized internal coordinates, local normal coordinates, cartesian coordinates, or natural internal coordinates can also be used.

*opt_coord* can be one of the following:

Optimization in delocalized internal coordinates using PQSOPT (available in Molpro 2023.1 or later). This is the default for full space (3N) optimisations. If the option is given, 3N optimization is enforced (but constraints can be defined).`DELO[CALIZED]`

Optimization in local normal coordinates.`NORMAL`

Don’t use local normal coordinates.`NONORM`

Use Pulay’s`BMAT`

[=*filename*]*natural internal coordinates*, see G. Fogarasi, X. Zhou, P. W. Taylor and P. Pulay*J. Am. Chem. Soc.***114**, 8191 (1992); P. Pulay, G. Fogarasi, F. Pang, J. E. Boggs*J. Am. Chem. Soc.***101**, 2550 (1979)). Optionally, the created coordinates as well as additional informations about this optimization are written to the specified file. These coordinates resemble in part the valence coordinates used by vibrational spectroscopist, and have the advantage of decreasing coupling between different modes. This often increases the speed of convergence. Similar coordinates are also used in the PQSOPT program (`DELO`

), and mostly this converges even better.

If the option [`NOROT`

] is given, the cartesian coordinates are not transformed to minimize rotations.

### Displacement coordinates (DISPLACE)

`DISPLACE`

,*displacement_type*

see section choice of coordinates (COORD) for details.

### Defining active geometry parameters (ACTIVE)

`ACTIVE`

,*param*;

Declares variable name *param* to be active in the optimization. By default, initially all variables on which the geometry depends are active; inclusion of an `ACTIVE`

card makes all parameters inactive unless explicitly declared active (see also `INACTIVE`

).

### Defining inactive geometry parameters (INACTIVE)

`INACTIVE`

,*param*;

Declares variable name *param* to be inactive in the optimization. If any `ACTIVE`

card appears in the input, this card is ignored! (see also `ACTIVE`

)

### Hessian approximations (HESSIAN)

By default, the Molpro geometry optimization utilizes a force field approximation to the hessian (“Model Hessian”, see R. Lindh, A. Bernhardsson, G. Karlström and P. Malmqvist Chem. Phys. Lett. **241**, 423 (1995)), which speeds up convergence significantly. The Model Hessian is parameterized for the elements up to the third row. Alternatively, the model Hessian of Schlegel can be used, or the Hessian can be computed numerically (see also section numerical Hessian (NUMHESS)).

`HESSIAN`

,*options*

where *options* can be

Use Lindh’s Model Hessian in optimization (default).`MODEL`

Use Schlegel’s Model Hessian.`MODEL=SCHLEGEL`

Add vdW terms to Lindh’s Model Hessian.`MODEL=VDW`

same as`SCHLEGEL`

`MODEL=SCHLEGEL`

.same as`VDW`

`MODEL=VDW`

.Don’t use Model Hessian approximation to the hessian.`NOMODEL`

Recompute Hessian after`NUMERICAL`

=*hstep**hstep*iterations. This disables the use of a model hessian. If*hstep*=0, the Hessian is only computed in the first iteration. Default parameters are used for computing the numerical Hessian, unless modified using options as described for the`NUMHESS`

directive, see Sect. numerical Hessian (NUMHESS). Any option valid for the`NUMHESS`

directive may also follow the`NUMERICAL`

option on the`HESSIAN`

directive.Read Hessian from given`READ|RECORD|HESSREC`

=*record**record*. If*record*is not given or zero, the last computed hessian will be read. See section numerical Hessian (NUMHESS) for more details about numerical Hessians.Method used for hessian update. See section hessian update (UPDATE) for possibilities and details.`UPDATE`

=*type*Max number of hessian updates. The count is reset to zero each time a hessian is computed.`MAXUPD`

=*maxupd*

If the Model Hessian is disabled (`NOMODEL`

) and no Hessian is read or computed, the initial hessian is assumed to be diagonal, with values 1 hartree*bohr**(-2) for all lengths, 1 hartree*radian**(-2) for all angles. Additional matrix elements of the hessian can be defined using the HESSELEM directive, see section hessian elements (HESSELEM).

In transition state searches the Hessian is evaluated numerically in the first iteration by default. Alternatively, if `READ`

is specified, a previously computed hessian is used.

### Numerical Hessian (NUMHESS)

`NUMHESS`

,*options*

or

`NUMHESS`

,*hstep,options* If this directive is present a numerical Hessian is computed using finite differences. If analytical gradients are available, one can use forward gradient differences (needs one gradient calculation for each coordinate) or central differences (more accurate, needs two gradient calculations for each coordinate). For transition state optimizations it is usually sufficient to use forward differences. If analytical gradients are not available for the optimized method, the energy is differentiated twice. In this case only central differences are possible.

The following options can be given:

`HSTEP`

=*hstep*
*hstep*=-1: Don’t calculate numerical hessian (default for minimization);

*hstep*=0 Calculate numerical hessian only once at the start of the optimization (default for transition state searches).

*hstep*=*n* Calculate numerical hessian after each *n* optimization steps. This is useful for difficult transition state optimizations (e.g. if the eigenvalue structure of the hessian changes during the optimization).
`FORWARD`

Use forward differences (default).
`CENTRAL`

Use the more accurate central differences.
`RSTEP`

=*rstep*
Step length for distances (in bohr). The default is 0.01.
`ASTEP`

=*astep*
Step length for angles (in degree). The default is 0.5 or 1 for angles below and above 90 degree, respectively.
`DSTEP`

=*dstep*
Step length for symmetrical displacements (in bohr). The default is 0.01.
`VARIABLE`

=*varname*
Use given variable for numerical calculation of the Hessian. Note that this disables the use of gradients, and Hessian evaluation can be very expensive.
`PROCEDURE`

=*procname*
Procedure to be used for computing Hessian. This procedure must be define a complete energy calculation (orbital optimization and correlation treatment). A different method can be used than for the optimized energy. For instance, an MP2 hessian can be used for CCSD(T) optimizations, or a CASPT2 hessian for MRCI optimizations. By default, the same procedure is used for the hessian as for the optimized energy.
`DISPLACE`

=*type*
*type* can be one of the following:

`SYMM`

Use symmetric displacement coordinates (default). This is the only recommended option.
`CART`

Use $3N$ cartesian displacements (not recommended). This requires many more energy calculations than necessary and does not preserve the molecular symmetry.
`UNIQUE`

Use symmetry-unique cartesian displacements (not recommended)

Note that the displacement type for gradient and hessian must be the same.
`CALC`

=*icalc*
*icalc*=0: Recalculate the complete Hessian matrix numerically after each *hstep* optimization steps (default).

*icalc*=1 (currently disabled): Recalculate selected Hessian matrix elements if the relative deviation of this element before and after update (see `UPDATE`

, section hessian update (UPDATE)) is larger than *thresh*. If *thresh* is not specified, a default value of $thresh=0.05$ (i.e. a maximum deviation of $5 \%$) is used.

*icalc*=2 (currently disabled): Recalculate complete Hessian matrix if the RMS deviation of the Hessian matrix before and after update is larger than $thresh$. If *thresh* is not specified a default value of $thresh=0.5$ $a.u.$ is used.
`THRESH`

=*thresh*
Threshold for partial or dynamical update of hessian, see above

### Hessian elements (HESSELEM)

[optgeo:hesselem]

`HESSELEM`

,*value*, *active1*,*active2*,…

sets the starting value for hessian matrix element between active variables *active1*, *active2* to *value*. If *active2* is omitted it defaults to *active1* (diagonal element). As many `HESSELEM`

directives as needed may be given.

### Hessian update (UPDATE)

`UPDATE`

,[`TYPE`

=]*type*,`MAXUPD`

=*maxupd*

This directive chooses the update type and limits the number of points used for the hessian update to *maxupd*. The default number of steps used in hessian update procedures is 5. If there are symmetry constraint in the coordinates of the optimization, the default number may be lower than five.

In minimizations *type* may be

Use BFGS update of hessian (default).`BFGS`

Use BFGS update of the inverse hessian.`IBFGS`

Use Conjugate Gradient update (see also`CGRD`

`CUT`

,`TRUST`

).Don’t do any update.`NONE`

In transition state optimizations *type* may be

Combined Powell/Murtagh-Sargent update of hessian (default).`PMS`

Use Powell’s update of the hessian.`POWELL`

Use update procedure of Murtagh and Sargent.`MS`

Don’t do any update.`NONE`

### Numerical gradients (NUMERICAL)

`NUMERICAL`

,*options*,*active$_1$=step$_1$, active$_2$=step$_2$ …*;

With this directive the gradients are computed by finite differences. *step$_i$* is the increment for the active geometry parameter $active_i$. For active parameters which are not specified, the default values are used. By default, the increment is 0.01 bohr for bond distances and 0.5 or 1 degree for angles less than or greater than 90 degrees, respectively. These defaults can be modified by specifying `RSTEP`

or `ASTEP`

. `DSTEP`

is the length of symmetrical displacements, which are used if the optimization is performed in 3N coordinates.

For each active variable, two energy calculations are necessary in each geometry optimization step – so numerical optimizations may be expensive! In optimizations of 3N coordinates symmetrical displacement coordinates are normally used to minimize the number of energy calculations. (see section choice of coordinates (COORD)).

For optimization of special energies see `VARIABLE`

section optimizing energy variables (VARIABLE).

The following options can be given:

Step length for distances (in bohr). The default is 0.01.`RSTEP`

=*rstep*Step length for angles (in degree). The default is 0.5 or 1 for angles below and above 90 degree, respectively.`ASTEP`

=*astep*Step length for symmetrical displacements (in bohr). The default is 0.01.`DSTEP`

=*dstep*Use central differences for gradient (default)`CENTRAL`

Use forward differences (not recommended for gradient).`FORWARD`

Use four-point formula for very accurate numerical gradients.`FOURPOINT`

Use given procedure for numerical calculation of the gradient. This procedure must define a complete energy calculation (orbital optimization and correlation treatment).`PROCEDURE`

=*procname*Use given variable for numerical calculation of the gradient.`VARIABLE`

=*varname*The displacement type. Note that the displacement type for gradient and hessian must be the same.`DISPLACE`

=*type**type*can be one of the following:Use symmetric displacement coordinates (default). This is the only recommended option.`SYMM`

Use $3N$ cartesian displacements (not recommended). This requires many more energy calculations than necessary and does not preserve the molecular symmetry.`CART`

Use symmetry-unique cartesian displacements (not recommended)`UNIQUE`

### Transition state (saddle point) optimization (ROOT)

`ROOT`

,*root* [root]

Specifies the eigenvector of the hessian to be followed.

specifies a minimization (default).*root*=1specifies a transition state (saddle point) optimization.*root*=2

In the present implementation a saddle point search is possible with the rational function method (`METHOD,RF`

), the geometry DIIS method (`METHOD,DIIS`

) and the quadratic steepest descent method of Sun and Ruedenberg (`METHOD,SRTRANS`

).

Note that convergence is usually much more difficult to achieve than for minimizations. In particular, a good starting geometry and a good approximation to the hessian is needed. The latter is achieved by evaluating the hessian numerically (see section numerical Hessian (NUMHESS)) or using a precomputed hessian (see section hessian approximations (HESSIAN)).

Example:

- examples/butane_opt_transition.inp
***, Bicyclo-[1.1.0]-butane Transition State basis=3-21G L1=1.495 ang ! Define Active Variables L2=1.418 ang L3=1.463 ang L4=1.093 ang L5=1.111 ang L6=1.098 ang L7=1.097 ang L8=1.110 ang L9=1.106 ang A1=92.1 degree A2=62.1 degree A3=136.0 degree A4=123.5 degree A5=122.4 degree A6=124.7 degree A7=126.7 degree A8=117.9 degree D1=-120.4 degree D2=4.4 degree D3=108.8 degree D4=-107.5 degree D5=84.2 degree D6=109.3 degree D7=-106.1 degree geometry={C1 ! Geometry Specification Z-Matrix C2 1 L1 C3 2 L2 1 A1 C4 1 L3 2 A2 3 D1 H5 1 L4 2 A3 3 D2 H6 2 L5 1 A4 4 D3 H7 3 L6 2 A5 1 D4 H8 3 L7 2 A6 1 D5 H9 4 L8 1 A7 2 D6 H10 4 L9 1 A8 2 D7} int rhf optg root,2 ! Transition State search method,qsd ! Use Quadratic Steepest Descent Method

### Setting a maximum step size (STEP)

`STEP`

,*steplength*,*drmax*,*damax*,*drmax1*,*damax1* [optgeo:step]

*steplength* is the initial step length in the scaled parameter space (default 0.3). In the AH-method this is dynamically adjusted, and can have a maximum value *ahmax* (see `TRUST`

).

*drmax* is the initial max change of distances (in bohr, default 0.3). In the AH-method this is dynamically adjusted up to a maximum value of *drmax1* (default 0.5 bohr).

*damax* is the initial max change of angles (in degree, default 2). In the AH-method this is dynamically adjusted up to a maximum value of *damax1* (default 10 degrees).

### Redefining the trust ratio (TRUST)

`TRUST`

,*ratio*,*ahmax*

*ratio* determines the radius around the current minimum in which points are used to update the Hessian with the conjugate gradient method (default 0.5; see also `UPDATE`

).

*ahmax* is the maximum step size allowed in the Augmented Hessian procedure. This refers to the scaled parameter space (default 0.5). The initial step size is *stepmx* (see `STEP`

card).

### Setting a cut parameter (CUT)

`CUT`

,*threshold*

Specifies a threshold for ortho-normalization used in conjugate gradient update of hessian (default 1.d-3; see also `UPDATE`

).

### Line searching (LINESEARCH)

`LINESEARCH`

,*iflag*,*thrlmin*,*thrlmax*

Interpolate the geometry of the stationary point (minimum or saddle point) by a quartic polynomial between the current and the previous geometry. If *iflag*=0 or no *iflag* is set, the next optimization step will be taken from the interpolated geometry using the interpolated energy and gradient. If *iflag*=1 the energy and gradient will be recalculated at the interpolated geometry before taking the new optimization step. Note though, that the additional effort of recalculating the energy and gradient is usually not met by the increase of the convergence rate of the optimization. *thrlmin* and *thrlmax* are min and max thresholds for the recalculation of the energy and the gradient in case *iflag*=1. I.e. the recalculation just takes place if the interpolated geometry isn’t too close to the actual geometry *thrlmin* and isn’t too remote from the actual geometry *thrlmax*. Default values are *thrlmin*=0.001 and *thrlmax*=0.05 in the scaled parameter space of the optimization.

### Reaction path following options (OPTION)

[optgeo:option] `OPTION`

,*key=param*;

where *key* can be

If starting at a transition state (or near a transition state) determine where to take the first step. If`IDIR`

`IDIR=0`

is chosen, the first step will be towards the transition state. This is the default. If`IDIR=1`

is given in the input the first optimization step will be along the ”transition vector” i.e. the hessian eigenvector to the smallest eigenvalue which points down towards the minimum. If using a larger`IDIR`

parameter, the first step will be larger; if using a negative value, the first step will be in the opposite direction.If using an updated hessian matrix, this parameter determines what update to take. If the step size between two subsequent points on which the steepest decent lines are puzzled together is smaller than`STPTOL`

*stptol*(i.e. if we are close to a minimum) the BFGS update is used, otherwise it is Powell update. The default value of*stptol*is $1.d-6$. Note that the stepsize is also affected by the`STEPMAX`

option.This option is only valid with the old version of the reaction path following algorithm (i.e.`SLMAX`

`METHOD,SRSTEEP`

). In this algorithm`slmax`

determines the frequency of the recalculation of the numerical hessian. If the total step size of the last steps exceeds*slmax*the hessian will be recalculated, otherwise it will be updated. By default*slmax*is two times the maximum step size of the optimization step*steplength*(see`STEP`

section setting a maximum step size (STEP)). If you are using`METHOD,QSD`

, the`SLMAX`

option is obsolete and the`NUMHES`

command (see above) should be used instead.

These options can also be given on the `OPTG`

command line.

### Optimizing energy variables (VARIABLE)

`VARIABLE`

,*name*; [optgeo:var]

Defines a variable *name* which holds the energy value to be optimized in using finite differences. By default, this is `ENERGY(1)`

as set by the most recent program. Other variables which can be used are

holds last energy for state`ENERGY`

(*i*)*i*.holds last reference energy for state`ENERGR`

(*i*)*i*.holds last Davidson corrected energy for state`ENERGD`

(*i*)*i*.holds last Pople corrected energy for state`ENERGP`

(*i*)*i*.holds CCSD (QCI, BCCD) energy in CCSD(T) [QCI(T), BCCD(T)] calculations (single state optimization).`ENERGC`

holds CCSD(T) energy in CCSD(T) calculations (single state)`ENERGT(1)`

holds CCSD[T] energy in CCSD(T) calculations (single state).`ENERGT(2)`

holds CCSD-T energy in CCSD(T) calculations (single state).`ENERGT(3)`

These variables are set automatically by the CI and/or CCSD programs. It is the user’s responsibility to use the correct variable name; an error exit occurs if the specified variable has not been defined by the last program or the user.

Note: The use of the `VARIABLE`

option triggers `NUMERICAL`

, so optimization can be very inefficient!

### Printing options (PRINT)

`PRINT`

,*code=level,…*;

Enables printing options. Usually level should be omitted or 0; values of $level>0$ produce output useful only for debugging. *code* can be

prints the updated hessian matrix. Note that its diagonal elements are printed anyway.`HESSIAN`

prints the complete set of previous geometries, gradients and energies.`HISTORY`

prints extended gradient information`GRADIENT`

prints detailed information about the optimization process (mainly for debugging).`OPT`

Several print options can be specified with one `PRINT`

command.

### Conical Intersection optimization (CONICAL)

[conical]

To optimize a conical intersection between two electronic states having the same spin, three vectors must be evaluated at SA-CPMCSCF level:

- Non-Adiabatic Derivative Coupling (DC).
- Gradient of the lower state (LSG).
- Gradient of the upper state (USG).

This requires three different `CPMCSCF`

directives in the MULTI input:

`CPMCSCF, NACM`

, $S_i$, $S_j$, `ACCU`

=*1.0d-7*, `record`

=*record1.file*

`CPMCSCF, GRAD`

, $S_i$, `SPIN`

=*Spin of state $S_i$*, `ACCU`

=*1.0d-7*, `record`

=*record2.file*

`CPMCSCF, GRAD`

, $S_j$, `SPIN`

=*Spin of state $S_j$*, `ACCU`

=*1.0d-7*, `record`

=*record3.file*

where $S_i$,$S_j$ are the electronic states in the usual format *istate.istsym*, and *record[n].file* specifies the name and the file number where CPMCSCF solutions should be stored. Parameter `SPIN`

is half of the value in the `WF`

card used to define the electronic state.

Things to remember:

- Specify always three different
*record.file*on the`CPMCSCF`

directives. - Evaluate the CPMCSCF for USG always last.
- Skip the DC evaluation if the conical intersection involves states with different spin (e.g., a Singlet/Triplet crossing) because the coupling is then zero.

Three sets of `FORCE`

commands (only two for Singlet/Triplet intersection) follow the `MULTI`

input. They will be like:

`FORCE`

`SAMC`

,*record[n].file*

`CONICAL`

,*record4.file*[,`NODC]`

where *record.file* is one of the records containing CPMCSCF info and *record4.file* points to a free record used for internal storage by the CONICAL code. *record4.file* must be the same on all the `CONICAL`

directives. Furthermore, the present implementation works properly only if *file=1* on the `CONICAL`

directive. The optional keyword `NODC`

must be used in case of different spins (e.g., S/T crossing) when DC is not needed.

The actual optimization is performed using `OPTG,STARTCMD=MULTI`

The example below optimizes the conical intersection in $LiH_2$ (ground and excited states are both doublets).

- examples/lih2_D0D1.inp
***, LiH2 basis=sto-3g print,orbitals,civector symmetry,x !use only molecular plane. Both states must be in the same symmetry. geometry={ Li; h1,Li,r; h2,Li,r,h1,theta} r=3.0 theta=35 {hf;wf,4,1,0} {multi;occ,6,1;wf,5,1,1;state,2 !state averaged casscf CPMCSCF,NACM,1.1,2.1,accu=1.0d-7,record=5100.1 !cpmcscf for non-adiabatic couplings CPMCSCF,GRAD,1.1,spin=0.5,accu=1.0d-7,record=5101.1 !gradient for state 1 CPMCSCF,GRAD,2.1,spin=0.5,accu=1.0d-7,record=5102.1} !gradient for state 2 {Force SAMC,5100.1 !compute coupling matrix element CONICAL,6100.1} !save information for optimization of conical intersection {Force SAMC,5101.1 !compute gradient for state 1 CONICAL,6100.1} !save information for optimization of conical intersection {Force SAMC,5102.1 !compute gradient for state 2 CONICAL,6100.1} !save information for optimization of conical intersection optg,startcmd=multi !find conical intersection

This second example optimizes the singlet-triplet intersection in $LiH_2(+)$ (ground state is Singlet, excited state is Triplet).

- examples/lih2+_S0T0.inp
***, LiH2 basis=sto-3g symmetry,nosym geometry={ Li; H1,Li,r; H2,Li,r,H1,theta} r=3.7 theta=160 {hf;wf,4,1,0} {multi; occ,7; wf,4,1,0; !singlet state wf,4,1,2; !triplet state CPMCSCF,GRAD,1.1,spin=0,accu=1.0d-7,record=5101.1 !cpmcscf for gradient of singlet state CPMCSCF,GRAD,1.1,spin=1,accu=1.0d-7,record=5100.1 !cpmcscf for gradient of triplet state } {Force SAMC,5101.1 !state averaged gradient for singlet state CONICAL,6100.1,NODC} !save information for OPTCONICAL {Force SAMC,5100.1 !state averaged gradient for triplet state CONICAL,6100.1,NODC} !save information for OPTCONICAL optg,startcmd=multi,gradient=1.d-6 !find singlet-triplet crossing point

## Using the SLAPAF program for geometry optimization

It is optionally possible to use the `SLAPAF`

program written by Roland Lindh for geometry optimizations. This is done by prepending the optimization method with ’SL’. The following methods are supported:

Use the rational function approximation;`SLRF`

Use the Newton-Raphson method;`SLNR`

Use the C1-DIIS method;`SLC1`

Use the C2-DIIS method.`SLC2`

When using DIIS methods (`SLC1`

or `SLC2`

), the DIIS parameters are specified in the same way as in standard molpro optimizer.

There are some differences when using the `SLAPAF`

program:

- It is not possible to use Z-matrix coordinates in the optimization.
- Instead, one can explicitly define internal coordinates to be varied or fixed.
- Additional constraints can be imposed on the converged geometry in a flexible way.

### Defining constraints

Constraints and internal coordinates (see below) can be linear combinations of bonds, angles etc. The latter, called here primitive internal coordinates, can be specified before the constraints definition, or directly inside. The general definition of a primitive coordinate is:

`PRIMITIVE`

,[`NAME=]`

*symbolic name*, *explicit definition;*

or

`PRIM`

,[`NAME=]`

*symbolic name*, *explicit definition;*

Here *symbolic name* is the name given to the primitive coordinate (if omitted, it will be generated automatically). This name is needed for further reference of this primitive coordinate.

*explicit definition* has the form:

*type,atoms*

*type* can be one of the following:

Bond length, defined by 2 atoms.`BOND`

Bond angle, defined by 3 atoms (angle 1–2–3).`ANGLE`

Dihedral angle, defined by 4 atoms (angle between the planes formed by atoms 1,2,3 and 2,3,4, respectively).`DIHEDRAL`

Out-of-plane angle, defined by 4 atoms (angle between the plane formed by atoms 2,3,4 and the bond 1–4).`OUTOFPLANE`

A dissociation coordinate, defined by two groups of atoms (Not permitted with constraints).`DISSOC`

Cartesian coordinates of an atom.`CARTESIAN`

For all types except `DISSOC`

and `CARTESIAN`

, atoms are given as:

`ATOMS=`

[*a1,a2,a3,…*]

where the number of atoms required varies with *type* as specified above, and the atomic names *a1,a2,a3,…* can be either atomic tag names from the Z-matrix input, or integers corresponding to Z-matrix rows. Note that the square brackets are required here and do not indicate optional input.

For `DISSOC`

the specification is as follows:

`DISSOC,GROUP1=`

[*a1,a2,…*],`GROUP2=`

[*b1,b2,…*];

The corresponding internal coordinate is the distance between the centres of mass of the two groups.

For `CARTESIAN`

the definition is

`CARTESIAN`

, *I, atom;*

where *I* can be one of `X,Y,Z`

or 1,2,3 and *atom* can be a z-matrix atom name or an integer referring to the z-matrix row.

With this definition, the constraints are defined as

`CONSTRAINT`

,[`VALUE=`

]*value*,[*unit*],[[`FACTOR=`

]*fac*,*prim*,[[`FACTOR=`

]*fac*],*prim,…;*

where *value* is the value imposed to the constraint, and *prim* is either the name of the primitive defined before this constraint, or an explicit definition; and *fac* is a factor of the corresponding primitive in the constraint. If *fac* is omitted it is taken to be 1.

If *value* is specified in `Angstrom`

or `Radian`

, *unit* must be given.

Examples for H$_2$O in $C_s$ symmetry:

Constraining the bond angle to 100 degrees:

`constraint,100,deg,angle,atoms=[h1,o,h2];`

which is equivalent to

`primitive,a1,angle,atoms=[h1,o,h2];`

`constraint,100,a1;`

Keeping the two OH distances equal:

`constraint,0,bond,atoms=[h1,o],-1.,bond,atoms=[h2,o];`

which is equivalent to

`primitive,b1,bond,atoms=[h1,o];`

`primitive,b2,bond,atoms=[h2,o];`

`constraint,0,b1,-1.,b2;`

### Defining internal coordinates

By default `SLAPAF`

optimizes in force-constant weighted normal coordinates that are determined automatically. However, the user can define his own coordinates. The definition of internal coordinates, similar to constraints, is based on primitive coordinates. The input is:

`INTERNAL`

,[[`NAME=`

]*name*],[[`FACTOR=`

]*fac*],*prim*,[[`FACTOR=`

]*fac*],*prim,…;*

`FIX`

, [[`NAME=`

]*name*],[[`FACTOR=`

]*fac*],*prim*,[[`FACTOR=`

]*fac*],*prim,…;*

Internal coordinates that are specified using `INTERNAL`

are varied and those using `FIX`

are fixed to their initial values.

An important point for the definition of internal coordinates is that their total number must be equal to the number of degrees of freedom of the molecule . Otherwise an error message is generated. Only symmetry independent coordinates need to be given.

### Additional options for SLAPAF

Some options can be passed to the SLAPAF program. Options are specified with SLOPT subdirective:

`{opt;method=slnr;{slopt;opt1;opt2,`

*par1,par2*`;opt3;...}}`

The available options are

Use eigenvectors of the approximate Hessian, expressed in cartesian coordinates, as the definition of internal coordinates;`CART`

Don’t impose any restrictions on the step size;`NOMA`

Order the gradients and displacement vectors according to Schlegel prior to the update of the Hessian. Default is no reordering;`UORD`

Use force field weighted internal coordinates (default);`HWRS`

Activate RS-P-RFO as default for transition state search; default is RS-I-RFO;`RS-P`

Use unweighted internal coordinates;`NOHW`

Print B-matrix;`PRBM`

Thresholds for redundant coordinate selection for bonds, bends and torsions, respectively. Default 0.2, 0.2, 0.2`RTHR,`

*Thra,Thrb,Thrt*Hessian vector index for mode following when calculating transition states.`MODE,`

*index*Enable unconstrained optimization for constrained cases, when looking for transition states (see MOLCAS manual).`FIND`

Threshold for FIND, default 0.2 (see MOLCAS manual).`GNRM,`

*thr*

For more information, please consult the MOLCAS manual.

## Examples

- examples/h2o_optmp2.inp
***,H2O !title include procedures basis=vtz !use cc-pVTZ basis geometry={O; !Z-matrix for water H1,O,R; H2,O,R,H1,THETA} R=0.96 Ang !start bond distance Theta=104 !start bond angle optmp2 !do MP2 geometry optimization show,R,Theta !show optimized geometry parameters

- examples/h2o_optmp2_runccsdt.inp
include procedures geometry={O; !Z-matrix for water H1,O,R; H2,O,R,H1,THETA} basis=vtz !use VTZ basis R=0.96 Ang !start bond distance Theta=104 !start bond angle hf mp2 optg !optimize MP2 energy mp4 !do single-point MP4 at MP2 geometry ccsd(t) !do single-point CCSD(T) at MP2 geometry

There is an additional example in the introductory examples.

### Simple HF optimization using Z-matrix

- examples/allene_optscf.inp
***, Allene geometry optimization using Z-Matrix basis=sto-3g rcc=1.32 ang rch=1.08 ang acc=120 degree Geometry={C1 !Z-matrix input C2,c1,rcc Q1,c1,rcc,c2,45 C3,c2,rcc,c1,180,q1,0 h1,c1,rch,c2,acc,q1,0 h2,c1,rch,c2,acc,h1,180 h3,c3,rch,c2,acc,h1,90 h4,c3,rch,c2,acc,h2,90} hf optg,saveact='allene.act',savexyz='allene.xyz' !default optimization !using model hessian. !Save optimized variables !in file allene.act !Save optimized geometry !in xyz style in !file allene.xyz

### Optimization using natural internal coordinates (BMAT)

- examples/allene_opt_bmat.inp
***, Allene geometry optimization using natural internal coordinates basis=sto-3g rcc=1.32 ang rch=1.08 ang acc=120 degree symmetry,nosym Geometry={C1; !Z-matrix input C2,c1,rcc Q1,c1,rcc,c2,45 C3,c2,rcc,c1,180,q1,0 h1,c1,rch,c2,acc,q1,0 h2,c1,rch,c2,acc,h1,180 h3,c3,rch,c2,acc,h1,90 h4,c3,rch,c2,acc,h2,90} hf; optg !default optimization using model hessian coord,bmat !use natural internal coordinates optg,coord=bmat !same as above

### MP2 optimization using a procedure

- examples/allene_optmp2.inp
***, Allene geometry optimization using Z-Matrix basis=vdz rcc=1.32 ang rch=1.08 ang acc=120 degree Geometry={C1 !Z-matrix input C2,c1,rcc Q1,c1,rcc,c2,45 C3,c2,rcc,c1,180,q1,0 h1,c1,rch,c2,acc,q1,0 h2,c1,rch,c2,acc,h1,180 h3,c3,rch,c2,acc,h1,90 h4,c3,rch,c2,acc,h2,90} optg,procedure=runmp2 !use procedure optmp2 runmp2={hf;mp2} !procedure definition

### Optimization using geometry DIIS

- examples/caffeine_opt_diis.inp
***, CAFFEINE cartesian coordinates (XYZ format) basis=sto-3g geomtyp=xyz geometry={ 24 CAFFEINE CARTESIAN COORDINATES C 0.8423320060 -0.3654865620 0.0000000000 C -0.2841017540 -1.1961236000 0.0000000000 N 2.0294818880 -1.1042264700 0.0000000000 N 0.0774743850 -2.5357317920 0.0000000000 N -1.6472646000 -0.6177952290 0.0000000000 C 1.4531962870 -2.3678913120 0.0000000000 C 0.6373131870 1.1735112670 0.0000000000 C -1.7812691930 0.7688916330 0.0000000000 N -0.6771444680 1.6306355000 0.0000000000 O 1.6106752160 1.9349693060 0.0000000000 O -2.9202890400 1.2510058880 0.0000000000 C -0.9202462430 3.1094501020 0.0000000000 C -2.8623938560 -1.4824503660 0.0000000000 C 3.4552156930 -0.6811094280 0.0000000000 H 2.0878150460 -3.2451913360 0.0000000000 H -1.4989252090 3.4222116470 -0.8897886280 H -1.4989252090 3.4222116470 0.8897886280 H 0.0071905670 3.7148499490 0.0000000000 H -3.4903070930 -1.2888938190 -0.8907763360 H -3.4903070930 -1.2888938190 0.8907763360 H -2.6289534570 -2.5638654230 0.0000000000 H 4.1360211370 -1.5529079440 0.0000000000 H 3.6817059520 -0.0685850980 0.8931597470 H 3.6817059520 -0.0685850980 -0.8931597470 } hf optg,savexyz=caffeine.xyz !save optimized geometry in file caffeine.xyz coord,bmat !Optimization in natural internal coordinates method,diis !Optimization method: Geometry DIIS optg,coord=bmat,method=diis,savexyz=caffeine.xyz !same as above

### Transition state of the HCN – HNC isomerization

The first example shows how to do a MP2 transition state optimization. The initial Hessian is taken from a previous HF frequency calculation.

- examples/hcn_mp2_ts.inp
***, HCN <--> NHC Isomerization - Transition State Optimization and Frequencies l1=1.18268242 ang l2=1.40745082 ang a1=55.05153416 degree basis=3-21G symmetry,nosym geometry={ C N,1,l1 H,2,l2,1,a1} hf ! HF-SCF frequencies,analytical ! Vibrational frequencies for HF-SCF (analytical Hessian) mp2 ! MP2 optg,root=2,method=rf,readhess ! Transition State Search using Rational Function ! Optimizer and HF hessian frequencies ! Vibrational frequencies for MP2 (numerical Hessian) ---

The second example shows how to do a CCSD(T) optimization with an MP2 hessian. Note that currently the CCSD(T) gradient is computed numerically using finite energy differences, and this can take long time for larger molecules. The calculation of the MP2 hessian finite differences of analytical gradients.

- examples/hcn_ccsd_ts.inp
***, HCN <--> NHC Transition State Optimization and Frequencies rcn=1.18 ang rnh=1.40 ang alpha=55 degree basis=vtz geometry={ C N,1,rcn H,2,rnh,1,alpha} hf ccsd(t) optg,numerical,root=2,hessproc=runmp2 !Transition state optimization for ccsd(t) !using mp2 hessian frequencies !CCSD(T) frequencies (using numerical second derivatives) runmp2={hf;mp2} !procedure definition ---

The last example shows how to do a MRCI+Q (MRCI with Davidson correction) optimization with an CASPT2 hessian. As for CCSD(T), the MRCI+Q gradient as computed numerically, while the CASPT2 hessian is obtained using finite differences of analytical CASPT2 gradients.

- examples/hcn_mrci_ts.inp
***, HCN <-> NHC Isomerization - Transition State Optimization and Frequencies print,orbitals,civector rcn=1.18 ang rnh=1.40 ang alpha=55 degree basis=vtz geometry={ C N,1,rcn H,2,rnh,1,alpha} closed,4 ! global setting for casscf inactive space hf ! HF-SCF multi ! CASSCF mrci ! MRCI optg,root=2,variable=energd,hessproc=runrs2 !optimize mrci+q transition state !and caspt2 for hessian runrs2={multi;rs2} !procedure definition for caspt2 ---

### Reaction path of the HCN – HNC isomerization

The following input first optimizes the transition state, and then performs reaction path calculations in both directions. The results are plotted.

- examples/hcn_isomerization.inp
***, HCN <---> NHC Isomerization Reaction Path basis=3-21G rcn=1.18282 ang ! Starting geometry is transition state rnh=1.40745 ang alpha=55.05 degree symmetry,x ! Cs Symmetry geometry={ C N,1,rcn H,2,rnh,1,alpha} int rhf,gap_min=0 !set default shift to zero to avoid convergence to wrong state optg,root=2,saveact=hcn_ts,rewind ! Find and store the TS {optg,method=qsdpath,dir=1, numhess=5,hesscentral,saveact=hcn_path} ! find IRC in positive direction readvar,hcn_ts.act ! Reset geometry to TS {optg,method=qsdpath,dir=-1,numhess=5,hesscentral,saveact=hcn_path,append} !find IRC in negative direction readvar,hcn_path.act alpha=alpha*pi/180 !convert angle to radian table,irc,rcn,rnh,alpha,e_opt !tabulate results {table,irc,e_opt !plot energy profile as function of irc plot,file='hcn_eopt.plot'} {table,irc,rcn,rnh,alpha !plot distances and angle as function of irc plot,file='hcn_dist.plot'}

This produces the plots

### Optimizing counterpoise corrected energies

Geometry optimization of counterpoise corrected energies is possible by performing for the total system as well as for each individual fragment separate `FORCE`

calculations. The gradients and energies are added using the `ADD`

directive. This requires that `NOORIENT`

has been specified in the geometry input, in order to avoid errors due to unintended rotation of the system. This default can be disabled using the `NOCHECK`

option, see `ADD`

above.

The way a counterpoise corrected geometry optimization works is shown in the following example. Note that the total counterpoise corrected energy must be optimized, not just the interaction energy, since the interaction energy depends on the monomer geometries and has a different minimum than the total energy. The interaction energy could be optimized, however, if the monomer geometries were frozen. In any case, the last calculation before calling `OPTG`

must be the calculation of the total system at the current geometry (in the example below the dimer calculation), since otherwise the optimizer gets confused.

- examples/hfdimer_cpcopt1.inp
***,HF dimer MP2/CP optimization with relaxed monomers basis=avtz gthresh,energy=1.d-8 ! INITIAL VALUES OF GEOMETRY VARIABLES RFF= 5.3 R1= 1.76 R2 = 1.75 THETA1 = 7.0 THETA2 = 111 symmetry,x orient,noorient geometry={ f1 f2 f1 rff h1 f1 r1 f2 theta1 h2 f2 r2 f1 theta2 h1 180.} label: text, CALCULATION AT LARGE SEPARATION rff_save=rff !save current rff distance rff=1000 !dimer calculation at large separation text, HF1 dummy,f2,h2; !second hf is now dummy {hf;accu,16} !scf for first monomer mp2; !mp2 for first monomer ehf1inf=energy !save mp2 energy in variable forces; !compute mp2 gradient for first monomer text, HF2 dummy,f1,h1; !first hf is now dummy {hf;accu,16} !scf for second monomer mp2; !mp2 for second monomer ehf2inf=energy !save mp2 energy in variable forces; !compute mp2 gradient for second monomer add,1 !add to previous gradient einf=ehf1inf+ehf2inf !total energy of unrelaxed momomers rff=rff_save !reset HF - HF distance to current value text, CP calculation for HF1 MONOMER dummy,f2,h2; !second hf is now dummy {hf;accu,16} !scf for first monomer mp2; !mp2 for first monomer ehf1=energy !save mp2 energy in variable forces; !compute mp2 gradient for first monomer add,-1 !subtract from previous gradient text, CP calculation for HF2 MONOMER dummy,f1,h1; !first hf is now dummy {hf;accu,16} !scf for second monomer mp2; !mp2 for second monomer ehf2=energy !save mp2 energy in variable forces; !compute mp2 gradient for first monomer add,-1 !subtract from previous gradient text, DIMER CALCULATION dummy !reset dummies {hf;accu,16} !scf for dimer mp2; !mp2 for dimer edimer=energy !save mp2 energy in variable forces; !compute mp2 gradient for dimer add,1 !add to previous gradient optg,gradient=1.d-4,startcmd=label: !find next energy text, compute optimized monomer energy rhf=r1 geometry={h1 F1,H1,rhf} {hf;accu,16} !scf for relaxed monomer mp2; !mp2 for relaxed monomer ehf=energy !save mp2 energy in variable optg !optimize HF structure text,optimized geometry parameters show,r1,r2,rhf,rff,theta1,theta2 text,computed interaction energies decpc=(einf-ehf1-ehf2)*tocm !counter poise correction de=(edimer-ehf1-ehf2)*tocm !CPC corrected interaction energy relative to unrelaxed monomers erelax=(2*ehf-einf)*tocm !relaxation energy derel=de-erelax !CPC corrected interaction energy relative to relaxed monomoers

The next example shows how the same calculations can be done using numerical gradients. In this case, first the total counter-poise corrected energy is formed and then optimized. Note that the `ADD`

command does not work for numerical gradients.

- examples/hfdimer_cpcopt1_num.inp
***,HF dimer MP2/CP optimization with relaxed monomers basis=avtz gthresh,energy=1.d-8 ! INITIAL VALUES OF GEOMETRY VARIABLES RFF= 5.3 R1= 1.76 R2 = 1.75 THETA1 = 7.0 THETA2 = 111 symmetry,x orient,noorient geometry={ f1 f2 f1 rff h1 f1 r1 f2 theta1 h2 f2 r2 f1 theta2 h1 180.} label: text, CALCULATION AT LARGE SEPARATION rff_save=rff !save current rff distance rff=1000 !dimer calculation at large separation text, HF1 dummy,f2,h2; !second hf is now dummy {hf;accu,16} !scf for first monomer mp2; !mp2 for first monomer ehf1inf=energy !save mp2 energy in variable text, HF2 dummy,f1,h1; !first hf is now dummy {hf;accu,16} !scf for second monomer mp2; !mp2 for second monomer ehf2inf=energy !save mp2 energy in variable einf=ehf1inf+ehf2inf !total energy of unrelaxed momomers rff=rff_save !reset HF - HF distance to current value text, CP calculation for HF1 MONOMER dummy,f2,h2; !second hf is now dummy {hf;accu,16} !scf for first monomer mp2; !mp2 for first monomer ehf1=energy !save mp2 energy in variable text, CP calculation for HF2 MONOMER dummy,f1,h1; !first hf is now dummy {hf;accu,16} !scf for second monomer mp2; !mp2 for second monomer ehf2=energy !save mp2 energy in variable text, DIMER CALCULATION dummy !reset dummies {hf;accu,16} !scf for dimer mp2; !mp2 for dimer edimer=energy !save mp2 energy in variable etot=edimer-ehf2-ehf1+ehf1inf+ehf2inf !total BSSE corrected energy optg,numerical,variable=etot,gradient=1.d-4,startcmd=label: !optimize geometry text, compute optimized monomer energy rhf=r1 geometry={h1 F1,H1,rhf} {hf;accu,16} !scf for relaxed monomer mp2; !mp2 for relaxed monomer ehf=energy !save mp2 energy in variable optg !optimize HF structure text,optimized geometry parameters show,r1,r2,rhf,rff,theta1,theta2 text,computed interaction energies decpc=(einf-ehf1-ehf2)*tocm !counter poise correction de=(edimer-ehf1-ehf2)*tocm !CPC corrected interaction energy relative to unrelaxed monomers erelax=(2*ehf-einf)*tocm !relaxation energy derel=de-erelax !CPC corrected interaction energy relative to relaxed monomoers

In the last example the monomer structures are kept fixed, and the interaction energy is optimized.

- examples/hfdimer_cpcopt2.inp
***,HF dimer MP2/CP optimization without monomer relaxation basis=avtz gthresh,energy=1.d-8 ! INITIAL VALUES OF GEOMETRY VARIABLES RFF= 5.3 THETA1 = 7 THETA2 = 111 symmetry,x orient,noorient geometry={ f1 f2 f1 rff h1 f1 1.74764059 f2 theta1 h2 f2 1.74764059 f1 theta2 h1 180.} !using fixed HF distances of isolated HF label: text, CP calculation for HF1 MONOMER dummy,f2,h2; !second hf is now dummy {hf;accu,16} !scf for first monomer mp2; !mp2 for first monomer ehf1=energy !save mp2 energy in variable forces; !compute mp2 gradient for first monomer scale,-1 !multiply gradient by -1 text, CP calculation for HF2 MONOMER dummy,f1,h1; !first hf is now dummy {hf;accu,16} !scf for second monomer mp2; !mp2 for second monomer ehf2=energy !save mp2 energy in variable forces; !compute mp2 gradient for first monomer add,-1 !subtract from previous gradient text, DIMER CALCULATION dummy !reset dummies {hf;accu,16} !scf for dimer mp2; !mp2 for dimer edimer=energy !save mp2 energy in variable forces; !compute mp2 gradient for dimer add,1 !add to previous gradient optg,gradient=1.d-5,startcmd=label: !find next energy text,optimized geometry parameters show,rhf,rff,theta1,theta2 text,computed interaction energies de=(edimer-ehf1-ehf2)*tocm !CPC corrected interaction energy with fixed monomers

### Adding CABS singles correction in CCSD(T)-F12 geometry optimization

Geometry optimization using CCSD(T)-F12 method including CABS singles corrections is possible by combining numerical gradients for HF+CABS singles method and analytical gradients for DF-CCSD(T)-F12 method obtained in separate `FORCE`

calculations. The gradients are added using the `ADD`

directive. Similarly to the example optimizing counterpoise corrected energies, `NOORIENT`

is required. The following example shows geometry optimization for the hydrogen fluoride molecule using DF-CCSD(T)-F12b method without and with CABS singles correction.

- examples/dfccsdtf12_3stc_tight_opt.inp
***, geometry optimization (tight thresholds) for HF using DF-CCSD(T)-F12b/3*C(FIXC,HY1) without and with CABS singles gthresh,energy=1.d-10,gradient=1.d-08 nosym noorient angstrom geometry={ 2 F 0.0000000000 0.0000000000 -0.0463705158 H 0.0000000000 0.0000000000 0.8745302823 } basis=avdz ! perform geometry optimization using CCCSD(T)-F12b/3*C(FIXC,HY1) without CABS singles by using analytical derivatives {df-hf;accu,18;} {df-ccsd(t)-f12,ansatz=3*C(fixc,hy1),cabs_singles=0,cabs=0,ri_basis=optri,thrcabs=1.d-10,thrcabs_rel=1.d-10,thrabs=1.d-10,thrabs_rel=1.d-10,gem_beta=1.0d0;cphf,thrmin=1.d-9;} {optg,energy=1.d-08,gradient=1.d-06;} ! save geometry parameters for printing struct,distvar=ccf12distwocabs,angvar=ccf12angwocabs,dihedvar=ccf12dihedwocabs ! perform geometry optimization using CCCSD(T)-F12b/3*C(FIXC,HY1) with CABS singles by using numerical and analytical derivatives label1 ! compute CABS singles forces by using numerical derivatives {df-hf;accu,18} {df-mp2-f12,cabs_singles=-1,ri_basis=optri,thrcabs=1.d-10,thrcabs_rel=1.d-10,thrabs=1.d-10,thrabs_rel=1.d-10,gem_beta=1.0d0} {forces,variable=ef12_singles,varsav=on,fourpoint} ! compute CCCSD(T)-F12b/3*C(FIXC,HY1) without CABS singles forces by using analytical derivatives {df-hf;accu,18;} {df-ccsd(t)-f12,ansatz=3*C(fixc,hy1),cabs_singles=0,cabs=0,ri_basis=optri,thrcabs=1.d-10,thrcabs_rel=1.d-10,thrabs=1.d-10,thrabs_rel=1.d-10,gem_beta=1.0d0;cphf,thrmin=1.d-9;} {forces,varsav=on;add,1} ! add the forces {optg,startcmd=label1,energy=1.d-08,gradient=1.d-06;} ! save geometry parameters for printing struct,distvar=ccf12dist,distdef=distance,angvar=ccf12ang,angdef=angle,dihedvar=ccf12dihed,diheddef=dihedral ! compute CABS singles correction for bond length distcorr=ccf12dist-ccf12distwocabs ! print results table,distance,ccf12distwocabs,ccf12dist,distcorr