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

The 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
applications, see the following pages in the manual:

* IR frequencies

* PES point character


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

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

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

  3. ... 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 criterion 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:

* IR frequencies

  * Thermodynamics

  * Gibbs free energy change for a gas phase reaction


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:

   \sigma_\alpha= \frac{1}{V_0} \left. \frac{\partial E}{\partial
   \epsilon_\alpha} \right|_\text{constant atomic fractional
   coordinates}

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:

   c_{\alpha, \beta} = \frac{1}{V_0} \frac{\partial^2 E}{\partial
   \epsilon_\alpha \partial \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: Example: Elastic tensor

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