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}$.




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$


 15

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](https://arxiv.org/abs/1306.5499v1) of George Phillies, the classic [Wilson, Decius & Cross](https://books.google.com/books?id=pfi7AQAAQBAJ), and the method as described on the [Gaussian.com website](http://gaussian.com/vib/). 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](https://en.wikipedia.org/wiki/Normal_mode#Coupled_oscillators) 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](https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process).


By construction, $\mbf D$ is an [orthogonal matrix](https://en.wikipedia.org/wiki/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](https://en.wikipedia.org/wiki/Atomic_units) 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](http://physics.nist.gov/cgi-bin/cuu/Value?meu%7Csearch_for=atomic%20mass).


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](https://en.wikipedia.org/wiki/Quantum_harmonic_oscillator#N-dimensional_harmonic_oscillator)), 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](https://chemistry.stackexchange.com/a/74686/11367), the speed of light [expressed in atomic units](https://en.wikipedia.org/wiki/Atomic_units#Related_physical_constants) 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](https://chemistry.stackexchange.com/a/74686/11367). 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](http://dx.doi.org/10.1063/1.1824881)


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


$\endgroup$


 14