In this part, we will solve the Schrödinger equation for molecular vibrations using the variational method (see Variational method). We will use the vibrational Hamiltonian, composed of the kinetic energy operator and potential energy surface, calculated in previous sections. We will learn how to create and optimize multivariate basis sets, as well as how to utilize Gaussian quadratures effectively.
The kinetic and potential energy operators for water molecule ($\text{H}_2\text{O}$) are expressed in a sum-of-products form in terms of valence-bond coordinates - two O—H distances $r_1$ and $r_2$, and H—O—H angle $\alpha$. The total Hamiltonian is a sum of kinetic and potential energy operators
$$ \hat{H} = \hat{T} + V(r_1,r_2,\alpha) = -\frac{\hbar^2}{2}\sum_{\lambda,\mu=1..3}\frac{\partial}{\partial q_\lambda}G_{\lambda\mu}(r_1,r_2,\alpha)\frac{\partial}{\partial q_\mu}+V(r_1,r_2,\alpha), $$
where
$$ G_{\lambda\mu} = \sum_{pqr} g_{pqr}^{(\lambda\mu)} (r_1-r_1^\text{eq})^p (r_2-r_2^\text{eq})^q (\alpha-\alpha^\text{eq})^r, \\ V = \sum_{pqr} f_{pqr} (1-e^{-a_1(r_1-r_1^\text{eq})})^p (1-e^{-a_2(r_2-r_2^\text{eq})})^q (\cos\alpha-\cos\alpha^\text{eq})^r, $$
and $q\in\{r_1, r_2,\alpha\}$. The kinetic energy $G$ matrix expansion coefficients $g_{pqr}^{(\lambda\mu)}$ may be obtained using forward-mode automatic differentiation, as demonstrated in section Forward-mode automatic differentiation. The potential energy $V$ expansion coefficients $f_{pqr}$ are obtained by fitting to the electronic energies computed on a grid of different molecular geometries using one of the quantum chemical methods.
We will be seeking the solution of the Schrodinger equation $\hat{H}\Psi_l = E_l\Psi_l$ by expanding the wave function $\Psi_l$ into linear combination of products of one-dimensional functions of chosen vibrational coordinates, i.e.,
$$ \Psi_l(r_1,r_2,\alpha)=\sum_{ijk}C_{ijk}^{(l)}\phi_i(r_1)\chi_j(r_2)\psi_k(\alpha). $$
The variational principle yields the following eigenvalue equation for determining the unknown coefficients $C_{ijk}^{(l)}$ (eigenvalues) and energies $E_l$ (eigenvectors):
$$ HC=SCE, \\ H_{p,q} \equiv H_{i'j'k',ijk} = \langle \phi_{i'}\chi_{j'}\psi_{k'}| \hat{H}|\phi_i\chi_j\psi_k \rangle, \\ S_{p,q} \equiv S_{i'j'k',ijk} = \langle \phi_{i'}\chi_{j'}\psi_{k'}|\phi_i\chi_j\psi_k \rangle = \langle \phi_{i'}|\phi_i \rangle \langle \chi_{j'}|\chi_j \rangle \langle \psi_{k'}|\psi_k \rangle .
$$
The number of terms in the expansion of total wave function is determined by the quality of the basis functions. Optimized functions result in a lesser number of terms, whereas non-optimized functions result in a higher number of expansion terms.
In the spirit of the perturbation theory, the whole Hamiltonian may be divided into two parts: the first largest zero-order Hamiltonian $H_0$, for which the solution of the eigenvalue problem is known, and the second smallest $\hat{H}'$, which is viewed as perturbation, i.e.
$$
\hat{H} = \hat{H}_0 + \hat{H}'~\text{and}~\hat{H}_0 \gg \hat{H}'.
$$
The wave function of the total Hamiltonian $\hat{H}$ may be expected to converge quickly if expanded in terms of eigenfunctions of $\hat{H}_0$. The same perturbation method may be applied to the zero-order Hamiltonian $\hat{H}_0$ itself by dividing it into a simpler zero-order component and a perturbation. In principle, this approach may be continued until the zero-order Hamiltonian with analytical eigenvalue solutions is reached, such as, for example, Harmonic oscillator, or particle in a box. Applied our water molecule, this may look like the following:
Three-step basis set contraction (optimization) for triatomic molecule
Thus, instead of solving a huge eigenvalue problem for the whole Hamiltonian using a basis of primitive functions (e.g., harmonic-oscillator functions), this technique solves multiple eigenvalue problems of smaller sizes, at each step making the basis more compact, or contracted. Given that the solution of the eigenvalue problem scales as the third power of the number of basis functions, basis set contraction often results in increased computational efficiency even for small molecules.