How to calculate wavenumbers of normal modes from the eigenvalues of the Cartesian force constant matrix?

I would like to get the wavenumbers (Frequencies -- in Gaussian) from the eigenvalues of the force constant matrix (Force constants in Cartesian coordinates: in Gaussian).

Could someone please explain how to perform this calculation? The Cartesian force constants are given in Hartrees per Bohr squared, while the wavenumbers are in $\pu{cm^-1}$.

• In short: sirius.chem.vt.edu/wiki/… – Wildcat May 16 '17 at 14:53
• @Wildcat That's a good start, but it leaves out a lot of details -- why mass-weight? How does the units conversion work? Answer pending as soon as I have the time to write it. – hBy2Py May 16 '17 at 14:57
• Am I being just too naive to suggest that for a harmonic oscillator $E=k\Delta x^2/2$ and so k in units of Hartree/Bohr^2 are to be expected as the same as $\pu{J/m^2}$ (with conversion constants) and that $2\pi c\nu = \sqrt{k/\mu}$ has $\nu$ in $\pu{cm^{-1}}$ ? – porphyrin May 16 '17 at 17:53
• @porphyrin No, you're right, the analogy is fairly direct. The mechanics of the calculation and the specifics of the units transformation are highly non-obvious, though. As noted, answer is pending. Possibly not until the weekend -- things are kind of hectic ATM. – hBy2Py May 16 '17 at 18:09

Short answer: (1) mass weight each force constant (see below). (2) Diagonalize the matrix. (3) Convert the eigenvalues ($\{ k_1^{\prime\prime} \}$) to wavenumbers $\tilde{\nu}=\sqrt{k_i^{\prime\prime}}/(2\pi \times 137 \times 5.29\times 10^{-9})$. 137 is the speed of light in atomic units and $5.29\times 10^{-9}$ converts bohr to cm.

More detail

Harmonic vibrational frequencies come from solving $F=ma$ assuming $E(x)=\frac{1}{2}k(x-x_0)^2 = \frac{1}{2}k^\prime l^2$ where $l=\sqrt{m}(x-x_0)$ and $k^\prime = k/m$. $F=-(\partial E/\partial x)$ so $F=ma$ becomes $$\frac{\partial^2 l}{\partial t^2}+k^\prime l = 0$$ The solution to this differential equation is $$l = K \cos(\sqrt{k^\prime}t) +\varepsilon$$ where $K$ and $\varepsilon$ are arbitrary constants.

The solution tells us that the displacement coordinate will increase and decrease in a repeating fashion every $\sqrt{k^\prime}/2\pi$ times, which is the definition of frequency ($\nu$). The units of $k^\prime$ are energy/(mass length$^2)$ or 1/time$^2$, because energy = (mass)(length/time)$^2$. So the units of $\nu$ is 1/time. This can be converted to a wavenumber ($\tilde{\nu}$, 1/length) by dividing $\nu$ by the speed of light.

Diatomic molecule

You have 6 Cartesian coordinates $\{q_i\}$ and 6x6 Cartesian force constants (some have the same value), $$k_{ij} = \left( \frac{\partial^2 E}{\partial q_i \partial q_j} \right)$$ and it's not clear what to do with the cross terms. So you diagonalize the matrix to get rid of the cross terms, but you need to mass-weight first for the reasons described above. $$k^\prime_{ij} = \frac{1}{\sqrt{m_i m_j}}k_{ij}$$ After diagonalization of $\mathbf{k^\prime}$ you have 3N-5(6) non-zero force constants $\{ k_i^{\prime \prime}\}$ , which you can convert to wavenumbers by $\sqrt{k_i^{\prime \prime}}/2\pi c$

• Thank you! Could you please help me, how to determine the normal mode which converts the minimum into the transition state numerically? (This normal mode corresponds to the one of the imaginary frequency.) – Roloka May 17 '17 at 7:50

This calculation is extraordinarily common in computational chemistry packages, but despite that I had an extraordinarily hard time finding a clear exposition of the method. I ended up having to piece it together myself, drawing from a couple of papers by Vincenzo Barone$^{1,2}$, an arXiv submission of George Phillies, the classic Wilson, Decius & Cross, and the method as described on the Gaussian.com website. The below exposition is written for the analysis of a nonlinear molecule; suitable adjustments to the formulae must be made for linear systems.

The first thing to do is to convert the Cartesian Hessian (force constant matrix) $\require{begingroup}\begingroup\def\mbf#1{\mathbf{#1}}\mbf{F}$ to vibrational normal coordinates; the first step of this process is to mass-weight it. This is necessary in order to establish the proper physical relationship between the frequencies of motion and the force constants in the normal-mode system of displacements. A useful analogy can be constructed to a one-dimensional mechanical system with two equal masses and three springs of equal force constant; for this system, the defining equation is:

$$\left(\omega^2m - 2k \right)^2-k^2 = 0, \tag{1}\label{det-eqn}$$

where $\omega$ is the vibrational frequency, each object is of mass $m$, and each spring has force constant $k$. As can be seen in Eq. $\eqref{det-eqn}$, the units of $\omega^2$ are related to those of $k$ by a factor of $m^{-1}$.

For vibrational analysis, this adjustment is accomplished through the use of the expanded atomic mass matrix, $\mbf M$. As there are $3N$ degrees of freedom to a molecular system with $N$ atoms, the dimension of $\mbf F$ is $\left(3N\times3N\right)$. Accordingly, the matrix $\mbf M$ used for mass-weighting must have the same dimensions. The elements of $\mbf F$ are usually clustered such that:

$$\def\pardsq#1#2{{\partial^2 #1 \over \partial #2 ^2}} \def\pardxy#1#2#3{{\partial^2 #1 \over \partial #2 \partial #3}} \mbf F = \left[ \begin{array}{cccc} \pardsq{E}{x_1} & \pardxy{E}{x_1}{y_1} & \cdots & \pardxy{E}{x_1}{z_N} \\ \pardxy{E}{y_1}{x_1} & \pardsq{E}{y_1} & \cdots & \pardxy{E}{y_1}{z_N} \\ \vdots & \vdots & \ddots & \vdots \\ \pardxy{E}{z_N}{x_1} & \pardxy{E}{z_N}{y_1} & \cdots & \pardsq{E}{z_N} \end{array} \right], \tag{2}\label{f-mtx}$$

where $x_i$ is the $x$-coordinate of atom $i$ in the system. Therefore, $\mbf M$ is usually constructed from the atomic masses $m_i$ as the diagonal matrix:

$$\mbf M = \left[ \begin{array}{ccccccc} m_1 & 0 & 0 & 0 & \cdots & 0 & 0\\ 0 & m_1 & 0 & 0 & \cdots & 0 & 0\\ 0 & 0 & m_1 & 0 & \cdots & 0 & 0 \\ 0 & 0 & 0 & m_2 & \cdots & 0 & 0\\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & m_N & 0 \\ 0 & 0 & 0 & 0 & \cdots & 0 & m_N \\ \end{array} \right] \tag{3}\label{m-mtx}$$

The mass-weighted Hessian $\mbf f$ is then calculated as the following:

$$\mbf f = \mbf M^{-1/2} \mbf F \mbf M^{-1/2} \tag{4}\label{weighted-f-mtx}$$

This results in each element $f_{ij}$ being divided by the product $\sqrt{m_am_b}$, where $m_a$ and $m_b$ are the atomic masses of the atoms whose geometric coordinates are involved in defining that element $f_{ij}$. In particular, it's important to note here that each factor of $m_a^{1/2}$ is weighting a matching length unit associated with atom $a$. To see this, note that each $F_{ij}$ is a second derivative of the energy with respect to the displacements of two atoms $a$ and $b$, where $a=b$ is possible, with units of Hartrees per Bohr squared $\left(\pu{E_h}/\pu{B2}\right)$. The weighting approach effected by Eq. $\eqref{weighted-f-mtx}$ means that $f_{ij} = \pardxy{E}{w_im_a^{1/2}}{w_jm_b^{1/2}}$, where atoms $a$ and $b$ are associated with displacements $w_i$ and $w_j$, respectively, and $\overline w=\left[x_1, y_1, z_1, x_2, \ldots \right]$. Thus, it is actually the displacements that are being mass-weighted, each by the square root of the mass of its associated atom.

At this point, in an ideal world, we could just diagonalize $\mbf f$ to obtain the eigenvalues we need (analogous to $\omega^2$ in the above example) to calculate the vibrational frequencies. This diagonalization yields the diagonal matrix $\mbf \Lambda$ containing the eigenvalues, with the accompanying eigenvectors matrix $\mbf L$:

$$\mbf \Lambda= \left[ \begin{array}{cccc} \lambda_1 & 0 & \cdots & 0 \\ 0 & \lambda_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \lambda_{3N} \end{array} \right] ~;~~ \mbf L= \left[ \begin{array}{cccc} | & | & & | \\ \overline \ell_1 & \overline \ell_2 & \cdots & \overline \ell_{3N} \\ | & | & & | \end{array} \right] \tag{5}\label{eigen-noproj}$$

Since molecules are isolated systems, in theory the energy should be independent of both translational and rotational motion. In this ideal circumstance, for a frequency analysis of a species at its equilibrium geometry, the first six eigenvalues $\lambda_1$ through $\lambda_6$ would be identically zero, corresponding to the absence of an effect on the second derivative of the energy from translational and rotational motions.

However, in reality, energy, gradient and Hessian calculations can never be perfectly converged and so small amounts of "contaminating" translational and rotational dependence will be present in the original Hessian $\mbf F$. This leads to nonzero, though usually rather small, values for $\lambda_1$ to $\lambda_6$. In order to "purify" the Hessian of these non-physical numerical artifacts, a projection operation is typically applied. A projection matrix $\mbf D$ is constructed for this purpose, composed of defined "normal modes" for translation and rotation, and trial normal mode vectors for the vibrational motion.

The "normal modes" of translational motion are straightforward to define. If a molecule is translating, then all of the atoms are moving in the same direction. The simplest way to express this in 'normal modes' form is by using translation along the coordinate axes $e_x, e_y, e_z$:

$$\def\vbuf#1{\left(\overset{~}{\underset{~}{#1}}\right)} \def\vecex{\left(\begin{array}{c}1 \\ 0 \\ 0 \end{array}\right)} \overline d_1= \left[\begin{array}{c}\vbuf{\vec e_x}_1 \\ \vbuf{\vec e_x}_2 \\ \vdots \\ \vbuf{\vec e_x}_N \end{array}\right] = \left[\begin{array}{c}\vecex_1 \\ \vdots \\ \vecex_N \end{array}\right] ~;~~ \overline d_2= \left[\begin{array}{c}\vbuf{\vec e_y}_1 \\ \vbuf{\vec e_y}_2 \\ \vdots \\ \vbuf{\vec e_y}_N \end{array}\right] ~;~~ \overline d_3= \left[\begin{array}{c}\vbuf{\vec e_z}_1 \\ \vbuf{\vec e_z}_2 \\ \vdots \\ \vbuf{\vec e_z}_N \end{array}\right] \tag{6}\label{trans-d-vecs}$$

Any set of three orthonormal vectors $e_1, e_2, e_3$ would work, however. For reasons I won't go into in detail here, the "normal mode" vectors for rotational motion can be constructed in a roughly similar fashion, using the cross products of the three orthonormal vectors (representing the three independent axes of rotation) with the displacement vectors $\vec r_a$ of each atom $a$ from the center of mass of the system:

$$\def\vbuf#1{\left(\overset{~}{\underset{~}{#1}}\right)} \def\vecex{\left(\begin{array}{c}1 \\ 0 \\ 0 \end{array}\right)} \overline d_4= \left[\begin{array}{c}\vbuf{\vec e_x \times \vec r_1} \\ \vbuf{\vec e_x \times \vec r_2} \\ \vdots \\ \vbuf{\vec e_x \times \vec r_N} \end{array}\right] ~;~~ \overline d_5= \left[\begin{array}{c}\vbuf{\vec e_y \times \vec r_1} \\ \vbuf{\vec e_y \times \vec r_2} \\ \vdots \\ \vbuf{\vec e_y \times \vec r_N} \end{array}\right] ~;~~ \overline d_6= \left[\begin{array}{c}\vbuf{\vec e_z \times \vec r_1} \\ \vbuf{\vec e_z \times \vec r_2} \\ \vdots \\ \vbuf{\vec e_z \times \vec r_N} \end{array}\right] \tag{7}\label{rot-d-vecs}$$

The corresponding columns of $\mbf D$ are obtained by normalizing each mass-weighted $\overline d_i$:

$$\overline D_i = {\mbf M^{1/2}\overline d_i \over \left| \mbf M^{1/2}\overline d_i \right|} \tag{8}\label{D-mtx-assy}$$

Note that the mass-weighting has a power of positive one-half here, since each element of a $\overline D_i$ indicates the contribution of the associated atomic displacement to the changes in the inertial configuration of the system when it undergoes deformation along mode i. To note, the normalization of Eq. $\eqref{D-mtx-assy}$ eliminates the dimensionality of this subset of $\mbf D$, converting its elements into unitless quantities.

The remainder of $\mbf D$ is filled with $3N-6$ arbitrary 'trial' normal vibrational mode vectors, which must be mutually orthonormal as well as orthogonal to the translational and rotational "modes" $\overline D_1$ through $\overline D_6$. This can be accomplished by seeding the remainder of the matrix with random values, checking to ensure linear independence, and performing any suitable orthonormalization procedure such as Gram-Schmidt.

By construction, $\mbf D$ is an orthogonal matrix, and thus can be used as a coordinate transform (a 'projection') from the Cartesian force constant matrix $\mbf f$ to a normal-mode force constant matrix $\mbf{\hat f}$:

$$\mbf{\hat f} = \mbf D^\mathrm T \mbf f \mbf D \tag{9}$$

At this point, all elements in rows 1-6 and columns 1-6 of $\mbf{\hat f}$ should be small, reflecting the translational and rotational artifacts present in the original Hessian. Only the elements $\hat f_{i,j}$ where both $i>6$ and $j>6$ should have appreciable magnitudes. We can now coerce the spurious translational and rotational content to zero and extract the lower-right sub-matrix of $\mbf{\hat f}$ that contains only the vibrational information of interest:

$$\mbf{\hat f}\Rightarrow \left[ \begin{array}{cc} \mbf{0} & \mbf{0} \\ \mbf{0} & \mbf{\hat f'} \end{array} \right]~;\quad \mbf{\hat f'}= \left[ \begin{array}{cccc} \hat f_{7,7} & \hat f_{7,8} & \cdots & \hat f_{7,3N} \\ \hat f_{8,7} & \hat f_{8,8} & \cdots & \hat f_{8,3N} \\ \vdots & \vdots & \ddots & \vdots \\ \hat f_{3N,7} & \hat f_{3N,8} & \cdots & \hat f_{3N,3N} \\ \end{array} \right] \tag{10}$$

Now we can perform the same diagonalization on $\mbf{\hat f'}$ as was used to generate Eq. $\eqref{eigen-noproj}$, obtaining the translation- and rotation-purified eigenvalues and eigenvectors for the vibrational modes:

$$\mbf \Lambda'= \left[ \begin{array}{cccc} \lambda_7' & 0 & \cdots & 0 \\ 0 & \lambda_8' & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \lambda_{3N}' \end{array} \right] ~;~~ \mbf L'= \left[ \begin{array}{cccc} | & | & & | \\ \overline \ell'_7 & \overline \ell'_8 & \cdots & \overline \ell'_{3N} \\ | & | & & | \end{array} \right] \tag{11}\label{eigen-proj}$$

With eigenvalues $\mbf \Lambda$ and/or $\mbf \Lambda'$ in hand, we now can turn to the actual calculation of the wavenumbers of the vibrational frequencies. Consider a representative eigenvalue $\def\hala{\hat\lambda}\hala$, which carries units of:

$$\def\uniteq{\enspace\left[=\right]\enspace} \hala \uniteq {\pu{E_h}\over\pu{amu B^2}}.\tag{12}\label{eigen-units}$$

The Hartree is the atomic unit of energy. Just as we can decompose the Joule into its composite fundamental units, we can do the same with the Hartree $(\pu{m_e}$ is the electron mass, and $\pu{T_a}$ is the atomic time unit, approximately $\pu{2.4E-17 s})$:

$$\pu J = {\pu{kg m2}\over\pu{s2}}~;~~\pu{E_h} = {\pu{m_e B2}\over\pu{T_a^2}} \tag{13}\label{energy-units}$$

One of the most subtle aspects of figuring out how to perform this wavenumber calculation correctly was to note, as shown above, that the atomic unit of mass is the electron mass, and not the atomic mass unit, despite the near identity of the terminology. Thus, the units of $\hat\lambda$ are actually:

$$\require{cancel} \hat\lambda \uniteq {\pu{E_h}\over\pu{amu B2}} = {\pu{m_e} \cancel{\mathrm{B^2}} \over \pu{T_a^2}} \cdot {1\over \pu{amu} \mathrm{\cancel{B^2}}} = {\pu{m_e}\over\pu{amu}} \cdot {1\over\pu{T_a^2}} \tag{14}\label{lambda-units}$$

For example, say one of the eigenvalues exhibited by a system has a value of $\hala = \pu{0.1000 E_h/amu B2}$. The value of $\hala$ in fundamental atomic units is thus:

\begin{align}\hala = 0.1000\, {\pu{E_h}\over\pu{amu B2}} &= 0.1000 \cdot {\pu{m_e}\over\pu{amu}} \cdot {1\over\pu{T_a^2}} \\ &= 0.1000 \cdot \pu{5.4858E-4} \cdot {1\over\pu{T_a^2}} = \pu{5.4858E-5 T_a^-2}\tag{15}\label{hala-Ta-val}\end{align}

The factor of $\pu{5.4858E-4}$ is the ratio $\pu{m_e / amu}$ per, e.g., NIST reported values.

The angular frequency of oscillation of this hypothetical vibrational mode is then just $\def\hanu{\hat\nu}\def\haom{\hat\omega}\def\brenu{\breve\nu} \haom = \sqrt{\hala}$, or $\haom = \pu{0.007407 T_a^-1}$. The cyclic oscillation frequency can be obtained via dividing by the usual factor of $2\pi$:

$$\hanu = {\haom\over2\pi} = 0.001179\,{\pu{cyc}\over\pu{T_a}}, \tag{16}\label{hanu-ta-val}$$

where $\pu{cyc}$ ("cycles") is included in the units as a reminder that $\hanu$ is a cyclic frequency.

Per standard physical chemistry derivations (e.g., at Wikipedia), the spacing between the energy levels of a quantum harmonic oscillator is equal simply to $\hbar\haom$ (or, equivalently, $h\hanu$). This represents the energy of photon required to excite the system along the associated vibrational mode $\hat\ell$. Per the normal expression for photons of $E=h\brenu$, the straightforward relationship of $\hanu=\brenu$ is obtained.

Thus, all that remains is to convert $\brenu$ to the desired wavenumber units of $\pu{cyc/cm}$. As noted by Jan Jensen, the speed of light expressed in atomic units has the value $c = \pu{137 B/T_a}$; combining this value for $c$ with a conversion from $\pu B$ to $\pu{cm}$ yields:

\begin{align} \brenu &= 0.001179\,{\pu{cyc}\over\pu{T_a}} \\ {\brenu\over c} &= \pu{8.604E-6}\,{\pu{cyc}\over\pu{B}} = 1626\,{\pu{cyc}\over\pu{cm}} \end{align}\tag{17}\label{final-result}

Thus, in the end, the conversion is as simple as that described by Jan. The above provides deeper insight into why that calculation gives the correct answer.

EPILOGUE: If I understand correctly, Gaussian prints the wavenumber values for the vibrational frequencies both before $(\lambda_i)$ and after $(\lambda_i')$ the rotational and translational noise is projected out of $\mbf F$. This is why two similar but non-identical sets of frequencies may appear in the output.

$^1$ Barone, V. J Chem Phys 122(1): 014108 (2005), doi:10.1063/1.1824881

$^2$ Barone et al. Int J Quantum Chem 112(9): 2185 (2012), doi:10.1002/qua.23224

$\endgroup$

• Thanks for the insights. However, when I read the book by Frank Jensen, he describes the approach by mentioning that we need to solve the nuclear schrodinger equation. He did mention the diagonalisation of the mass-weighted Hessian. But because it is an SE of the molecule, it contains the KE operator as well and a the end, the approach involves solving 3N 1-dimensional SEs and then obtains the eigenvalues, which are related to the vibrational frequencies by v = 1/2pi * root(eignenvalue). – Tan Yong Boon Jan 1 at 1:10
• But in your approach, there is no nuclear schrodinger equation mentioned anywhere. Is it because his approach is based on QM while yours is classical? – Tan Yong Boon Jan 1 at 1:11
• @TanYongBoon The nuclear vibrational SE is that stated by Jan Jensen in his answer, and both his answer and mine are based on it. Thus, both Jan's answer and mine do provide quantum solutions. – hBy2Py Jan 1 at 15:09
• You are incorrect: this nuclear vibrational SE does not contain a kinetic energy term, as it is usually assumed that the nuclear $\frac{1}{2}mv^2$ quantities are small enough to be neglected and thus only the Hooke's Law potential energy term is present. – hBy2Py Jan 1 at 15:12
• I'm not sure how to convey this to you. But perhaps, if you have the book by Frank Jensen, you can look at p. 312, under Ch. 13, Change of Coordinate System. – Tan Yong Boon Jan 1 at 15:15