Vibrational analysis#

The following steps are carried out mostly by the geomeTRIC module [WS16], hence they are described in less detail. For more details on the topic, the reader is referred to Molecular Vibrations by Wilson, Decius and Cross [WDC80].

Cartesian and mass-weighted Hessian#

The starting point for the vibrational analysis of molecules is the Hessian matrix in Cartesian coordinates, \(\mathbf{H}^{\text{Cart}}\), calculated either numerically or analytically as described above. Generally, the elements of \(\mathbf{H}^{\text{Cart}}\) are given by second derivatives of the energy \(E\) with respect to nuclear displacement,

\[%:label: eq:Hessian_matrix_elements H_{ij}^{\text{Cart}} = \bigg( \frac{\mathrm{d}^2 E}{\mathrm{d} \xi_i \mathrm{d} \xi_j} \bigg)_0 \, .\]

Hence, \(\mathbf{H}^{\text{Cart}}\) is a \(3N \times 3N\) matrix (where \(N\) is the number of atoms), and \(\xi_1, \xi_2, \ldots, \xi_{3N}\) is used for the displacements of the Cartesian coordiates, \(\Delta x_1, \Delta y_1, \Delta z_1, \ldots, \Delta z_N\). The ‘\(0\)’ subscript of the parentheses refers to the equilibrium geometry of the atoms at which the derivatives are taken, and that the first derivatives vanish.

As a first step, the Hessian is converted to mass-weighted Cartesian coordinates (MWC), \(q_1 = \sqrt{m_1} \xi_1 = \sqrt{m_1} \Delta x_1\), \(q_2 = \sqrt{m_1} \xi_2 = \sqrt{m_1} \Delta y_1\), \(\ldots\), \(q_{3N} = \sqrt{m_N} \xi_{3N} = \sqrt{m_N} \Delta z_{N}\), where \(m_i\) is the mass of atom \(i\), such that \(\mathbf{H}^{\text{MWC}}\) is given by

\[%:label: eq:Hessian_mwc H_{ij}^{\text{MWC}} = \frac{H_{ij}^{\text{Cart}}}{\sqrt{m_i m_j}} = \bigg( \frac{\mathrm{d}^2 E}{\mathrm{d} q_i \mathrm{d} q_j} \bigg)_0 \, .\]

Diagonalizing this Hessian gives \(3N\) eigenvalues are the fundamental frequencies of the molecule, which still include the translation and rotational modes. However, these should be close to zero.

Translating and rotating frame#

In order to remove translational and rotational degrees of freedom, one first determines the center of mass (COM) \(\mathbf{R}^{\text{COM}}\) in the usual way,

\[%:label: eq:center_of_mass \mathbf{R}^{\text{COM}} = \frac{\sum_{K} m_{K} \mathbf{R}_K}{\sum_{K} \mathbf{R}_K} \, ,\]

where the sum runs over all atoms \(K\), and the origin is then shifted to the COM, \(\mathbf{R}_{K}^{\text{COM}} = \mathbf{R}_K - \mathbf{R}^{\text{COM}}\). Subsequently, one determines the inertia tensor and diagonalizes it to obtain principal moments and axes of inertia. Next, one needs to find the transformation from mass-weighted Cartesian coordinates to a set of \(3N\) coordinates, where the molecule’s translation and rotation are separated out, leaving \(3N - 6\) (or \(3N-5\) for linear molecules) vibrational modes.

This is achieved by applying the so-called Eckart conditions [Eck35]. While the three vectors of length \(3N\) corresponding to translation are simply given by \(\sqrt{m_i}\) times the coordinate axis, the vectors corresponding to rotational motion of the atoms are obtained from the coordinates of the atoms with respect to the COM and the corresponding row of the matrix used to diagonalize the moment of inertia tensor. In the next step, these vectors are normalized and a Gram–Schmidt orthogonalization is carried out to create \(N_\text{vib} = 3N-6\) (or \(3N-5\)) remaining vectors, which are orthogonal to the five or six translational and rotational vectors. Thus, one obtains a transformation matrix \(\mathbf{D}\) which allows for the transformation of the mass-weighted Cartesian coordinates \(\mathbf{q}\) to internal coordinates \(\mathbf{S} = \mathbf{Dq}\), where translation and rotation have been projected out.

Hessian in internal coordinates and harmonic frequencies#

Now the Hessian \(\mathbf{H}^\text{MWC}\), which is still given in mass-weighted Cartesian coordinates, is transformed the the internal coordinate system,

\[%:label: eq:Hessian_to_internal \mathbf{H}^{\text{Int}} = \mathbf{D}^\dagger \mathbf{H}^\text{MWC} \mathbf{D} \, ,\]

yielding a representation in \(N_\text{vib}\) internal coordinates from the full \(3N\) Cartesian coordinates. The Hessian in internal coordinates \(\mathbf{H}^{\text{Int}}\) is successively diagonalized,

\[%:label: eq:Hessian_internal_diagonalized \mathbf{L}^\dagger \mathbf{H}^{\text{Int}} \mathbf{L} = \mathbf{\Lambda} \, ,\]

where \(\mathbf{\Lambda}\) is the diagonal matrix of \(N_{\text{vib}}\) eigenvalues \(\lambda_i\) which are related to the harmonic vibrational frequencies \(\nu_i\) and \(\mathbf{L}\) is the transformation matrix composed of the eigenvectors.

Finally, the eigenvalues \(\lambda_i = 4 \pi^2 \nu_i^2\) can be converted from frequencies \(\nu_i\) to wavenumbers \(\tilde{\nu}_i\) in reciprocal centimeters by using the relationship \(\nu_i = c \tilde{\nu}_i\), where \(c\) is the speed of light. The wavenumbers are thus obtained from

\[%:label: eq:vibrations_wavenumbers \tilde{\nu}_i = \sqrt{\frac{\lambda_i}{4\pi^2 c^2}} \, ,\]

and successively appropriate conversion factors are applied to obtain the wavenumbers in inverse centimeters (cm\(^{-\text{1}}\)).

Cartesian displacements, reduced masses, and force constants#

The Cartesian normal modes \(\mathbf{l}^{\text{Cart}}\) are obtained by combining this and this equation together with a diagonal matrix \(\mathbf{M}\) defined by \(M_{ii} = \frac{1}{\sqrt{m_i}}\) to undo the mass-weighting, \(\mathbf{l}^{\text{Cart}} = \mathbf{M D L}\), with the individual elements of this matrix being given by

\[ l_{ij}^{\text{Cart}} = \sum_{k=1}^{3N} \frac{D_{ik} L_{kj}}{\sqrt{m_i}} \, .\]

The (normalized) column vectors of \(\mathbf{l}^{\text{Cart}}\) correspond to the normal-mode displacements in Cartesian coordinates, which are used together with property gradients for the calculation of spectral intensities as described here.

From the Cartesian normal modes \(\mathbf{l}^{\text{Cart}}\), the reduced mass \(\mu_i\) of vibration \(i\) can be calculated as

\[%:label: eq:reduced_masses \mu_i = \frac{1}{\sum_{k=1}^{3N} \big( l_{ki}^{\text{Cart}} \big)^2} \, ,\]

and from those the corresponding force constants \(k_i\) are calculated as

\[%:label: eq:force_constants k_i = 4 \pi^2 \tilde{\nu}_i^2 \mu_i \, ,\]

since \(\tilde{\nu}_i = \frac{1}{2 \pi} \sqrt{\frac{k_i}{\mu_i}}\). The force constants are then converted from atomic units to milli-dyne per ångström.

The spectral intensities are calculated from the dipole moment gradient (IR spectroscopy), or the polarizability gradient (Raman spectroscopy).