Calculate Hessian for configurations that are not at the minimum of the Potential Energy Surface

Question

I want to calculate the Hessian matrix for molecular configurations that are not at the minimum of the Potential Energy Surface, using cp2k. But it seems that cp2k requires the configuration to be geometrically optimized. In the Vibrational Analysis section of the cp2k manual, it says

The analysis assumes a stationary state (minimum or TS), i.e. tight geometry optimization (MAX_FORCE) is needed as well.

What would happen if I just use cp2k to calculate Hessian for the not-geometrically-optimized configurations? Will it give the correct answer?

Answer 1

It is difficult to understand what you mean by "correct" answer, especially when we are discussing a "not-geometrically-optimized configuration" as you have called it, so I will try to go through the principles.

This matrix of all possible second derivatives of the energy with respect to the nuclear displacements of the atoms in your structure is referred to as the Hessian matrix. The Hessian can be diagonalized to obtain harmonic force constants, vibrational frequencies, and normal modes.

The vibrational frequencies of the normal modes and the eigenvectors give the amplitudes of motion along each the mass-weighted Cartesian coordinates that belong to each mode. If this same kind of analysis were performed at a geometry corresponding a "not-geometrically-optimized configuration" you would get negative modes which would essentially characterize a "downward" curvature of the energy surface. There is a peculiar case when you have exactly one negative mode, this would usually correspond to a transition state in a configuration or reaction analysis. Here is a good reference on this : https://simons.hec.utah.edu/ITCSecondEdition/chapter3.pdf

Finally, this means that if you do your calculations on your proposed non-optimized structure, you would get several negative frequencies. Does that make them correct or incorrect ? They would technically be correct for that structure itself. However, the structure itself is non-optimized. So maybe the best answer would be that the values are correct but have no real scientific value.

Another good reference : https://kthpanor.github.io/echem/docs/hessians.html

Answer 2

You could in principle compute the entire Hessian for any configuration of a molecule. The problem is that we are typically not interested in the entire Hessian: we want the Hessian for the internal coordinates/vibrations, with rotations and translations projected out.

To separate out translations and rotations, we can use the Eckhart conditions, which define a reference frame that translates/rotates with the molecule. These Eckhart conditions are not satisfied for a nonstationary point, i.e. a point where the nuclear gradients are not all zero. In this case, we can't cleanly separate rotations/translations from vibrations.

So you could compute the Hessian at a nonstationary point and compute it's eigenvectors/values, but away from a stationary point they won't have a clear connection to the vibrational modes/frequencies you are likely interested in computing.

Answer 3

The Hessian at a non-stationary point is well-defined, since it is the second derivative of the energy wrt. nuclear coordinates. The transformation to vibrational normal coordinates, however, assumes that the gradient is zero. In practice, this is never the case, but "sufficiently small" is normally good enough. At non-stationary points, one can define 3N-7 generalized vibrations by projecting out the translational and rotational degrees of freedom and the gradient. This have been used in "reaction path methods" for including the effects of nuclear motion perpendicular to the reaction coordinate, for example defined by the IRC path.