Table of Contents

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.

Options for optimization methods

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

The following options are not available in PQSOPT:

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

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:

IGEN:

ITYPE:

ICONS:

Cluster optimizations:

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

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 I J value (Angstrom)
* BEND I J K value (degrees)
* TORS I J K L value (degrees)
* OUTP I J K L value (degrees)
* LINC I J K L value (degrees)
* LINP I J K L value (degrees)

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.].

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.

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.

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.

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.

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.

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

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:

Miscellaneous options:

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:

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.

For transition state searches (invoked with the ROOT option, see section transition state (saddle point) optimization (ROOT)) key can be

For reaction path following the input key is

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:

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:

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

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

In transition state optimizations type may be

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:

Transition state (saddle point) optimization (ROOT)

ROOT,root [root]

Specifies the eigenvector of the hessian to be followed.

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

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

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

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:

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:

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:

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:

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:

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

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