Harmonic Approximations

Let $\mathbf{x} \in \mathbb{R}^{N\times 3}$ denote the Cartesian coordinates of a molecular structure, at thermodynamic equilibrium, the probability of observing such a particular structure follows a Boltzmann distribution

where $U: \mathbb{R}^{N \times 3} \to \mathbb{R}$ is the potential energy of the system. Since, an accurate description of $U$ is generally not feasible a full description is generally not possible. Around modes, also known as conformations, however the potential energy surface (PES) looks locally Gaussian. This essentially follows form a Taylor expansion. Let $\bm{\mu}$ denote a local minimum of the PES, then around $\bm{\mu}$ we can approximate $U(\mathbf{x})$ by

where the second term vanishes because we are at a mode of the PES.

For ease of notation we will from now on assume that $\bm{\mu}$ is at the center of the coordinate system. That is we transform the coordinates such that $\mathbf{x} = \mathbf{x}' - \bm{\mu}$, where $\mathbf{x}'$ are the original coordinates.

As we can see, the results potential is a simple harmonic potential, giving rise to the name Harmonic Approximation. We will now see, that because of the approximation, we can find a linear transformation that decouples the coordinates of the system, i.e. makes them independent. As a result we can analyze the whole system by looking at individual dimensions.

Let $\mathbf{c} = \mathbf{P}\mathbf{x}$, for the yet to be found transformation. The full energy of the approximated system can then be written as

where we have used $\mathbf{M} = \text{diag}(\mathbf{m})$ to denote the mass matrix with entries $\mathbf{m}_i$ on its diagonal.

For the dimensions of $\mathbf{c}$ to be independent we require $\mathbf{P}^T \mathbf{M} \mathbf{P}$ and $\mathbf{P}^T H \mathbf{P}$ to be simultaneously diagonalizable.

In other words, we want to find $\mathbf{P}$, such that

for some diagonal $\mathbf{D}, \bm{\lambda}$. Since $\mathbf{M}$ is already diagonal, we can make a choice for $\mathbf{D}$, which in turn will determine $\mathbf{P}$. A common choice is simply $\mathbf{D} = \mathbf{I}$, such that $\mathbf{P} = \mathbf{M}^{-\frac{1}{2}}$. Note that this choice however does not necessarily give rise to diagonal $\bm{\lambda}$. Observe, that we can satisfy $\mathbf{P}^T \mathbf{M} \mathbf{P} = \mathbf{I}$, also if we choose $\mathbf{P} = \mathbf{M}^{-\frac{1}{2}} \mathbf{Q}$ for any orthogonal $\mathbf{Q}$.

With this choice we derive the condition

The matrix $\tilde{\mathbf{H}}$ is known as the mass-weighted Hessian or dynamical matrix in the literature. Because the weights are positive, and the Hessian is always symmetric, this is a real symmetric matrix and by the spectral theorem guarantees that a $\mathbf{Q}$, satisfying this constraint exists. Moreover, we know that $\mathbf{Q}$ will contain the eigenvectors of $\tilde{\mathbf{H}}$ as rows and $\bm{\lambda}$ has the eigenvalues on the diagonal. Hence to find $\mathbf{Q}$, we need to solve the eigenvalue problem

Using this insight we can write the energy of the system as

This is exactly the energy of a system of $N$ independent harmonic oscillators, which are probably to most well understood system in all physics.

For example, observe that the eigenvalues of the mass-weighted Hessian also tell us which modes (eigenvectors) correspond to stiff directions (large change in potential energy) and soft directions (large displacements are possible). Hence, we can make use of this knowledge, to ignore stiff directions when moving on the PES.

Rotations, Translations and Symmetries¶

Something that is not immediately obvious from the above analysis, is that in general 6 of the eigenvalues will be exactly zero. This follows from the fact that rotations and translations do not change the potential energy of the system.

Since in three dimensions these groups have only 3 elements each,

More generally, as we have some function $f(\mathbf{x})$ that is invariant under some transformation $g$, i.e. $f(g(\mathbf{x})) = f(\mathbf{x})$.