The aim of this post is to determine and visualize the vibrational modes of an isolated ammonia molecule. After a discussion of the theory of normal modes, the frequencies will be calculated directly using DMol3 in Materials Studio. Materials Studio also provides animations of the normal modes to help with visualization of the modes.
Normal Modes
Consider a classical system with generalized coordinates \( \{q_i \} \) for \( i=1\dots n\) moving under the influence of a potential energy \( V(q_i) \). Classical mechanics tells us that the equations of motion for this system are
\begin{equation} m_i \ddot{q}_i=-\frac{\partial V}{\partial q_i} \end{equation}
where \( m_i \) is the mass associated with the generalized coordiante \( q_i \). The fixed points of this system occur when \( \ddot{q_i}=0 \) for all \(i\) or equivalently,
\begin{equation} \frac{\partial V}{\partial q_i}=0 \end{equation}
Of course, this means that the fixed points of the system occur when the system is configured such that \(V\) is at a critical point. At such a point, the system can sit at equilibrium. Denote the configuration of the generalized coordinates that \( (2) \) holds as \( \{q^0_i\} \). If we suppose that our system is in a configuration close to \( \{q^0_i\} \), we can expand the potential in a Taylor series centered at \( \{q^0_i\} \) and truncate it at second order in \( q\):
\begin{equation} V=V_0+\frac{1}{2} \sum_{jk} \frac{\partial^2 V}{\partial q_j \partial q_k} (q_j – q^0_j)(q_k – q^0_k) \end{equation}
where the derivatives are evaluated at \( \{q^0_i\} \). Under such circumstances, the equations of motion reduce to a simple and familiar form:
\begin{equation} \ddot{q}_i=-\sum_j \frac{1}{m_i} \frac{\partial^2 V}{\partial q_i \partial q_j}(q_j – q^0_j) \end{equation}
These equations of motion describe \(n\) coupled simple harmonic oscillators. It is useful to introduce the following change of coordinates:
\[x_i=\frac{1}{\sqrt{m_i}}(q_i – q^0_i)\]
The equations of motion then become:
\begin{equation} \ddot{x}_i=-\sum_j \frac{1}{\sqrt{m_i m_j}} \frac{\partial^2 V}{\partial q_i \partial q_j}x_j \end{equation}
The solutions of the above equation form an \(n\)-dimensional vector space, and it is possible to construct an orthogonal basis for this vector space. Using the ansatz \( x_i(t)=x_i(0) \sin(\omega t+\delta) \), where \( \delta \) is a constant that depends on initial conditions, the equation becomes:
\begin{equation} \omega^2 x_i(0)=\sum_j \frac{1}{\sqrt{m_i m_j}} \frac{\partial^2 V}{\partial q_i \partial q_j}x_j(0) \end{equation}
If we think of \(x_i\) as vectors, this equation says that \(x_i(0)\) is an eigenvector of the \(n \times n\) symmetric matrix \(\frac{1}{\sqrt{m_i m_j}} \frac{\partial^2 V}{\partial q_i \partial q_j}\) with eigenvalue \(\omega^2\). A theorem of linear algebra guarantees that the eigenvectors of a symmetric matrix can be made pairwise orthogonal. Therefore the \(n\) eigenvectors of \(\frac{1}{\sqrt{m_i m_j}} \frac{\partial^2 V}{\partial q_i \partial q_j}\) form an orthogonal basis for the space of solutions of \((5)\). Denote these solutions as
\begin{equation} x^a_i(t)=e^a_i \sin(\omega_a t+\delta_a) ~~(a=1,\dots,n) \end{equation}
where \(a\) is a label for the eigenvalues and eigenvectors. Any solution of \((5)\) can be written as a unique linear combination of the \(x^a_i(t)\):
\begin{equation} y_i(t)= \sum_{a=1}^n C_a x^a_i(t) \end{equation}
The solutions \(x^a_i(t)\) are called normal modes of the system and the \(\omega_a\) are called normal frequencies. Equation \((8)\) highlights for us the use of normal modes: given initial conditions of the system, we can expand in terms of the normal modes to find a solution compatible with the initial conditions. It is worth remembering, however, that this approach is only valid when the coordinates are close to the fixed points of \(V\).
Determing the normal frequencies
To determine the normal frequencies, we need to know the partial derivatives of the potential energy evaluated at the critical points. This can be accomplished using DFT calculated energies for a molecular system as follows:
\begin{equation} \frac{\partial^2 V}{\partial q_i \partial q_j} \approx \frac{E(q^0_i+\delta q_i,q^0_j+\delta q_j)-E(q^0_i+\delta q_i,q^0_j-\delta q_j)-E(q^0_i-\delta q_i,q^0_j+\delta q_j)+E(q^0_i-\delta q_i,q^0_j-\delta q_j)}{4(\delta q_i)(\delta q_j)} \end{equation}
This formula is exact in the limit \(\delta q_i,\delta q_j \rightarrow 0\). After determining the derivatives of \(V\) we can diagonalize the matrix \(\frac{1}{\sqrt{m_i m_j}} \frac{\partial^2 V}{\partial q_i \partial q_j}\) to determine its eigenvalues. Note that nothing restricts the sign of these eigenvalues and in fact only \(\omega_a^2>0\) correspond to oscillatory motion. Materials Studios has a function to compute these positive eigenvalues from a geometry optimization calculation. These are the modes of ammonia we will focus on in this post.
Modes of ammonia
Figure 1: A relaxed ammonia molecule NH3. The atoms are arranged in a trigonal pyramidal shape.
We performed a DMol3 geometry optimization calculation with vibrational analysis in Materials Studio to determine the oscillatory frequencies of an isolated ammonia molecule. The functional used was GGA-PBE. The energy convergence tolerance was \(10~\mu Ha\), the max. force was .002 Hartrees per Angstrom, and the max. displacement was .005 Angstroms. For this system, \(n=12\) and there are 6 positive eigenvalues of interest. The results, along with the comparisons with experimental values, are shown below:
| Mode | DFT Result (1/cm) |
|---|---|
| N-H wagging | 1064.55 |
| H-N-H scissoring (1) | 1631.18 |
| H-N-H scissoring (2) | 1631.18 |
| N-H symmetric stretch | 3390.54 |
| N-H asymmetric stretch (1) | 3524.06 |
| N-H asymmetric stretch (2) | 3524.06 |
Materials Studio also gives animations of these modes. Click on any of the images below to view an animation of the normal modes.
Figure 2: N-H wagging
Figure 3: H-N-H scissoring (1)
Figure 4: H-N-H scissoring (2)
Figure 5: N-H symmetric stretch
Figure 6: N-H asymmetric stretch (1)
Figure 7: N-H asymmetric stretch (2)
References
[1] PBE functional: Perdew Burke Ernzerhof: Phys. Rev. Lett. 77, 3865 (1996)
[2] The Generation and Use of Delocalized Internal Coordinates in Geometry optimization;
Baker Kessi Delley: J. Chem. Phys., 105, 192 (1996)
Andzelm Fitzgerald King-Smith: Chem. Phys. Lett., 335, 321 (2001)
[3] Parallel solution of partial symmetric eigenvalue problems from electronic structure calculations:
Auckenthaler, Blum, Bungartz, et al.: Parallel Computing 37, 783 (2011)
