Introduction
In this post we will report on our efforts to determine the normal modes of atomic vibrations in a water molecule and compute the associated frequencies. For this we will explicitly calculate the mass-weighted Hessian matrix elements by calculating the change in energy when atoms are displaced by finite distances and thus numerically calculating the second derivatives of the energy of the water molecule with respect to the atomic coordinates (i.e. the Hessian matrix elements). We will compare our results for vibrational frequencies of normal modes in H\(_2\)O with the normal mode vibrational frequencies calculated by Materials Studio.
All energy calculations were performed using DFT calculations using the DMol3 software package as implemented by Materials Studio. DMol3 employs numerical radial functions centered at atomic positions as basis sets. The GGA base PBE functional was used to treat exchange correlational effects and the double numerical precision – DNP 3.5 basis set was used. All energy calculations were deemed to be converged once an SCF convergence of \(10^{-6}\) was achieved.
In the Harmonic approximation, expanding the energy of the water molecule about the equilibrium positions is,
\begin{equation}E=E_0+\frac{1}{2}\sum_{i=1}^{9}\sum_{j=1}^{9}\left[\frac{\partial^2E}{\partial x_i\partial x_j}\right]_{\{x_i\}_{eq}}x_ix_j.\label{HarmonicApproximation}\end{equation}
Where, \(x_i\) is the \(\left(i-3\lfloor\frac{i-1}{3}\rfloor\right)\)-th coordinate of the \(\lceil\frac{i}{3}\rceil\)-th atom. Hence the sum over \(i\) runs to 9 for the 3 possible coordinates of the 3 atoms in H\(_2\)O.
The Hessian matrix elements are defined as,
\begin{equation}H_{ij}=\left[\frac{\partial^2E}{\partial x_i\partial x_j}\right]_{\{x_i\}_{eq}}.\label{Hessian}\end{equation}
The matrix form of the the equations of motion for the atomic vibrations can be written as, (according to equation 5.7)
\begin{equation}\frac{d^2\boldsymbol{x}}{dt^2}=-\boldsymbol{Ax}.\label{EOM}\end{equation}
Where \(\boldsymbol{A}\) is the mass-weighted Hessian matrix whose matrix elements are given by,
\begin{equation}A_{ij}=\frac{1}{m_i}H_{ij}.\label{MassWeightedHessian}\end{equation}
Once A is determined the normal modes and corresponding vibrational frequencies can be determined using the eigenvectors {e\(_i\)} and eigenvalues {\(\lambda_i\)} respectively. The i-th eigenvector e\(_i\) define thes i-th normal mode of atomic vibrations and the corresponding eigenvalue (\(\lambda_i\)) defines the angular frequency of that vibrational mode by, \(\omega_i=\sqrt{\lambda_i}\).
Calculations
Figure 1: Relaxed H\(_2\)O molecule. The molecule has been oriented such that the O atom is at the origin of the coordinate system and the entire unperturbed molecule is in the xy-plane and symmetric about the yz-plane.
The notation \(\delta x_i=\delta r\) is used to denote an atomic displacement. Here \(\delta r\) gives the magnitude of the displacement, \(\alpha=\lceil \frac{i}{3} \rceil \) specifies the atom that is moving (\( \alpha \in \{1,2,3\}\) where the atoms are numbered as shown in figure 1) and \(\left(i-3\lfloor\frac{i-1}{3}\rfloor\right)\) denotes the Cartesian coordinate of the displacement. For example \(\delta x_6=0.1\mathring{A}\) means the H atom assigned #2 has been moved \(0.1\;\mathring{A}\) in the z direction, from its equilibrium position.
The Hessian matrix elements were calculated using numerical differentiation for a finite displacements of the atomic positions. In principle the numerical differentiation approaches the actual derivative for infinitesimally small displacements, however we need to consider the accuracy of our energy calculations and choose a displacement that is sufficiently large to produce a difference in energy between the unperturbed and perturbed molecule.
Table 1: Energy of a perturbed H\(_2\)O molecule and the corresponding change in energy from the unperturbed energy \(E_0\) for different \(\delta x_1\) values
| \(\delta x_1\) (Å) | \(E(\delta x_1) \) (Ha) | \(E(\delta x_1)\) -\(E_0\) (Ha) |
| 0.02 | -76.378143 | 0.000462 |
| 0.03 | -76.377565 | 0.00104 |
| 0.04 | -76.376754 | 0.001851 |
| 0.05 | -76.37571 | 0.002895 |
| 0.06 | -76.374432 | 0.004173 |
Table 1 above lists the energy of a perturbed H\(_2\)O molecule for different \(\delta x_1\) values and the corresponding change in energy from the unperturbed energy \(E_0 = -76.378605 \)Ha. From the values in the table we decided to use 0.04 Å displacements to perturb the molecule in order to obtain a perturbation in energy that is \(>\;10^{-3}\;\mathring{A}\). For the rest of the post an atomic displacement will be specified by \(\pm\delta x_i\), since the magnitude of any displacement is always chosen to be 0.04 Å.
When evaluating the derivatives using numerical differentiation we chose to use the symmetric difference quotient (SDQ) method for improved accuracy. In this method the derivative of a function \(f(x)\) at \(x=x_0\) is calculated, (using a small displacement from \(x_0\) , h)
\begin{equation}f^{\prime}=\frac{f(x_0+h)-f(x_0-h)}{2h}.\label{SDQ}\end{equation}
The alternative methods to compute the derivatives are the forward difference quotient (FDQ) and the backward difference quotient (BDQ), which are defined in as,
\begin{equation}f^{\prime}=\pm\frac{f(x_0\pm h)-f(x_0)}{h}.\label{FDQBDQ}\end{equation}.
Where choosing the + ve sign gives the FDQ and – sign gives the BDQ.
Now to show that SDQ is the best choice; shown below in table 2 is the gradient of the total energy with respect to \(\delta \boldsymbol{x}_3\) evaluated at the equilibrium position using the SDQ, FDQ and BDQ. Since we are at equilibrium the actual gradient with respect to any coordinate should vanish by definition. None of our numerical results gradients vanish exactly but the gradients obtained using SDQ are consistently lower than the gradients from the other methods by atleast an order of magnitude, indicating this is the best choice to compute derivatives.
The gradient of the total energy with respect to \(\delta x_7\) evaluated at the equilibrium position using the SDQ, FDQ and BDQ.
| \(\delta x_7\) (Å) | \(E(+\delta x_7)\) (Ha) | \(E(-\delta x_7)\) (Ha) | \(\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{dE}{dx_7}\) \(\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\)(Ha/Å) |
||
| FDQ | BDQ | SDQ | |||
| 0.02 | -76.378353 | -76.378351 | 0.0126 | -0.0127 | -5E-05 |
| 0.03 | -76.37805 | -76.378023 | 0.0185 | -0.0194 | -0.00045 |
| 0.04 | -76.377634 | -76.377553 | 0.024275 | -0.0263 | -0.0010125 |
| 0.05 | -76.377111 | -76.376164 | 0.02988 | -0.04882 | -0.00947 |
| 0.06 | -76.376485 | -76.376935 | 0.035333333 | -0.027833333 | 0.00375 |
With our choice to use the SDQ method to compute numerical derivatives, the Hessian matrix elements defined in equation \eqref{Hessian} become,
\begin{equation}H_{ij}=\left[\frac{\partial^2E}{\partial x_i\partial x_j}\right]_{\{x_i\}_{eq}}=\frac{E(\delta x_i +\delta x_j)+E(-\delta x_i -\delta x_j)-E(\delta x_i -\delta x_j)-E(-\delta x_i +\delta x_j)}{4 \delta x_i \delta x_j}.\label{Hessian_num}\end{equation}
Results
Table 3 below shows the Hessian matrix calculated using equation \eqref{Hessian_num}. The units of the matrix elements are \(Ha/\mathring{A}^2\). And the corresponding mass-weighted Hessian matrix is shown in table 4, with matrix elements having units of \(\frac{Ha}{amu.\mathring{A}^2}\).
Table 3. The Hessian of the water molecule, matrix elements have units of \(Ha/\mathring{A}^2\)
| \(\delta x_1\) | \(\delta x_2\) | \(\delta x_3\) | \(\delta x_4\) | \(\delta x_5\) | \(\delta x_6\) | \(\delta x_7\) | \(\delta x_8\) | \(\delta x_9\) | |
| \(\delta x_1\) | 2.321875 | 2.22E-12 | 0 | -1.16109375 | 0.9196875 | 0.00E+00 | -1.16109375 | -0.9196875 | 0 |
| \(\delta x_2\) | 2.22E-12 | 1.55828125 | 0 | 0.7003125 | -0.77921875 | 0 | -0.7003125 | -0.77921875 | 0 |
| \(\delta x_3\) | 0 | 0 | 0.0078125 | -2.22E-12 | 2.22E-12 | -0.0034375 | -2.22E-12 | 2.22E-12 | -0.0034375 |
| \(\delta x_4\) | -1.16109375 | 0.7003125 | -2.22E-12 | 1.26921875 | -0.8096875 | 0 | -0.10796875 | 0.11 | 2.22E-12 |
| \(\delta x_5\) | 0.9196875 | -0.77921875 | 2.22E-12 | -0.8096875 | 0.74421875 | 2.22E-12 | -0.11 | 0.03546875 | 0 |
| \(\delta x_6\) | 0 | 0 | -0.0034375 | 0 | 2.22E-12 | 0.0040625 | 2.22E-12 | 0 | -0.0009375 |
| \(\delta x_7\) | -1.16109375 | -0.7003125 | -2.22E-12 | -0.10796875 | -0.11 | 2.22E-12 | 1.26921875 | 0.8096875 | 0 |
| \(\delta x_8\) | -0.9196875 | -0.77921875 | 2.22E-12 | 0.11 | 0.03546875 | 0 | 0.8096875 | 0.74421875 | 2.22E-12 |
| \(\delta x_9\) | 0 | 0 | -0.0034375 | 2.22E-12 | 0 | -0.0009375 | 0 | 2.22E-12 | 0.0040625 |
Table 04 : The mass-weighted Hessian of the water molecule. [Units : \(\left(\frac{Ha}{amu.\mathring{A}^2}\right)\) ]
| \(\delta x_1\) | \(\delta x_2\) | \(\delta x_3\) | \(\delta x_4\) | \(\delta x_5\) | \(\delta x_6\) | \(\delta x_7\) | \(\delta x_8\) | \(\delta x_9\) | |
| \(\delta x_1\) | 0.14512263 | 1.39E-13 | 0 | -0.072571081 | 0.057482624 | 0.00E+00 | -0.072571081 | -0.057482624 | 0 |
| \(\delta x_2\) | 1.39E-13 | 0.09739623 | 0 | 0.043771173 | -0.048702998 | 0 | -0.043771173 | -0.048702998 | 0 |
| \(\delta x_3\) | 0 | 0 | 0.0004883 | -1.39E-13 | 1.39E-13 | -0.000214852 | -1.39E-13 | 1.39E-13 | -0.000214852 |
| \(\delta x_4\) | -1.151993005 | 0.694823395 | -2.20E-12 | 1.259270513 | -0.803341105 | 0 | -0.107122482 | 0.109137811 | 2.20E-12 |
| \(\delta x_5\) | 0.912478917 | -0.773111172 | 2.20E-12 | -0.803341105 | 0.738385505 | 2.20E-12 | -0.109137811 | 0.035190743 | 0 |
| \(\delta x_6\) | 0 | 0 | -0.003410557 | 0 | 2.20E-12 | 0.004030658 | 2.20E-12 | 0 | -0.000930152 |
| \(\delta x_7\) | -1.151993005 | -0.694823395 | -2.20E-12 | -0.107122482 | -0.109137811 | 2.20E-12 | 1.259270513 | 0.803341105 | 0 |
| \(\delta x_8\) | -0.912478917 | -0.773111172 | 2.20E-12 | 0.109137811 | 0.035190743 | 0 | 0.803341105 | 0.738385505 | 2.20E-12 |
| \(\delta x_9\) | 0 | 0 | -0.003410557 | 2.20E-12 | 0 | -0.000930152 | 0 | 2.20E-12 | 0.004030658 |
Now that we have the mass-weighted Hessian matrix, we calculate the eigenvalues and eigenvectors to determine the vibrational frequencies and corresponding normal modes. The eigenvalues are listed in table 5. The first three eigenvalues are non-negligible and are the vibrational frequencies of the normal modes of vibrations in a water molecule. The eigenvectors corresponding to these vibrational frequencies, which define the normal modes of these vibrational frequencies are listed in table 6. Inspecting the results table 6 we see that the normal modes given by \(\boldsymbol{e}_1\), \(\boldsymbol{e}_2\) and \(\boldsymbol{e}_3\) respectively depict the anti-symmetric stretching, symmetric stretching and bending normal modes shown respectively in the animated figures 2, 3 and 4 (click on figures for animated view).
Figure 2 : The Anti-symmetric stretching normal mode
Figure 3 : The Symmetric stretching normal mode
Figure 2 : The Bending normal mode
Table 05: Eigenvalues of the mass weighted Hessian matrix
| i | \(\lambda_i\) \((10^{29} \;s^{-2})\) |
|---|---|
| 1 | 5.285202922 |
| 2 | 4.977711427 |
| 3 | 0.896378227 |
| 4 | -0.032981179 |
| 5 | 0.013024606 |
| 6 | 0.009386832 |
| 7 | 2.21E-07 |
| 8 | 0.000113396 |
| 9 | 3.56E-05 |
Table 06 : The normal modes corresponding to the non negligible eigenvalues of the mass-weighted Hessian matrix (i.e. the normal modes of vibration in H\(_2\)O)
| \(x_i\) | \(\boldsymbol{e}_1\) | \(\boldsymbol{e}_2\) | \(\boldsymbol{e}_3\) |
|---|---|---|---|
| x\(_1\) | -0.069920136 | -1.19E-12 | -4.15E-13 |
| x\(_2\) | 7.82E-13 | -0.049910964 | -0.070473846 |
| x\(_3\) | -7.67E-14 | 5.79E-14 | 4.52E-13 |
| x\(_4\) | 0.555004088 | -0.584720384 | 0.427667384 |
| x\(_5\) | -0.435345897 | 0.396051158 | 0.560907592 |
| x\(_6\) | 1.31E-13 | 1.14E-12 | 8.48E-13 |
| x\(_7\) | 5.55E-01 | 0.584720384 | -0.427667384 |
| x\(_8\) | 0.435345897 | 0.396051158 | 0.560907592 |
| x\(_9\) | 1.09E-12 | -2.20E-13 | 6.45E-12 |
Now we compute the vibrational frequencies of the normal modes and compare them with the vibrational frequencies calculated by Material Studio in table 7.
Table 07 : The vibrational frequencies (reported in inverse wavelength) of normal modes calculated in this work compared with the results from Materials Studio.
| Normal Mode | Vibrational Frequency (cm\(^{-1}\)) | Discrepancy between results \(\frac{\left|This\;work-Materials\;Studio\right|}{This\;work}\times100\) (%) |
|
| This Work | Materials Studio Results | ||
| Anti-symmetric stretching | 3859.40 | 3848.97 | 0.27 |
| Symmetric stretching | 3745.44 | 3734.41 | 0.29 |
| Bending | 1589.40 | 1615.69 | 1.65 |
As we can see from the results in table 7 our explicit calculation is in very good agreement with the results from Materials Studio. This is not surprising since Materials Studio employs a similar method to evaluate the Hessian matrix. In Materials Studio the Hessian matrix elements are obtained by displacing each atom in the system and computing a gradient vector, which is then used to build the second derivative matrix.
References
