Structure and Reactivity¶

To study structure and reactivity common tasks are geometry optimization, transition state search and a frequency run. A frequency run can be done to obtain a good initial Hessian for a geometry optimization, or to check whether a proper minimum or transition state has been found.

The basis for locating stationary points on the PES are the nuclear gradients. (The user may assume that the program does that automatically when needed.)

Nuclear energy gradients¶

If you want the program to calculate the gradient of the energy w. r. t. nuclear displacements you should add the Gradients keyword. If the SelectedAtoms key is present, only the gradient of the selected atoms will be evaluated, and the rest will be set to zero.

Details about the theory of the analytical gradients within the BAND program can be found in the work of Kadantsev et al. [19].

Gradients (block-type)

Key to control the geometry optimization:

Gradients
  {LatticeGradients [false | true]}
  {LatticeStepSize step}
  {UseRelativeLatticeStep  [true | false]}
End

LatticeGradients: (Default: false) This key handles whether the lattice gradients are evaluated numerically (false) or analytically (true).
latticeStepSize: (Default: 0.012) With this keyword the step size for the numerical differentiation can be changed. When UseRelativeLatticeStep is set to true the step will be relative to the size of the displacement. (Unit is Bohr)

The user rarely needs to specify the Gradients keyblock as it is automatically added for geometry optimizations and frequency runs.

This option honors the SelectedAtoms key, in which case only the forces will be calculated for the selected atoms.

Lattice gradients¶

The step size for numerical differentiation is controlled with Gradients%latticeStepSize, and Gradients%UseRelativeLatticeStep.

LatticeGradients (block-type)

Key to handle special options for the calculation of analytical and numerical lattice gradients.

LatticeGradients
   {Analytical [true | false]}
   {AnalyticalPulay [true | false]}
   {AnalyticalKinetic [true | false]}
   {AnalyticalXC [true | false]}
   {AnalyticalElectrostatic [true | false]}
End

Analytical: (Default: false) If set to true, the following options AnalyticalElectrostatic, AnalyticalKinetic, AnalyticalPulay, and AnalyticalXC will be set to true automatically.
AnalyticalElectrostatic: (Default: true) If true the electrostatic energy contribution to the lattice gradients is calculated analytically, else it is calculated numerically.
AnalyticalKinetic: (Default: true) If true the kinetic energy contribution to the lattice gradients is calculated analytically, else it is calculated numerically.
AnalyticalPulay: (Default: true) If true the Pulay contribution to the lattice gradients is calculated analytically, else it is calculated numerically.
AnalyticalXC: (Default: true) If true the exchange-correlation energy contribution to the lattice gradients is calculated analytically, else it is calculated numerically.

In case you want to calculate forces acting on lattice vectors without performing an geometry optimization, you should use LatticeGradients%Analytical true!

Geometry optimization (GeoOpt)¶

The geometry can be optimized by adding the GeoOpt keyblock. Within the block you can specify the maximum number of cycles and the convergence criterion. The optimizer works with Cartesian coordinates and a Quasi Newton stepper. The initial Hessian is by default the unit matrix, but can also be loaded from a previous geometry optimization or frequency run, also from faster methods as DFTB and UFF. Various constraints are possible such as bond distances and angles. The GDIIS convergence accelerator is available but is disabled by default. A geometry optimization can be restarted from a previous result file. The geometry can be optimized non-relativistically, with scalar relativistic ZORA, or with spin-orbit coupling.

Convergence¶

A geometry optimization is considered converged when all the following criteria are met (this information is also printed to the logfile):

the difference between the bond energy at the current geometry and at the previous geometry step is smaller than \(\text{TolE}\)
the maximum Cartesian nuclear gradient is smaller than \(\text{TolG}\)
the root mean square (RMS) of the Cartesian nuclear gradients is smaller than \(\frac{2}{3} \text{TolG}\)
the maximum Cartesian step is smaller than \(\text{TolR}\)
the root mean square (RMS) of the Cartesian steps is smaller than \(\frac{2}{3} \text{TolR}\)

Note: if the maximum and RMS gradients are 10 times smaller than the convergence criterion, then criteria 4 and 5 are ignored

The convergence criteria can be changed via the Converge key:

GEOMETRY
   Converge {E=TolE} {Grad=TolG} {Rad=TolR} {Angle=TolA}
END

TolE: The criterion for changes in the energy, in Hartree. Default: 1e-3.
TolG: The criterion for the nuclear gradients, in Hartree/Angstrom. Default: 1e-3.
TolR: Refers to changes in the Cartesian coordinates or bond lengths, depending on in what coordinates you optimize, in Angstrom. Default: 1e-2.
TolA: Refers to changes in bond- and dihedral angles, in degrees. This is only meaningful if optimization takes place in Z-matrix coordinates. Default: 0.5 degree.

If only a numerical value is supplied as argument for Converge, rather than a specification by name, it is considered to apply to the gradients (only):

GEOMETRY
   Converge 1E-4
End

is equivalent to:

GEOMETRY
   Converge grad=1E-4
End

Numerical frequencies (Hessian)¶

By specifying the RunType Frequencies you can calculate the Hessian which in turn yields the harmonic frequencies and vibrational modes (q → 0 limit of a phonons calculation). The Hessian is an important product that can be used in a geometry optimization, and in particular to start a transition state search. The second derivative of the energy (i.e. the Hessian) is obtained by numerically differentiating the analytical gradients. When your unit cell has N atoms in the unit cell and no symmetry (other than translational symmetry) the number of displacements is 2x3xN. For large systems this can be very time consuming, and you could consider to calculate only a partial Hessian. If you have a very symmetric unit cell, but are interested in asymmetric modes, set the subkeys useA1Displacements and doA1Projection to false. Note: the lattice vectors are assumed to be fixed.

A frequency run is invoked with

RunType
  Frequencies
End

The rest is controlled by the Frequencies key.

Frequencies (block-type)

Key to control the numerical frequency run:

Frequencies
  {Step size}
  {useA1Displacements [true | false]}
  {doA1Projection [true | false]}
  {StaticCosmoSurface [true | false]}
End

Step: (Default: 0.01) The step size (in Å) for the numerical differentiation.
useA1Displacements: (Default: true) Determine only the total symmetric modes. When disabled, NOSYM calculations will be performed.
doA1Projection: (Default: true) Symmetrize the Hessian. For most users the only sensible setting is to make it equal to the value of useA1Displacements.
StaticCosmoSurface: (Default: true) Keep Cosmo surface (if any) constant. Although this is an approximation it is numerically much more stable. If you set it to false there is the risk that the number of Cosmo points is not constant for all finite displacements.

This option honors the SelectedAtoms key, in which case only the Hessian will be calculated for the selected atoms.

Transition state search¶

A transition state search is technically a geometry optimization. It is a search in the 3xN space for a point with vanishing gradients. The only difference is that it will try to go uphill in the direction of the lowest vibrational mode. The normal procedure is to guess a point that is as close as possible to the transition state. In that point you do a frequency run, which should hopefully produce a lowest vibrational mode in the direction of the transition state. You then start the Transition state search from that point, using the Hessian of the frequency run. At the end you can do again a frequency run to check that there is indeed only one single negative vibrational mode.

The transition state search is invoked like this

RunType
   TS
End

The rest is controlled by the GeoOpt key as in a normal geometry optimization.

Partial Hessian and (pre-)optimizations¶

If you consider the interaction of a molecule with a surface it may be initially convenient to freeze some of the lower layers in the slab, making calculations much cheaper. This can be done for the evaluation of the Hessian (see SelectedAtoms) and during a geometry optimization or transition state search (see Constraints).

Caution

When using a partial Hessian to start a geometry optimization note that the atoms that you select during the frequency run, are the ones that should not be frozen in the geometry optimization. In a way the selection needs to be inverted.

Constrained optimization¶

During a geometry optimization certain constraints can be enforced with the Constraints key. The constraints refer to the coordinates of the unit cell as specified on input.

Constraints (block-type)

Key to control constraints during a geometry optimization:

Constraints
  {Atom aI {x y z}
  {Coord aI coordI {coordIvalue}
  {Dist aI aJ distance}
  {Angle aI aJ aK angle}
  {Dihed aI aJ aK aL angle}
End

Atom: Freeze the position of atom a1 to the current coordinate, or if specified to (x y z). This constraint will make the calculation cheaper, because the gradients of frozen atoms need not be calculated.
Dist: Constrain the distance between atoms a1 and a2 to distance.
Angle: Constrain the angle between atoms a1, a2, and a3 to angle.
Dihed: Constrain the dihedral angle between atoms a1, a2, a3, and a4 to angle.

Selected atoms¶

When doing a frequency or geometry optimization run you can restrict it to a limited number of selected nuclei. The selection may not break the symmetry

SelectedAtoms a1 a2 ... an

The list of atoms with indices ‘a1 a2 ... an’ are considered selected.

Note: Currently this option only has effect in a numerical frequency run and with the Gradients keyword.

Run Types (GeoOpt, TS, Frequencies and Phonons)¶

RunType (block-type)

Key to invoke the following calculation types:

Geometry Optimization

RunType
  GeoOpt
End

Transition State Search

RunType
  TS
End

Frequency Calculation

RunType
  Frequencies
End

Phonon Calculation

RunType
  Phonons
End