Vibrational Analysis and Thermochemistry

This document describes the theoretical background and implementation details for computing vibrational frequencies, IR intensities, and thermochemistry from the Hessian matrix in Metalquicha.

From Hessian to Vibrational Frequencies

The Cartesian Hessian

The Hessian matrix contains second derivatives of the energy with respect to atomic displacements:

\[H_{ij} = \frac{\partial^2 E}{\partial x_i \partial x_j}\]

For a system with N atoms, the Hessian is a 3N x 3N matrix in Cartesian coordinates, with units of Hartree/Bohr².

Mass-Weighting

The Hessian must be mass-weighted before diagonalization. Each element is divided by the square root of the product of the corresponding atomic masses:

\[\tilde{H}_{ij} = \frac{H_{ij}}{\sqrt{m_i \times m_j}}\]

where \(m_i\) is the mass of the atom associated with coordinate i (in atomic mass units converted to atomic units via AMU_TO_AU = 1822.888).

Projection of Translational and Rotational Modes

For an isolated molecule, 6 degrees of freedom (5 for linear molecules) correspond to overall translation and rotation. These must be projected out to obtain pure vibrational frequencies.

Translation Vectors (3 modes)

\[T_x = (1,0,0, 1,0,0, \ldots) / \sqrt{N}\]

for all atoms, and similarly for \(T_y\) and \(T_z\).

Mass-weighted: \(\tilde{T} = T \times \sqrt{m}\)

Rotation Vectors (3 modes for nonlinear, 2 for linear)

For rotation about axis \(\alpha\), the displacement of atom i at position \(r_i\) from center of mass:

\[R_{\alpha,i} = (r_i - r_{COM}) \times \hat{e}_\alpha\]

where \(\hat{e}_\alpha\) is the unit vector along axis \(\alpha\).

Projection Procedure

Orthonormalize the translation/rotation vectors using Gram-Schmidt
Construct projection matrix: \(P = I - \sum_k |v_k\rangle\langle v_k|\)
Apply to mass-weighted Hessian: \(\tilde{H}_{proj} = P \times \tilde{H} \times P\)

Diagonalization

Diagonalize the projected mass-weighted Hessian to obtain eigenvalues \(\lambda_i\) and eigenvectors (normal modes):

\[\tilde{H}_{proj} \times L = L \times \Lambda\]

where \(\Lambda\) is diagonal with eigenvalues \(\lambda_i\).

Conversion to Frequencies

Eigenvalues have units of Hartree/(Bohr² x m_e). Convert to cm⁻¹:

\[\nu_i \text{ (cm}^{-1}\text{)} = \text{sign}(\lambda_i) \times \sqrt{|\lambda_i|} \times \text{AU\_TO\_CM1}\]

where AU_TO_CM1 = 2.642461 x 10^7 is the conversion factor.

Negative eigenvalues indicate imaginary frequencies (saddle points), reported as negative cm⁻¹ values.

Additional Vibrational Properties

Reduced Masses

\[\mu_i = \frac{1}{\sum_j (L_{ij})^2}\]

where the sum is over mass-weighted normal mode components. Units: amu.

Force Constants

\[k_i = 4\pi^2 c^2 \nu_i^2 \mu_i\]

Can be converted to mDyne/Å using AU_TO_MDYNE_ANG = 15.569141.

Cartesian Displacements

Un-mass-weight the eigenvectors:

\[d_{ij} = \frac{L_{ij}}{\sqrt{m_j}}\]

IR Intensities from Dipole Derivatives

Theory

IR intensity is proportional to the square of the transition dipole moment, which depends on how the dipole moment changes along the normal mode coordinate:

\[I_i \propto \left|\frac{\partial \mu}{\partial Q_i}\right|^2\]

Dipole Derivative Matrix

The dipole derivative matrix (3 x 3N) contains derivatives of the dipole moment with respect to Cartesian displacements:

\[\left(\frac{\partial \mu}{\partial x}\right)_{\alpha i} = \frac{\partial \mu_\alpha}{\partial x_i}\]

where \(\alpha \in \{x, y, z\}\) and i runs over all 3N Cartesian coordinates.

This is computed via finite differences of the dipole moment during the Hessian calculation.

Transformation to Normal Coordinates

Transform dipole derivatives from Cartesian to normal mode coordinates:

\[\frac{\partial \mu}{\partial Q_i} = \sum_j \frac{\partial \mu}{\partial x_j} \times \frac{\partial x_j}{\partial Q_i}\]

The transformation matrix \(\partial x / \partial Q\) comes from the mass-weighted normal modes:

\[\frac{\partial x_j}{\partial Q_i} = \frac{L_{ji}}{\sqrt{m_j}}\]

IR Intensity Calculation

For each normal mode i:

\[I_i = \frac{N_A \pi}{3 \times 4\pi\varepsilon_0 \times c^2} \times \left|\frac{\partial \mu}{\partial Q_i}\right|^2\]

In atomic units, the conversion factor to km/mol is:

AU_TO_KMMOL = 1.7770969e6

The intensity is:

\[I_i \text{ (km/mol)} = \text{AU\_TO\_KMMOL} \times \sum_\alpha \left(\frac{\partial \mu_\alpha}{\partial Q_i}\right)^2\]

Thermochemistry (RRHO Approximation)

The Rigid Rotor Harmonic Oscillator (RRHO) approximation computes thermodynamic properties by treating:

Translation as an ideal gas
Rotation as a rigid rotor (classical limit)
Vibration as quantum harmonic oscillators

Input Requirements

Molecular geometry (coordinates, atomic numbers)
Vibrational frequencies (cm⁻¹)
Temperature T (default: 298.15 K)
Pressure P (default: 1 atm)
Symmetry number σ (default: 1)
Spin multiplicity (default: 1)

Moments of Inertia

Compute the inertia tensor in the center-of-mass frame:

\[I_{\alpha\beta} = \sum_i m_i \times (r_i^2 \delta_{\alpha\beta} - r_{i\alpha} \times r_{i\beta})\]

Diagonalize to get principal moments \(I_A, I_B, I_C\) (in amu·Å²).

Linear molecule detection: \(I_A \approx 0\) (smallest moment below threshold).

Rotational Temperatures

\[\Theta_{rot} = \frac{h^2}{8\pi^2 I k_B} = \frac{24.2637}{I} \text{ [K, for I in amu·Å²]}\]

Zero-Point Energy (ZPE)

\[\text{ZPE} = \frac{1}{2} h c \sum_i \nu_i = \frac{1}{2} \sum_i \nu_i \times \text{CM1\_TO\_KELVIN} \times k_B\]

where CM1_TO_KELVIN = 1.4387773538277 K/cm⁻¹.

Only positive (real) frequencies are included. Imaginary frequencies are skipped with a warning.

Translational Contributions

For an ideal gas:

\[\begin{split}E_{trans} &= \frac{3}{2} RT \\ H_{trans} &= \frac{5}{2} RT \quad \text{(includes PV = RT)} \\ S_{trans} &= R \left[\frac{5}{2} + \ln(q_{trans})\right] \quad \text{(Sackur-Tetrode)} \\ C_{p,trans} &= \frac{5}{2} R\end{split}\]

Partition function:

\[q_{trans} = \frac{V}{\lambda^3} = \left(\frac{2\pi m k T}{h^2}\right)^{3/2} \times \frac{kT}{P}\]

Rotational Contributions

Nonlinear molecules (3 rotational DOF):

\[\begin{split}E_{rot} &= \frac{3}{2} RT \\ S_{rot} &= R \left[\frac{3}{2} + \ln(q_{rot})\right] \\ C_{v,rot} &= \frac{3}{2} R \\ q_{rot} &= \frac{\sqrt{\pi}}{\sigma} \times \frac{T^{3/2}}{\sqrt{\Theta_A \times \Theta_B \times \Theta_C}}\end{split}\]

Linear molecules (2 rotational DOF):

\[\begin{split}E_{rot} &= RT \\ S_{rot} &= R \left[1 + \ln\left(\frac{T}{\sigma \Theta_{rot}}\right)\right] \\ C_{v,rot} &= R \\ q_{rot} &= \frac{T}{\sigma \Theta_{rot}}\end{split}\]

Vibrational Contributions

For each mode with vibrational temperature \(\theta_v = 1.4388 \times \nu\) (K):

Define \(u = \theta_v / T\):

\[\begin{split}E_{vib} &= R \sum_i \frac{\theta_{v,i}}{e^{u_i} - 1} \quad \text{[thermal only, excludes ZPE]} \\ S_{vib} &= R \sum_i \left[\frac{u_i}{e^{u_i}-1} - \ln(1-e^{-u_i})\right] \\ C_{v,vib} &= R \sum_i \frac{u_i^2 e^{u_i}}{(e^{u_i}-1)^2} \\ q_{vib} &= \prod_i \frac{1}{1 - e^{-u_i}}\end{split}\]

Electronic Contributions

Ground state only:

\[\begin{split}E_{elec} &= 0 \\ S_{elec} &= R \ln(2S+1) \quad \text{[spin multiplicity]}\end{split}\]

Total Thermodynamic Functions

\[\begin{split}\text{Thermal correction to Energy:} \quad E_{corr} &= \text{ZPE} + E_{trans} + E_{rot} + E_{vib} \\ \text{Thermal correction to Enthalpy:} \quad H_{corr} &= E_{corr} + RT \\ \text{Thermal correction to Gibbs:} \quad G_{corr} &= H_{corr} - T \times S_{total}\end{split}\]

\[\begin{split}\text{Total Enthalpy:} \quad H &= E_{elec} + H_{corr} \\ \text{Total Free Energy:} \quad G &= E_{elec} + G_{corr}\end{split}\]

Implementation Notes

Matching Reference Programs

During development, the thermochemistry output was validated against xtb (extended tight-binding). Several adjustments were needed to match reference values:

Rotational Temperature Constant

Issue: Rotational partition function and entropy were off by a factor of ~1.7
Root cause: Incorrect rotational temperature constant
Fix: Changed from 16.8576 / I to 24.2637 / I (K, for I in amu·Å²)

Heat Capacity Convention

Issue: TR (translational) heat capacity showed 2.981 cal/K/mol instead of 4.968
Root cause: Reported Cv instead of Cp
Fix: For the TR row, report \(C_p = C_v + R = \frac{5}{2}R\) instead of \(C_v = \frac{3}{2}R\)

Enthalpy Table (TOT row)

Issue: TOT enthalpy showed ~15000 cal/mol instead of ~2372
Root cause: ZPE was included in the TOT row of the contribution table
Fix: TOT row shows only thermal contributions (VIB + ROT + TR), ZPE is reported separately

Thermal Correction Display

Issue: H(0)-H(T)+PV value was ~6x too large
Root cause: ZPE was included in this quantity
Fix: H(0)-H(T)+PV represents only the thermal correction without ZPE: \(E_{trans} + E_{rot} + E_{vib} + RT\)

Key Constants

Physical Constants Used
Constant	Value	Units	Description
`AMU_TO_AU`	1822.888	m_e/amu	Mass conversion
`AU_TO_CM1`	2.642461 x 10^7	cm⁻¹	Frequency conversion
`CM1_TO_KELVIN`	1.4387773538277	K/cm⁻¹	Vibrational temperature
Θ_rot constant	24.2637	K·amu·Å²	Rotational temperature
`R_CALMOLK`	1.98720425864	cal/(mol·K)	Gas constant
`AU_TO_KMMOL`	1.7770969 x 10^6	km/mol	IR intensity conversion