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 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
  {oldio [0 | 1]}
  {LatticeGradients [false | true]}
  {LatticeStepSize step}
  {UseRelativeLatticeStep  [true | false]}
End

oldio: (Default: 1) There are two slightly different implementations. The oldio=1 implementation has a better scaling behavior for large systems, although it can be a bit slower for small ones. When used with MetaGGAs oldio=1 is required.
LatticeGradeints: (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 this key as it is automatically added for geometry optimizations and frequency runs.

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

Lattice gradients¶

The calculation of the forces, acting on the lattice vectors, is invoked by performing a geometry optimization, including lattice vectors. (see the GeoOpt%OptimizeLattice)

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
   {AnalyticalPulay [true | false]}
   {AnalyticalKinetic [true | false]}
   {AnalyticalXC [true | false]}
   {AnalyticalElectrostatic [true | false]}
End

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.

Geometry optimization (GeoOpt)¶

The geometry can be optimized by adding the GeoOpt key block. 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. 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.

GeoOpt (block-type)

Key to control the automatic geometry optimization:

GeoOpt
   {OptimizeLattice [false | true]}
   {Iterations n}
   {Converge {Grad=tolgrad} {E=tolE} {Step=tolStep}}
   {TrustRadius radius}
   {UseVariableTrustRadius var}
   {InitialHessian [UnitMatrix | filename | Swart]}
   {RestartSCF [true | false]}
   {RestartFit [true | false]}
   {GDIIS [true | false]}
   {StaticCosmoSurface [true | false]}
   {TSMode mode}
End

OptimizeLattice: (Default: false) Whether or not to optimize the lattice vectors. The lattice gradients are obtained by numerical differentiation, see GRADIENTS%LatticeStepSize.
Iterations: (Default: 50) maximum number of cycles.
Converge {Grad=tolgrad} {E=tolE} {Step=tolStep}: various criteria to determine the convergence can be specified on line. Here tolGrad is the maximum gradient allowed (Default: 1.0e-2 Hartree per Angstrom), tolE is the maximum energy change allowed (Default: 1.0e-3 Hartree), tolStep is the maximum step size (Default: 3e-2 Bohr)
TrustRadius: (Default: 0.2 Bohr) The step sizes taken by the optimizer will be limited to this value. If the proposed step is larger it will be scaled back.
UseVariableTrustRadius: (Default: false) Automatic adjustment of the trust radius based on the observed energy changes.
InitialHessian: (Default: UnitMatrix) The initial Hessian is read from a “RUNKF” file filename of a previous frequency calculation. Another option, Swart, is to use the empirical force field derived Swart Hessian.
RestartSCF: (Default: true) During the optimization try to restart the SCF from the previous geometry step.
RestartFit: (Default: false) During the optimization try to restart the SCF from fit coefficients of the previous geometry step. Can not be used simultaneously with RestartSCF.
GDIIS: (Default: false) Use the GDIIS convergence accelerator.
StaticCosmoSurface: (Default: false) Keep the Cosmo surface (if any) constant.
TSMode: Hessian normal mode index to follow during transition state search. The specified mode is selected on the first iteration and is followed regardless of its index on subsequent iterations. Modes with zero frequencies are not counted.

Note that a (failed) geometry optimization can be continued with the Restart key by specifying as option GeometryOptimization.

Numerical frequencies (Hessian)¶

By specifying the RunType Frequencies you can calculate the Hessian which in turn yields the harmonic frequencies and vibrational modes. 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 options 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 a 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). When using a partial Hessian to start a geometry optimization note this. Note: 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.

Phonons and thermodynamics¶

Atoms vibrate about their equilibrium positions, giving rise to lattice vibrations, called phonons. BAND can calculate phonon dispersion curves within standard harmonic theory, implemented with a finite difference method in the spirit of the program PHON [45]. Within the harmonic approximation we can calculate the partition function and from that thermodynamic properties, such as the specific heat and the free energy.

Our implementation can also handle 1D and 2D periodic systems.

To run a phonon calculation you need to do the following steps (see also: Phonons example):

Optimize the structure, including the lattice vectors. The resulting geometry is the starting geometry for the phonon run.
Specify a super cell transformation. In principle this should be as large as possible. In practice one may want to start with a 2x2x2 cell. Set the RunType to Phonons:
```
RunType
   Phonons
End
```
Examine with the GUI the dispersion curves, and check for negative frequencies. If these exist then either:
- the geometry was not close enough to the minimum. Try specifying a smaller gradient convergence criteria in your geometry optimization
- the super cell transformation was too small
- the numerical precision of the phonon run was not good enough.

PhononConfig (block-type)

Key to control the phonon run:

PhononConfig
  {SuperCell (Sub block type)}
  {StepSize step}
  {nSides [1|2]}
  {MinTempKelvin tmin}
  {MaxTempKelvin tmax}
  {NumTemperatures ntemp}
End

SuperCell: This sub block key should hold the supercell transformation, which expresses the lattice vectors in terms of the vectors of the primitive cell.
StepSize: (Default: 0.076) The step size taken to obtain all force constants.
nSides: (Default: 2) By default a two-sided (or quadratic) numerical differentiation of the nuclear gradients is used. Using a single-sided (or linear) numerical differentiation is computationally faster but much less accurate (Note: in older versions of the program only the single-sided option was available).
MinTempKelvin: (Default: 0.0 Kelvin) Minimum temperature for thermodynamic plots.
MaxTempKelvin: (Default: 1000.0 Kelvin) Maximum temperature for thermodynamic plots.
NumTemperatures: (Default: 1000) Number of temperatures for thermodynamic plots.

Here is an example input:

PhononConfig
  StepSize 0.0913
  SuperCell
    2 0 0
    0 2 0
    0 0 2
  SubEnd
End

In the standard output the thermodynamic functions are printed:

========================
ThermoDynamic Properties
========================
Zero-point Energy (Hartree)            0.00448537
Zero-point Energy (eV)                 0.12205313

==================================================================================================================
Temperature(K) Internal energy(Hartree)              Entropy(kB)     Free Energy(Hartree)        Specific Heat(kB)

  0.000000E+00             0.448537E-02             0.000000E+00             0.448537E-02             0.000000E+00
  0.100000E+01             0.448540E-02             0.116921E-01             0.448536E-02             0.329896E-01
  0.200000E+01             0.448566E-02             0.643295E-01             0.448525E-02             0.138532E+00
  ...                      ...                      ...                      ...                      ...