NormalModes.jl

This is a simple package that provides utility functions to compute normal modes from a Hessian matrix.

Additionally, this package gives me a space to ramble about normal modes and I am shamelessly use this README for this purpose.

A proper introduction to normal modes

I find the subject confusing and the introductions describing the problem were of little help to me. Therefore, I think it is worth laying down the basics here.

We are interested in the ground state of an equilibrium Hamiltonian $\mathcal{H} = \mathcal{T} + \mathcal{V}$, whith $\mathcal{T}$ the kinetic energy operator and $\mathcal{V}$ the potential energy operator. We note ${\rm \bf x}$ the vector of concatenated position, that is ${\rm \bf x} = {\rm vcat}({\rm \bf r}_1, {\rm \bf r}_2, {\rm \bf r}_3, \dots)$¹, with ${\rm \bf r}_A$ the position of atoms $A$.

Then we do a Taylor expansion of the potential, leading to

$$ \mathcal{V}({\rm \bf x}) = V_0 + \nabla \mathcal{V}({\rm \bf x}) + \frac{1}{2} {\rm \bf x}^T {\rm \bf H} {\rm \bf x} + \mathcal{O}({\rm \bf x}^3), $$

where ${\rm \bf H}$ is the Hessian matrix of the potential², and $V_0$ a constant shift that we set to zero without loss of generality. Moreover, we are assuming to be at equilibrium, therefore the gradient $\nabla \mathcal{V}({\rm \bf x})$ is zero (we are a local minimum of the potential energy). This leads to the approximation

$$ \mathcal{V}({\rm \bf x}) = \frac{1}{2} {\rm \bf x}^T {\rm \bf H} {\rm \bf x}. $$

With this potential the Hamiltonian becomes

$$ \mathcal{H} = \sum_i \frac{-\hbar^2}{2 m_i} \frac{\partial^2}{\partial x_i^2} + \frac{1}{2} {\rm \bf x}^T {\rm \bf H} {\rm \bf x}, $$

where $m_i$ is the mass of dimension $i$³.

This is the Hamiltionian of a system of coupled harmonic oscillators. The standard procedure is here to change basis and rewrite it as a system of uncoupled harmonic oscillators. To this end, we introduce the diagonal mass weighting matrix ${\rm \bf M}$ with ${\rm \bf M}_{ii} = m_i^{-\frac{1}{2}}$. Then we perform an eigendecomposition of the weighted Hermitian $\rm \bf M H M$ as

$$ {\rm \bf M H M} = {\rm \bf U \Omega U}^T, $$

where $\rm \bf U$ is the matrix of eigenvectors and $\rm \bf \Omega$ the diagonal matrix of eigenvalues. Since we are looking at a stable equilibrium, all eigenvalues must be positive and we write them as ${\rm \bf \Omega}_{ii} = \omega_i^2$⁴.

Using ${\rm \bf M}$ and ${\rm \bf U}$ we can define the new coordinates $\rm \bf z$ as

$$ {\rm \bf x} = {\rm \bf M U z}. $$

Inserting it in the Hamiltonian we get the following decoupled system of harmonic oscillators (see the annex for the derivation)

$$ \mathcal{H} = \sum_j \frac{-\hbar^2}{2} \frac{\partial^2}{\partial z_j^2} + \frac{1}{2} \omega_j^2 z_j^2. $$

From this expression, we can get the normal modes and all associated quantities.

The package

The package aims at providing an interface to normal modes, but what do we mean by normal modes exactly? There are unfortunately multiple reasonnable candidates, namely

1. The eigen vectors ${\rm \bf u}_j$ which are the columns of $\rm \bf U$. They are normed and orthogonal, but live in the mass weighted space.
2. The columns of $\rm \bf MU$. They live in real space, but they are neither normed nor orthogonal.
3. The columns of $\rm \bf MU$ normalized. They are still not orthogonal, but at least they have norm 1. In addition, all information about $\rm \bf M$ is removed, which is often needed for further calculations. So those are ofter return together with so-called mode masses, the norm of the columns of $\rm \bf MU$⁵.

As far as I understand, the last ones are the quantities known as normal modes. Thankfully, to avoid users a slow descent into madness, the package provides a specific object type NormalDecomposition(hessian), that stores the useful information and encapsulate the inner confusing parts. The available methods acting on a NormalDecompoisition are

• normal_modes : The normal modes (in real space). They are useful to see the effect off the vibration on the geometry, for example when plotting them.
• frequencies : Frequencies of the modes (in GHz).
• wave_numbers : Wave numbers of the modes (in inverse cm). Sometimes those are also called "frequencies", as they are directly proportional to each other.
• reduced_mass : Reduced masses of the modes (in AMU), following this question and (I believe) similar to what Gamess does.
• sample : Perform Wigner sampling, and get the displacements from the mean for positions and momenta.

All the values are given as Unitful quantities.

Tools

Data can be loaded from an ORCA hess file, using load_orca(filename, NormalDecomposition).

Appendix

Derivation of the uncoupled equation

We want to substitute ${\rm \bf x} = {\rm \bf M U z}$ into

$$ \mathcal{H} = \sum_i \frac{-\hbar^2}{2 m_i} \frac{\partial^2}{\partial x_i^2} + \frac{1}{2} {\rm \bf x}^T {\rm \bf H} {\rm \bf x}. $$

By the definition of the diagonal matrix ${\rm \bf M}$ and the orthogonal matrix ${\rm \bf U}$ we have

$$ \frac{1}{2} {\rm \bf x}^T {\rm \bf H x} = \frac{1}{2} {\rm \bf z}^T {\rm \bf U}^T {\rm \bf M} {\rm \bf HMUz} = \frac{1}{2} {\rm \bf z}^T {\rm \bf \Omega z}, $$

where ${\rm \bf \Omega}$ is diagonal.

For the other term, we use Leibnitz chain rule, starting with the first derivative only

$$ \frac{\partial}{\partial x_i} = \sum_j \frac{\partial z_j}{\partial x_i} \frac{\partial}{\partial z_j}. $$

Next we use ${\rm \bf z} = {\rm \bf U}^T {\rm \bf M}^{-1} {\rm \bf x}$ to compute

$$ \frac{\partial z_j}{\partial x_i} = \frac{\partial}{\partial x_i} {\rm \bf e}_j^T {\rm \bf U}^T {\rm \bf M}^{-1} {\rm \bf x} = {\rm \bf e}_j^T {\rm \bf U}^T {\rm \bf M}^{-1} {\rm \bf e}_i, $$

where ${\rm \bf e}_i$ are the unit vectors of the canonical basis. Since no new dependency in ${\rm \bf x}$ or ${\rm \bf z}$ appear, we can apply it twice. We transpose one of them however, to make the vector with index $i$ appear together⁶

$$ \sum_i \frac{-\hbar^2}{2 m_i} \frac{\partial^2}{\partial x_i^2} = -\frac{\hbar^2}{2} \sum_{i, j, k} \frac{1}{m_i} \frac{\partial z_k}{\partial x_i} \left[\frac{\partial z_j}{\partial x_i} \right]^T = -\frac{\hbar^2}{2} \sum_{j, k} {\rm \bf e}_k^T {\rm \bf U}^T {\rm \bf M}^{-1} \left[ \sum_i \frac{1}{m_i} {\rm \bf e}_i {\rm \bf e}_i^T \right] {\rm \bf M}^{-1} {\rm \bf U} {\rm \bf e}_j \frac{\partial}{\partial z_j} \frac{\partial}{\partial z_k}. $$

Now it would be really beautiful if

$$ \sum_i \frac{1}{m_i} {\rm \bf e}_i {\rm \bf e}_i^T = {\rm \bf M}^2, $$

as it would make everything collapse. Thankfully with our choice of ${\rm \bf M}$ it is exactly what happens.

Putting it in and using the fact that ${\rm \bf e}_j^T {\rm \bf e}_k$ is a Kroenecker delta, we get, as expected,

$$ \sum_i \frac{-\hbar^2}{2 m_i} \frac{\partial^2}{\partial x_i^2} = -\frac{\hbar^2}{2} \sum_j \frac{\partial^2}{\partial z_j^2}. $$

Footnotes

Footnotes

1. I am indeed using the julia function vcat to describe this mathematical operation, I have never found a better way to express it. ↩

2. Which is well defined because in the position representation the potential is just a scalar (and not an operator). ↩

3. Which is the mass of the atom from which this component comes from. In other words, in the mass vector, each mass is repeated 3 times, once for each spatial dimension. ↩

4. Eigenvalue are in theory positive or zero, because the Hamiltonian is semi-positive definite and symmetric. Numerically the smallest eigenvalues can be negative instead of exactly zero. ↩

5. Which is pure lunacy if you ask me. ↩

6. We transpose a scalar, which is a trick I like very much and which is surprisingly useful. ↩