# Gradients, Hessian, Stress tensor, Elasticity¶

No matter what application the AMS driver is used for, in one way or another it
always explores the potential energy surface (PES) of the system. One can
furthermore ask AMS to calculate additional properties of the PES in the points
that are visited. These are mostly (but not exclusively) derivatives of the energy, e.g. we can ask
AMS to calculate the gradients or the Hessian in the visited points. In general
all these PES point properties are requested through the `Properties`

block in
the AMS input:

```
Properties
Gradients Yes/No
StressTensor Yes/No
Hessian Yes/No
SelectedRegionForHessian string
PESPointCharacter Yes/No
ElasticTensor Yes/No
End
```

This properties described on this page in the AMS manual are all related to derivatives of the energy.

Note that because these properties are tied to a particular point on the
potential energy surface, they are found on the engine output
files. Note also that the properties are not always
calculated in every PES point that the AMS driver visits during a calculation.
By default they are only calculated in *special* PES points, where the
definition of special depends on the task AMS is performing: For
a geometry optimization properties would for
example only be calculated at the final, converged geometry. This behavior can
often be modified by keywords special to the particular running task.

## Nuclear gradients¶

The first derivative with respect to the nuclear coordinates can be requested as follows:

```
Properties
Gradients Yes/No
End
```

`Properties`

`Gradients`

- Type
Bool

- Default value
No

- GUI name
Nuclear gradients

- Description
Calculate the nuclear gradients.

Note that these are gradients, not forces, the difference being the sign. The
gradients are printed in the output and written to the engine result
file belonging to the particular point on the PES in the
`AMSResults%Gradients`

variable as a \(3 \times n_\mathrm{atoms}\) array
in atomic units (Hartree/Bohr).

## Hessian¶

The calculation of the second derivative of the total energy with respect to the nuclear coordinates is enabled by:

```
Properties
Hessian Yes/No
SelectedRegionForHessian string
End
```

`Properties`

`Hessian`

- Type
Bool

- Default value
No

- Description
Whether or not to calculate the Hessian.

`SelectedRegionForHessian`

- Type
String

- GUI name
Hessian only for

- Description
Compute the Hessian matrix elements only for the atoms in a particular region. If not specified, the Hessian will be computed for all atoms.

The Hessian is not printed to the text output, but is saved in the engine result
file as variable `AMSResults%Hessian`

. Note that this is just the plain
partial second derivatives (no mass-weighting) of the total energy. The
\(3 \times N_\mathrm{atoms}\) columns/rows of the matrix are grouped by
atom: the first three rows/columns correspond to the first atom, etc..

Note that the AMS driver also supports the Mobile Block Hessian (MBH) method, which allows treating parts of the system as rigid blocks.

Often one is not interested in the Hessian matrix itself, but just in using it for the calculation of IR frequencies or to characterize a PES point (as e.g. a local minimum or a saddle point). For these application, see the following pages in the manual:

## PES point character¶

A PES point can according to the slope and curvature of the PES at that point be classified in the following categories:

A local minimum on the PES with vanishing nuclear gradients and no negative frequencies.

A transition state with vanishing nuclear gradients and exactly one negative frequency, i.e. a first order saddle point on the PES.

A higher order saddle point, i.e. a stationary point on the PES with vanishing nuclear gradients but more than one imaginary frequency.

A non-stationary point on the PES. Here the gradients are non-zero.

This classification can easily be done if both the gradients and the normal
modes have already been calculated. However, calculating the full Hessian needed
for the entire set of normal modes is very expensive and undesirable if one only
wants to know the character of a PES point. The AMS driver can quickly, and
without calculating the full Hessian, characterize a PES point into one of the
above categories. This can be used to confirm the success of e.g. a
transition state search or geometry
optimization. A PES point can be characterized by
requesting `PESPointCharacter`

as a property:

```
Properties
PESPointCharacter Yes/No
End
```

`Properties`

`PESPointCharacter`

- Type
Bool

- Default value
No

- GUI name
Characterize PES point

- Description
Determine whether the sampled PES point is a minimum or saddle point. Note that for large systems this does not entail the calculation of the full Hessian and can therefore be used to quickly confirm the success of a geometry optimization or transition state search.

This will calculate the few lowest normal modes using an iterative diagonalization of the Hessian 1 based on a Davidson algorithm implemented in the PRIMME library 2. The procedure has been optimized for finding a small number of low-lying eigenvalues in as few matrix-vector multiplications (and thus single point calculations) as possible. This is facilitated by performing the iterative method using a pre-conditioning matrix based on an approximation of the Hessian. The approximate Hessian is obtained from the full Hessian at a lower level of theory. This calculation also provides the initial guesses for the desired normal modes. What the lower level of theory is depends on the main engine used in the calculation: DFTB with the GFN1-xTB model is used as the lower level of theory for relatively slow engines, e.g. DFT based engines. For semi-empirical engines like DFTB or MOPAC, the lower level of theory is currently UFF. It is currently not possible to change the engine used to obtain the preconditioning Hessian and the approximate modes.

Note that the iterative calculation of the normal modes is skipped when …

… the nuclear gradients are so large that the PES point is considered non-stationary. The calculation of the modes is then just not necessary for classifying it.

… the full normal modes or Hessian have also been requested. The iterative calculation is then not necessary, as all modes are already known.

… the molecule is very small. (For small systems the iterative calculation of the few lowest normal modes is not faster than the full calculation of all modes, so all modes are calculated instead.)

The classification as a stationary or non-stationary point uses the gradient convergence criterion from the geometry optimizer as the tolerance, see geometry optimization. This makes sure that the criterion for what is considered converged/stationary is always in sync between the optimizer and the PES point characterization.

For periodic systems the PES point characterization does not take the lattice degrees of freedom into account. A PES point where the nuclear gradients are small enough would for example be classified as a stationary point, even if the system is under stress.

Details of the iterative procedure can be configured in the
`PESPointCharacter`

block:

```
PESPointCharacter
Displacement float
NegativeEigenvalueTolerance float
NumberOfModes integer
Tolerance float
End
```

`PESPointCharacter`

- Type
Block

- Description
Options for the characterization of PES points.

`Displacement`

- Type
Float

- Default value
0.04

- Description
Controls the size of the displacements used for numerical differentiation: The displaced geometries are calculated by taking the original coordinates and adding the mass-weighted mode times the reduced mass of the mode times the value of this keyword.

`NegativeEigenvalueTolerance`

- Type
Float

- Default value
-0.0001

- Unit
Hartree/Bohr^2

- Description
The threshold in Hessian eigenvalue below which a mode is considered imaginary, i.e. indicating a transition state. This is a small negative number, as very small negative eigenvalues may be due to numerical noise on an essentially flat PES and do not indicate true transition states.

`NumberOfModes`

- Type
Integer

- Default value
2

- Description
The number of (lowest) eigenvalues that should be checked.

`Tolerance`

- Type
Float

- Default value
0.016

- Description
Convergence tolerance for residual in iterative Davidson diagonalization.

Note that the residual tolerance that can be achieved is limited by the numerical differentiation that is performed. The default values should apply in most cases, but if convergence becomes a problem one may choose to increase the tolerance or to increase the step size (slightly). Note that the default residual tolerance is lower than for the other mode selective methods. This is because PRIMME uses a different convergence criteria than mode tracking/refinement. The higher value used as a default will therefore not result in decreased levels of accuracy. The method will bail if the number of iterations exceeds the number of normal modes as at this point still achieving convergence becomes unlikely, in part due to the next point.

In order to avoid producing the known and irrelevant rigid modes, the method searches for normal modes orthogonal to six (or five) rigid modes. Imperfections due to the numerical differentiation may mean that the translational and rotational rigid modes are not exact eigenmodes of the Hessian that is constructed. As a result, some part of the lowest vibrational normal mode may lie in the span of the theoretical rigid modes and therefore be inaccessible to the Davidson method. This places a lower bound on the residual tolerance that can be achieved, which is directly related to the numerical differentiation accuracy. The take-away: do not set the tolerance too low, the default usually suffices.

Behind the scenes, the method actually computes a few more modes than requested. In the case of multiplicities, eigenvalues may still converge out of order. These additional eigenvalues essentially guarantee that the obtained modes are indeed the lowest ones.

- 1
P. Deglmann and F. Furche,

*Efficient characterization of stationary points on potential energy surfaces*, J. Chem. Phys. 117, 9535 (2002)- 2
A. Stathopoulos and J. R. McCombs,

*PRIMME: PReconditioned Iterative MultiMethod Eigensolver: Methods and software description*, ACM Transactions on Mathematical Software, Vol. 37, No. 2, (2010), 21:1–21:30.

## Thermodynamics, gas phase Gibbs free energy¶

At the end of a completed IR Frequencies (normal modes) calculation, a survey is given of thermodynamic properties: entropy, internal energy, constant volume heat capacity, enthalpy and Gibbs free energy, see:

## Stress tensor¶

For periodic systems (chains, slabs, bulk) one can also request the clamped-ion
stress tensor (note: the clamped-ion stress is only part of the *true* stress
tensor):

```
Properties
StressTensor Yes/No
End
```

`Properties`

`StressTensor`

- Type
Bool

- Default value
No

- GUI name
Stress tensor

- Description
Calculate the stress tensor.

The clamped-ion stress tensor \(\sigma_\alpha\) (Voigt notation) is computed via numerical differentiation of the energy \(E\) WRT a strain deformations \(\epsilon_\alpha\) keeping the atomic fractional coordinates constant:

where \(V_0\) is the volume of the unit cell (for 2D periodic system \(V_0\) is the area of the unit cell, and for 1D periodic system \(V_0\) is the length of the unit cell).

The clamped-ion stress tensor (in Cartesian notation) is written to the engine
result file in `AMSResults%StressTensor`

.

## Elastic tensor¶

The elastic tensor \(c_{\alpha, \beta}\) (Voigt notation) is computed via second order numerical differentiation of the energy \(E\) WRT strain deformations \(\epsilon_\alpha\) and \(\epsilon_\beta\):

where \(V_0\) is the volume of the unit cell (for 2D periodic system \(V_0\) is the area of the unit cell, and for 1D periodic system \(V_0\) is the length of the unit cell).

For each strain deformation \(\epsilon_\alpha \epsilon_\beta\), the atomic positions will be optimized. The elastic tensor can be computed for any periodicity, i.e. 1D, 2D and 3D. Computing the elastic tensor will also calculate the bulk modulus, Young’s modulus, and shear modulus, as well as the Poisson ratio.

See also

To compute the elastic tensor, request it in the `Properties`

input block of
AMS:

```
Properties
ElasticTensor Yes/No
End
```

Note

The elastic tensor should be computed at the fully optimized geometry. One
should therefore perform a geometry optimization of all degrees of freedom,
**including the lattice vectors**. It is recommended to use a tight gradient
convergence threshold for the geometry optimization (e.g. Quality “Good”).
Note that all this can be done in one job by combining the geometry
optimization task with the elastic tensor
calculation.

The elastic tensor (in Voigt notation) is printed to the output file and stored
in the engine result file in the `AMSResults`

section (for 3D system, the elastic tensor in Voigt notation is a 6x6 matrix;
for 2D systems is a 3x3 matrix; for 1D systems is just one number).

Options for the numerical differentiation procedure can be specified in the
`ElasticTensor`

input block:

```
ElasticTensor
ConvergenceQuality [Normal | Good | VeryGood]
StrainStepSize float
End
```

`ElasticTensor`

- Type
Block

- Description
Options for numerical evaluation of the elastic tensor.

`ConvergenceQuality`

- Type
Multiple Choice

- Default value
Good

- Options
[Normal, Good, VeryGood]

- GUI name
Convergence

- Description
The tightness of the convergence of the geometry optimizations for each strain deformation. This should not be set higher than the overall convergence quality of the preceeding geometry optimization configured by the

`GeometryOptimization%Convergence%Quality`

keyword.

`StrainStepSize`

- Type
Float

- Default value
0.001

- Description
Step size (relative) of strain deformations used for computing the elastic tensor numerically.

Pressure or non-isotropic external stress can be included in your simulation via the corresponding engine addons.

The elastic tensor calculation supports AMS’ double parallelization, which can perform the calculations for the individual displacements in parallel. This is configured automatically, but can be further tweaked using the keys in the `NumericalDifferentiation%Parallel`

block:

```
ElasticTensor
Parallel
nCoresPerGroup integer
nGroups integer
nNodesPerGroup integer
End
End
```

`ElasticTensor`

`Parallel`

- Type
Block

- Description
Options for double parallelization, which allows to split the available processor cores into groups working through all the available tasks in parallel, resulting in a better parallel performance. The keys in this block determine how to split the available processor cores into groups working in parallel.

`nCoresPerGroup`

- Type
Integer

- GUI name
Cores per group

- Description
Number of cores in each working group.

`nGroups`

- Type
Integer

- GUI name
Number of groups

- Description
Total number of processor groups. This is the number of tasks that will be executed in parallel.

`nNodesPerGroup`

- Type
Integer

- GUI name
Nodes per group

- Description
Number of nodes in each group. This option should only be used on homogeneous compute clusters, where all used compute nodes have the same number of processor cores.

## Numerical differentiation options¶

The following options apply whenever AMS computes gradients, Hessians or stress tensors via numerical differentiation.

```
NumericalDifferentiation
NuclearStepSize float
StrainStepSize float
End
```

`NumericalDifferentiation`

- Type
Block

- Description
Define options for numerical differentiations, that is the numerical calculation of gradients, Hessian and the stress tensor for periodic systems.

`NuclearStepSize`

- Type
Float

- Default value
0.005

- Unit
Bohr

- Description
Step size for numerical nuclear gradient calculation.

`StrainStepSize`

- Type
Float

- Default value
0.001

- Description
Step size (relative) for numerical stress tensor calculation.

AMS may use symmetry (key `NumericalDifferentiation%UseSymmetry`

) in case of numerical differentiation calculations.
If symmetry is used only symmetry unique atoms are displaced.
Symmetry is only recognized if the starting geometry has symmetry.
Symmetry is only used for molecules if the molecule has a specific orientation in space, like that the z-axis is the main rotation axis.
If the GUI is used one can click the Symmetrize button (the star), such that the GUI will (try to) symmetrize and reorient the molecule.
There are some cases that even after such symmetrization, the orientation of the molecule is not what is needed for the symmetry to be used in case of numerical differentiation calculations. If that is the case or if key `NumericalDifferentiation%UseSymmetry`

is set to ‘False’, then no symmetry will be used.

The numerical differentiation calculation supports AMS’ double parallelization, which can perform the calculations for the individual displacements in parallel. This is configured automatically, but can be further tweaked using the keys in the `NumericalDifferentiation%Parallel`

block:

```
NumericalDifferentiation
Parallel
nCoresPerGroup integer
nGroups integer
nNodesPerGroup integer
End
End
```

`NumericalDifferentiation`

`Parallel`

- Type
Block

- Description
Options for double parallelization, which allows to split the available processor cores into groups working through all the available tasks in parallel, resulting in a better parallel performance. The keys in this block determine how to split the available processor cores into groups working in parallel.

`nCoresPerGroup`

- Type
Integer

- GUI name
Cores per group

- Description
Number of cores in each working group.

`nGroups`

- Type
Integer

- GUI name
Number of groups

- Description
Total number of processor groups. This is the number of tasks that will be executed in parallel.

`nNodesPerGroup`

- Type
Integer

- GUI name
Nodes per group

- Description
Number of nodes in each group. This option should only be used on homogeneous compute clusters, where all used compute nodes have the same number of processor cores.