Next: 4.25 vibrot Up: 4. Programs Previous: 4.23 seward

4.24 `slapaf`

Provided with the first order derivative with respect to the nuclei the program is capable to optimize molecular structures with or without constraints for minima or transition states. This will be achieved with a quasi-Newton approach in combination with 2nd ranks updates of the approximate Hessian or with the use of an analytic Hessian.

4.24.1 Description

SLAPAF has three different ways in selecting the basis for the displacements during the optimization. First the old format (as in MOLCAS-3) which is user specified. The internal coordinates are here represented as linear combination of internal coordinates (such as bonds, angles, torsions, out of plane angles, cartesian coordinates and separation of centers of mass) and the linear combinations are totally defined by user input. This format does also require the user to specify the Hessian (default the unit matrix). This option does also allow for frozen internal coordinates. The two other options are totally black boxed. The first automatic option employs the Cartesian eigenvectors of the approximative Hessian (generated by the Hessian model functional [64]). The second automatic option utilizes linear compinations of some curvilinear coordinates (stretches, bends, and torsions). This implementation has two variations. The first can be viewed as the conventional use of non-redundant internal coordiantes [77,86,12]. The second variation is a force constant weighted (FCW) redundant space (the HWRS option) version of the former implementation [65]. All versions of internal coordinates can be used in combinations with constraints on the molecular parameters which do not have to fulfill the constraint at the start of the optimization.

The displacements are symmetry adapted and any rotation and translation if present is deleted. The relaxation is symmetry preserving.

4.24.2 Dependencies

SLAPAF depends on the results of ALASKA and also possibly on MCKINELY and MCLR.

4.24.3 Files

4.24.3.1 Input files

Apart from the standard input file SLAPAF will use the following input files.

File Contents
ONEINT Retrieve symmetry information
RELAX Current cartesian gradients are read from this file. In addition to this SLAPAF will retrieve the Hessian matrix, previous internal coordinates, gradients and shifts.
COMFILE Auxiliary data.

4.24.3.2 Output files

In addition to the standard output file SLAPAF will use the following output files.

File Contents
COMFILE Auxiliary data.
RELAX Updated Hessian matrix, shift vectors and gradient vectors in the internal coordinate basis.
MOLDEN Molden input file for geometry optimization analysis.

4.24.4 Input

4.24.5 `SLAPAF`

SLAPAF will as standard provided with an energy and a corresponding gradient update the geometry (optimize). Possible update methods include different quasi-Newton methods. The program will also provide for updates of the Hessian. The program has a number of different variable metric methods available for the Hessian update. This section describes the input to the SLAPAF program. The input is always preceded by the dummy namelist input

&SLAPAF &END

Compulsory keywords

Keyword Meaning
END of Input This marks the end of the input to the program.

Optional convergency controll keywords

Keyword Meaning
ITERations On the next lines follows the max number of iterations which will be allowed in the relaxation procedure. Default is 50 iterations.
THRShld This keyword is followed by two real numbers on the next line which specifies the convergency criterion with respect to the energy change and the norm of the gradient. The defaults are 0.5D-7 and 0.5D-5.
BAKEr Activate convergency criterions according to Baker [11]. Default is to use the convergency criterions as in Gaussian program [41].
MAXStep This keyword is followed by the value which defines the largest change of the internal coordinates which will be accepted. A change which is larger is reduced to the max value. The default value is 0.5 au or rad.
NOMAxstep This keyword indicats that there should be made no modifications to the value of large changes of the internal coordinates. The default is to reduce large changes.

Optional coordinate selection keywords

Keyword Meaning
CARTesian coordinates This keyword will make SLAPAF use the eigenvectors of the approximative Hessian expressed in cartesians as the definition of the internal coordinates. The default is to use the FCW non-reduntant internal coordinates. The Hessian will be modelled by the Hessian Model Functional.
CONStraints This marks the start of the definition of the constraints which the optimization is subject to. This section is always ended by the keyword End of Constraints. For a complete description of this keyword see the section 4.24.5. This option can be used in conjunction with any definition of the internal coordinates. This option will automatically turn of the line search, and change the Hessian update method (if used) from BFGS to MSP. The defaults is to apply no constraints to the optimization.
INTErnal coordinates This marks the start of the definition of the internal coordinates. This section is always ended by the keyword End of Internal. For a complete description of this keyword see the section 4.24.5. This option will also as default use a unit matrix as default for the Hessian matrix. The default is to use the FCW non-reduntant internal coordinates.
HWRS Use the force constant weighted (FCW) redundant space version of the nonredundant internal coodinates. This is the default. The Hessian will be modelled by the Hessian Model Functional.
NOHWrs Disable the use of the force constant weighted redundant space version of the nonredundant internal coodinates. The default is to use the HWRS option. The Hessian will be modelled by the Hessian Model Functional.
RTHR Change the thresholds for including redundant coordinates. Followed by one line with three real entries, corresponding to bonds, bends, and torsions. Default values are 0.2, 0.2, and 0.2.

Optional Hessian update keywords

Keyword Meaning
NOUPdate No update is applied to the Hessian matrix. Default is that the Broyden-Fletcher-Goldfarb-Shanno update is applied.
MEYEr Activate update of the Hessian matrix according to Meyer's method [83,39]. This method does not allow for any modifications of the suggested change of the geometry as suggested by the Hessian and the gradient. Default is that the Broyden-Fletcher-Goldfarb-Shanno update [21,38,46,112] is applied.
BPUPdate Activate update accoding to Broyden-Powell [28]. Default is that the Broyden-Fletcher-Goldfarb-Shanno update is applied. The Broyden-Powell update is recommended in searches for transition states.
BFGS Activate update according to Broyden-Fletcher-Goldfarb-Shanno. This is the default.
MSP-update Activate the Murtagh-Sargent-Powell update according to Bofill [20]. This update is prefered for the localization of transition states.

Optional optimization procedure keywords

Keyword Meaning
RATIonal Functional Activate geometry optimization using the Rational Functional approach [13], this is the default.
C1-Diis Activate geometry optimization using the c1-gdiis method [26,84,85]. The default is to use the Rational Functional approach.
C2-Diis Activate geometry optimization using the c2-gdiis method [109]. The default is to use the Rational Functional approach.
DXDX This option is associated to the use of the c1- and c2-gdiis procedures. This option will activate the computation of the so-called error matrix elements as $e=\delta x^{\dagger}\delta x$ , where $\delta x$ is the displacement vector.
DXG This option is associated to the use of the c1- and c2-gdiis procedures. This option will activate the computation of the so-called error matrix elements as $e=\delta x^{\dagger}g$ , where $\delta x$ is the displacement vector and is the gradient vector.
GDX See above.
GG This option is associated to the use of the c1- and c2-gdiis procedures. This option will activate the computation of the so-called error matrix elements as $e=g^{\dagger}g$ , where is the gradient vector.
NEWTon Activate geometry optimization using the standard quasi-Newton approach. The default is to use the Rational Functional approach.
NOLIne search Disable line search. Default is to use line search for minima.
TS Keyword for optimization of transition states. This flag will activate the use of the mode following rational functional approach [54]. The mode to follow can either be the one with the lowest eigenvalue (if positive it will be changed to a negative value) or by the eigenvector which index is specified by the MODE keyword (see below). The kewyword will also active the Murtagh-Sargent-Powell update of the Hessian and inactivate line search. This keyword will also enforce that the Hessian has the right index (i.e. one negative eigenvalue).
MODE Specification of the Hessian eigenvector index, this mode will be followed by the mode following RF method for optimization of transition states. The keyword card is followed by a single card specifying the eigenvector index.

Optional force constant keywords

Keyword Meaning
Schlegel The approximate Hessian is computed according to Schlegel [106]. The default is to compute the approximate Hessian with the Hessian model functional [64].
OLDForce constant matrix The Hessian matrix is read from the file RLXOLD.
FCONstant matrix Input of Hessian. There are two different syntaxes.

The keyword is followed by a line with the number of elements which will be set (observe that the update will preserve that the elements $H_{ij}$ and $H_{ji}$ are equal). The next lines will contain the value and the indices of the elements to be replaced.
The keyword if followed by the label Square or Triangular. The subsequent line specifies the rank of the Hessian. This is then followed by lines specifying the Hessian in square or lower triangular order.
NUMErical hessian This envokes as calculation of the force constant matrix by a two-point finite difference formula. The resulting force constant matrix is used for an analysis of the harmonic frequences. Observe that in case of the use of internal coordinates defined as cartesian coordinates that these has to be linear combinations which are free from translational and rotational components for the harmonic frequency analysis to be valid.
CUBIc force constants This envokes a calculation of the 2nd and the 3rd order force constant matrix by finite difference formula.
DELTa This keyword is followed by a real number which defines the step length used in the finite differentiation. Default: 1.0D-2.

Optional miscellaneous keywords

Keyword Meaning
EXTRact This keyword will make the program append a brief output to the file EXTRACT.
RTRN Max number of atoms for which bond lengths, angles and dihedral angles are listed, and the radius defining the maximum length of a bond follows on the next line. The latter is used as a threshold when printing out angles and dihedral angles. The length can be followed by Bohr or Angstrom which indicates the unit in which the length was specified, the default is Bohr. The default values are 15 and 3.0 au.

4.24.5.1 Definition of internal coordinates or constraints

The input section defining the internal coordinates always start with the keyword Internal coordinates and the definition of the constraints starts with the keyword Constraints.

The input is always sectioned into two parts where the first section defines a set of primitive internal coordinates and the second part defines the actual internal coordinates as any arbitrary linear combination of the primitive internal coordinates that was defined in the first section. In case of constraints the second part does also assign values to the constraints.

In the first section we will refer to the atoms by their atom label (SEWARD will make sure that there is no redundancy). In case of symmetry one will have to augment the atom label with a symmetry operation in parenthesis in order to specify a symmetry related center. Note that the user only have to specify distinct internal coordinates (ALASKA will make the symmetry adaptation). The symmetry adapted coordinates will be the sum of all symmetry related partners multiplied with one over the square root of the number of the degeneracy of the internal coordinate (for example, the symmetry adapted bond in water will be $\frac{1}{Sqrt{2}} ( R_1 + R_2 )$ ). This definition is of most importance in optimization with contraints when the value of a symmetry adapted constraint is computed.

In the specification below rLabel is a user defined label with no more than 8 (eight) characters. The specifications atom1, atom2, atom3, and atom4 are the unique atom labels as specified in the input to SEWARD.

The primitive internal coordinates are defined as

rLabel = Bond atom1 atom2 -- a primitive internal coordinate rLabel is defined as the bond between center atom1 and atom2.
rLabel = Angle atom1 atom2 atom3 -- a primitive internal coordinate rLabel is defined as the angle between the bonds formed from connecting atom1 to atom2 and connecting atom2 to atom3.
rLabel = LAngle(1) atom1 atom2 atom3 -- a primitive internal coordinate rLabel is defined as the linear angle between the bonds formed from connecting atom1 to atom2 and connecting atom2 to atom3. To define the direction of the angle the following procedure is followed.
1. - the three centers are linear,
  1. form a reference axis, R1, connecting atom1 and atom3,
  2. compute the number of zero elements, nR, in the reference vector,
    1. - nR=0, a first perpendicular direction to the reference axis is formed by
      
      $\begin{displaymath}{\bf R2}=(R1x,R1y,-R1z)\end{displaymath}$
      
      followed by the projection
      
      $\begin{displaymath}{\bf R2}={\bf R2}-\frac{ {\bf R2} \cdot {\bf R1} }{{\bf R1} \cdot {\bf R1}} {\bf R1}.\end{displaymath}$
      
      The second perpendicular direction completes the right-handed system.
    2. - nR=1, a first perpendicular direction to the reference axis is defined by setting the element in R2 corresponding to the zero entry in R1 to unity. The second perpendicular direction completes the right-handed system.
    3. - nR=2, a first perpendicular direction to the reference axis is defined by setting the element corresponding to the first zero entry in R1 to unity. The second perpendicular direction completes the right-handed system.
2. - the three centers are nonlinear, the first bending direction is the one which is in the plane formed by atoms atom1, atom2, and atom3. The second bending direction is taken as the direction perpendicular to the same plane.
The direction of the bend for LAngle(1) is taken in the direction of the first perpendicular direction.
rLabel = LAngle(2) atom1 atom2 atom3 -- a primitive internal coordinate rLabel is defined as the linear angle between the bonds formed from connecting atom1 to atom2 and connecting atom2 to atom3. The definition of the perpendicular directions is as described above. The direction of the bend for LAngle(2) is taken in the direction of the second perpendicular direction.
rLabel = Dihedral atom1 atom2 atom3 atom4 -- a primitive internal coordinate rLabel is defined as the angle between the planes formed of atom1, atom2 and atom3, and atom2, atom3 and atom4, respectively.
rLabel = OutOfP atom1 atom2 atom3 atom4 -- a primitive internal coordinate rLabel is defined as the angle between the plane formed by atom2, atom3, and atom4 and the bond formed by connecting atom1 and atom4. This option is not available for constraints.
rLabel = Dissoc (n1+n2) atom1 atom2 atom3 ... atomN -- a primitive internal coordinate rLabel is defined as the distance between the center of masses of two sets of centers. The first center has n1 members and the second has n2. The input contains the labels of the atoms of the first group followed immediately by the labels of the second group. This option is not available for constraints.
rLabel = Cartesian i atom1 -- a primitive internal coordinate rLabel is defined as the pure cartesian displacement of the center labelled atom1. The label i is selected to x, y, or z to give the appropriate component.

The second section starts with the label Vary or in the case of contstraints with the label Values. In this section the user specifies all symmetric internal coordinates excluding translation and rotation. In case of a defintion of internal coordinates the section contains either a direct reference to a rLabel being defined in the first section or will use expressions like

label = fac1 rLabel1 + fac2 rLabel2 + ....

which defines an internal coordinate as the linear combination of the primitive internal coordinates rLabel1, rLabel2, ... with the coefficients fac1, fac2, ..., respectively.

Observe, if some internal coordinates are chosen to be fixed they should be defined after the label Fix and the Vary section always preceeds the section in which the fixed coordinates are specified.

In case of a definition of constraints the sections contains either a direct reference to a rLabel as in

rLabel = rValue [Angstrom,Degrees]

where rValue is the desired value of the constraint in au or rad. Like in the case of defining internal coordinates one can also use expressions like

fac1 rLabel1 + fac2 rLabel2 + .... = Value [Angstrom,Degrees]

Example: A definition of user specified internal coordines of benzene. The molecule is in D $_{6h}$ and since MOLCAS-3 only uses up to D $_{2h}$ the Fix option is used to constrain the relaxation to the higher point group. Observe that this will only restrict the nuclear coordinates to D $_{6h}$ . The electronic wavefunction, however, can have lower symmetry.

Internal coordinates
r1 = Bond C1 C2
r2 = Bond C1 H1
r3 = Bond C2 H2
r4 = Bond C2 C2(x)
f1 = Angle H1 C1 C2
f2 = Angle H2 C2 C1
Vary
a = 1.0 r1 +  1.0 r4
b = 1.0 r2 +  1.0 r3
c = 1.0 f1 +  1.0 f2
Fix
a = 1.0 r1 + -1.0 r4
b = 1.0 r2 + -1.0 r3
c = 1.0 f1 + -1.0 f2
End of Internal

Example: A complete set of input decks for a CASSCF geometry optimization. These are the input decks for the optimization of the enediyne molecule.

 &SEWARD &END
Title
Enediyne MCSCF structure
Symmetry
x z

Basis set
C.ANO...5s4p2d.
C1             1.2869761127        2.0799281025         .0000000000
C2             2.8355091288        -.1380881195         .0000000000
C3             4.1954709187       -1.9656839604         .0000000000
End of basis

Basis set
H.ANO...3s2p.
H1             2.2478721352        3.8639049616         .0000000000
H2             5.3554366293       -3.5799988030         .0000000000
End of basis

End of input

 &SCF &END
Title
Enediyne

ITERATIONS
 30
Occupied
9 8 2 1
Thresholds
 1.d-8  .5d-8
IVO
End of input

 &RASSCF &END
Lumorb
Title
 Bergman reaction
NactEl
12 0 0
Spin
1
Inactive
7 7 0 0
Ras2
3 3 3 3
Iterations
50 50
CiRoot
1 1
1
Thrs
 1.0e-08 1.0e-05 1.0e-05
Symmetry
 1
End of input

 &ALASKA &END
End of input

 &SLAPAF &END
Iterations
20
End of input

Example: A input for the optimization of water constraining the structure to be linear at convergency.

 &Structure &End
method scf
End of Input 
  
 &SEWARD &END
Title
 H2O geom optim, using the ANO-S basis set.
Pkthre
1.0D-11 

Basis set
H.ANO-S...1s.
H1      1.43354233    .00000000    .95295406
H2     -1.43354233    .00000000    .95295406
End of basis

Basis set
O.ANO-S...2s1p.
O        .00000000    .00000000    .00000000
End of basis

End of input

 &SCF &END             
Title                 
 H2O 
ITERATIONS         
 40               
Occupied         
 5
End of input 

 &ALASKA &END
End of input

 &SLAPAF &END
Iterations
15
Constraints
a1 = langle(1) H1 O H2
Values    
a1 = 180.00 degrees
End of Constraints
End of input