Chemistry Stack Exchange is a question and answer site for scientists, academics, teachers, and students in the field of chemistry. It only takes a minute to sign up.
Sign up to join this community
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$
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$
