Iterative Methods for Computing Vibrational Spectra

Abstract

I review some computational methods for calculating vibrational spectra. They all use iterative eigensolvers to compute eigenvalues of a Hamiltonian matrix by evaluating matrix-vector products (MVPs). A direct-product basis can be used for molecules with five or fewer atoms. This is done by exploiting the structure of the basis and the structure of a direct product quadrature grid. I outline three methods that can be used for molecules with more than five atoms. The first uses contracted basis functions and an intermediate (F) matrix. The second uses Smolyak quadrature and a pruned basis. The third uses a tensor rank reduction scheme.

1. Introduction

Effective numerical methods for solving the time-independent Schroedinger equation to compute vibrational spectra of polyatomic molecules have been developed in the last thirty years [1,2,3,4,5]. They are important when approximations, often based on perturbation theory, are not accurate enough. Almost all methods begin by choosing a basis in which to represent both the wavefunctions and the Hamiltonian and then solve a linear algebra problem. These two basic tasks are not independent: a basis with structure favours iterative linear algebra methods (vide infra). Computing vibrational spectra is useful because it helps experimentalists to assign and interpret measured spectra.

In this article, I present a subjective review of several methods for solving the time-independent Schroedinger equation to calculate vibrational spectra [4,6,7]. It is possible to generalize the methods I describe so that they can also be used to compute ro-vibrational spectra [8,9,10,11,12,13,14]. All of the methods presented here obtain solutions to the Schroedinger equation from a space built by evaluating matrix-vector products (MVPs) and are called iterative methods [15]. I shall ignore Multimode-type methods (MM) [5,16,17,18,19,20], that work when the potential energy surface (PES) is a sum of terms that depend on a subset of the coordinates [17,21] (denoted an MM representation) and work quite well for semi-rigid molecules for which normal coordinates are appropriate. Although widely used, I shall also ignore multiconfiguration time-dependent Hartree (MCTDH) methods [22,23]. They can be used with a block power method [24], an “improved relaxation” method, [25,26,27], or a block Lanczos method [28] to compute accurate vibrational energy levels. However, improved relaxation, the most popular MCTDH approach for calculating spectra, converges poorly if the density of states is high and therefore cannot be used to compute a large number of levels of a large molecule [29].

When using iterative methods, it is better not to calculate a Hamiltonian matrix. Many calculations are done with a basis so large that it would not be possible to store the Hamiltonian matrix in memory. Iterative methods require the evaluation of MVPs. How is it possible to compute MVPs without building a matrix representing the Hamiltonian? I shall first outline ideas that make it possible to use a product basis to evaluate MVPs without a Hamiltonian matrix. They exploit the structure of the basis, the quadrature grid, and the kinetic energy operator (KEO). For molecules with more than five atoms, vectors representing wavefunctions and the vector representing the PES on the quadrature grid [4,30,31] are so large that they require too much memory. In the rest of this chapter, I therefore review methods that obviate the need to store large vectors. The first method (Section 4) uses a contracted basis. To make the contracted basis method useful, it is essential that it be possible to evaluate MVPs in the contracted basis without transforming to a huge product grid. The second method (Section 5) uses a pruned basis and a pruned grid. Pruning significantly reduces the size of the largest vectors one must store. The third method (Section 6) builds a basis from MVPs by using tensor rank reduction.

2. Direct Product Basis Sets

When there are D vibrational coordinates, a direct product basis function is

$${\mathsf{\Phi }}_{n_1,n_2,\dots ,n_D}={\varphi }_{n_1}\left(q_1\right){\varphi }_{n_2}\left(q_2\right)\dots {\varphi }_{n_D}\left(q_D\right),$$

(1)

where the indices $\left\{n_k\right\}$ are independent and $n_c=0,1,\dots ,n_c^{max}$. ${\varphi }_{n_c}\left(q_c\right)$ is a 1D basis function for coordinate c. If $n_c^{max}=n\forall c$, then the direct product basis set has $n^D$ functions. The univariate functions are often ${\varphi }_k\left(x\right)=h_k^{-1/2}{\left[w\left(x\right)\right]}^{1/2}{p}_k\left(z\right)$, where z is a function of x, $p_k\left(z\right)$ is a classical orthogonal polynomial, $w\left(x\right)$ is the corresponding weight function, and $h_k$ is a normalization factor. Such a basis is usually called a variational basis representation (VBR) [4,30].

Although there are problems for which a VBR basis is best, it is sometimes advantageous to use a discrete variable representation (DVR) basis [4,30,31,32]. In 1D, a standard DVR basis is a set of orthogonal but localized functions that spans the same space as a set of orthogonal de-localized functions, ${\varphi }_k\left(x\right)$. The 1D DVR Hamiltonian matrix eigenvalue problem is

$${\mathbf{T}}^{\mathbf{T}}(\mathbf{K}+{\mathbf{V}}^{\mathbf{FBR}})\mathbf{T}\mathbf{U}=\mathbf{U}\mathbf{E},$$

(2)

where $\mathbf{K}$ is an exact kinetic matrix in a basis of ${\varphi }_n\left(q\right)$ (VBR) functions and ${\mathbf{V}}^{\mathbf{FBR}}$ is either a product or a quadrature approximation for the exact potential matrix [4]. One way to obtain the transformation matrix $\mathbf{T}$ is to diagonalize the matrix representing x in the VBR,

$$\mathbf{x}\mathbf{T}=\mathbf{T}\mathbf{X},$$

(3)

where $\mathbf{x}$ is the matrix representing x in the ${\varphi }_n\left(q\right)$ basis and $\mathbf{X}$ is a diagonal matrix whose nonzero values are eigenvalues [33]. Equation (2) can be written

$$({\mathbf{T}}^{\mathbf{T}}\mathbf{K}\mathbf{T}+{\mathbf{V}}^{\mathbf{diag}})\mathbf{U}=\mathbf{U}\mathbf{E},$$

where ${\mathbf{V}}^{\mathbf{diag}}$ is a diagonal matrix whose diagonal elements are values of the potential at the quadrature (DVR) points. A potential optimised DVR (PO-DVR) [34,35] is made from 1D basis functions that are solutions of 1D Schroedinger equations.

3. Using a Direct Product Basis Set to Solve the Schroedinger Equation

Direct product DVR and VBR bases are popular [4,36,37,38,39,40]. Although direct product bases are huge, they can be used by exploiting the structure of the basis to efficiently evaluate the MVPs required to use iterative methods. The Lanczos and filter diagonalisation methods are popular iterative methods for solving the time-independent Schroedinger equation [40,41,42,43,44,45,46,47]. For a molecule with as many as five atoms, they make it possible, even with a direct product basis, to solve the vibrational Schroedinger equation with a general PES. The key ideas have been reviewed several times [4,5,36,37]. They are all based on doing sums sequentially [40,48].

In a direct product DVR, potential MVPs are trivial because the potential matrix is diagonal. When the KEO is a sum of products (SOPs) , with g terms each with D factors,

$$\widehat{K}=\sum _{l=1}^g\prod _{k=1}^D{\widehat{h}}^{(k,l)}\left(q_k\right),$$

(4)

then kinetic MVPs can be efficiently evaluated by doing sums sequentially,

$$\sum _{l=1}^g\sum _{n_1}{h}_{n_1^{\prime },n_1}^{(1,l)}\sum _{n_2}{h}_{n_2^{\prime },n_2}^{(2,l)}\dots \sum _{n_D}{h}_{n_D^{\prime },n_D}^{(f,l)}{u}_{n_1,n_2,\dots ,n_D}=u_{n_1^{\prime },n_2^{\prime },\dots ,n_D^{\prime }}^{\prime },$$

(5)

where $h_{n_k^{\prime },n_k}^{(k,l)}$ is an element of the $n\times n$ matrix representation of the factor ${\widehat{h}}^{(k,l)}\left(q_k\right)$. Matrix elements of the full KEO are never computed.

If there are important singularities in the KEO, then a VBR basis is better than a DVR basis [49]. At an important singularity, the KEO is singular and vibrational wavefunctions have significant amplitude. In general, singularities occur whenever one coordinate takes a limiting value and another is undefined [50]. In a VBR basis, it is possible to evaluate the potential MVP by doing sums sequentially [48,49]. This enables one to avoid calculating potential matrix elements, which would require computing many-dimensional integrals. Consider a 2D example. The matrix–vector product is

$$\sum _{n_1}\sum _{n_2}{V}_{n_1^{\prime }{n}_2^{\prime },n_1{n}_2}{u}_{n_1,n_2}=u_{n_1^{\prime },n_2^{\prime }}^{\prime },$$

(6)

where

$$V_{n_1^{\prime }{n}_2^{\prime },n_1{n}_2}=\int d{q}_{1}d{q}_2{\varphi }_{n_1^{\prime }}\left(q_1\right){\varphi }_{n_2^{\prime }}\left(q_2\right)V(q_1,q_2){\varphi }_{n_1}\left(q_1\right){\varphi }_{n_2}\left(q_2\right).$$

(7)

In terms of $\mathbf{T}$ matrices (see Equation (3)) ,

$$V_{n_1^{\prime }{n}_2^{\prime },n_1{n}_2}\approx \sum _{\alpha }\sum _{\beta }{\left(T\right)}_{n_1^{\prime },\alpha }{\left(T\right)}_{n_2^{\prime },\beta }V({\left(q_1\right)}_{\alpha },{\left(q_2\right)}_{\beta }){(T^{†})}_{\alpha ,n_1}{(T^{†})}_{\beta ,n_2}.$$

(8)

The matrix–vector product can be written,

$$\sum _{n_1}\sum _{n_2}\sum _{\alpha }\sum _{\beta }{\left(T\right)}_{n_1^{\prime },\alpha }{\left(T\right)}_{n_2^{\prime },\beta }V({\left(q_1\right)}_{\alpha },{\left(q_2\right)}_{\beta }){(T^{†})}_{\alpha ,n_1}{(T^{†})}_{\beta ,n_2}{u}_{n_1,n_2}=u_{n_1^{\prime },n_2^{\prime }}^{\prime }$$

(9)

and evaluated by doing sums sequentially,

$$\sum _{\alpha }{\left(T\right)}_{n_1^{\prime },\alpha }\sum _{\beta }{\left(T\right)}_{n_2^{\prime },\beta }V({\left(q_1\right)}_{\alpha },{\left(q_2\right)}_{\beta })\sum _{n_1}{(T^{†})}_{\alpha ,n_1}\sum _{n_2}{(T^{†})}_{\beta ,n_2}{u}_{n_1,n_2}=u_{n_1^{\prime },n_2^{\prime }}^{\prime }.$$

(10)

If there is an important singularity, good basis functions are always nondirect product functions which are products of functions of the coordinate which is undefined and the coordinate which takes a limiting value, with a shared index.

The ideas of this section make it possible to compute vibrational spectra without computing and storing a Hamiltonian matrix. However, for molecules with more than five atoms, the memory cost of storing vectors in a direct product basis is prohibitive. For molecules with more than five atoms, it is necessary to introduce other ideas to reduce the memory cost of calculations.

4. Using a DVR to Make a Contracted Basis

To include information about coupling in the basis functions, it is common to use basis functions that are products of factors that depend on more than one coordinate. I shall call the multi-dimensional factors contracted basis functions. It is important to devise good algorithms for evaluating MVPs in a contracted basis. Contracted bases are necessarily more complicated, i.e., they have less structure, and it is structure that is exploited to evaluate MVPS efficiently. An important advantage that contracted bases have is the reduced spectral range of the contracted-basis Hamiltonian matrix. Reducing the spectral range decreases the number of MVPs required to compute eigenvalues.

For molecules with more than three atoms, it is best to use contracted functions diagonalizing matrices that represent the Hamiltonian with one or more coordinates fixed. The basis functions are direct products of functions of different coordinates or groups of coordinates [6,51,52,53].

Evaluating Matrix-Vector Products without Storing a Vector as Large as the Direct Product DVR

The most obvious way to evaluate MVPs is to transform from the contracted basis to a primitive basis, in which the contracted functions are determined. Computing matrix elements of the potential in the primitive basis requires storing the potential on a large (direct product) grid of points. Because the grid is huge, this is impractical for molecules with more than five atoms. An alternative is to store an intermediate matrix [6,54]. To explain how this is done, consider a ($J=0$) Hamiltonian in polyspherical coordinates [55,56,57]

$$H=T_{ben}(\theta ,r)+T_{str}\left(r\right)+V(\theta ,r)$$

(11)

with

$$\begin{array}{ccc}\hfill T_{ben}(\theta ,r)& =& \sum _i{B}_i\left(r\right)T_b^{\left(i\right)}\left(\theta \right)\hfill \\ \hfill T_{str}\left(r\right)& =& \sum _i\frac{-1}{2{\mu }_i}\frac{{\partial }^2}{\partial r_i^2}.\hfill \end{array}$$

(12)

$\theta$ represents all of the bend coordinates and r represents all of the stretch coordinates. The functions $B_i\left(r\right)$ and the operators $T_b^{\left(i\right)}\left(\theta \right)$ are known [55,56,58]. One constructs contracted bend functions from a Hamiltonian obtained by fixing the stretch coordinates at some reference geometry and contracted stretch functions from a Hamiltonian obtained by fixing all the bend coordinates at reference values. Products of the bend contracted functions and stretch contracted functions are the final basis functions.

The reduced-dimension Hamiltonian for the bend contraction is,

$$H^{\left(b\right)}=T_{ben}(\theta ,r_e)+V(\theta ,r_e).$$

(13)

Its wavefunctions are denoted by

$$X_b\left(\theta \right)=\sum _l{C}_{lb}{f}_l\left(\theta \right)$$

(14)

and the energies by $E_b$. The $f_l$ are primitive bend basis functions (l is a composite index) and the number of retained bend wavefunctions is denoted by $n_b$. Similarly, the reduced-dimension Hamiltonian for the stretch contraction is,

$$H^{\left(s\right)}=T_{str}\left(r\right)+V({\theta }_e,r).$$

(15)

with the wavefunctions denoted by,

$$Y_s\left(r\right)=\sum _{\alpha }{D}_{\alpha s}{g}_{\alpha }\left(r\right)$$

(16)

and the energies by $E_s$. The $g_{\alpha }$ are primitive DVR stretch basis functions ($\alpha$ is a composite index representing a multidimensional DVR function) and the number of retained stretch wavefunctions is denoted by $n_s$. ${\theta }_e$ and $r_e$ represent reference (often equilibrium) values of all the bend coordinates and all the stretch coordinates. The final basis is a product of the retained stretch and bend eigenfunctions

$$|bs\rangle =|X_b\rangle |Y_s\rangle .$$

(17)

The full Hamiltonian is

$$H=H^{\left(b\right)}+H^{\left(s\right)}+\Delta T+\Delta V$$

(18)

where

$$\Delta V(\theta ,r)=V(\theta ,r)-V(\theta ,r_e)-V({\theta }_e,r)$$

(19)

and

$$\Delta T=\sum _i\Delta B_i\left(r\right)T_b^{\left(i\right)}\left(\theta \right)$$

(20)

with

$$\Delta B_i\left(r\right)=B_i\left(r\right)-B_i\left(r_e\right).$$

(21)

In the contracted basis, MVPs for $\Delta T$ and $H^{\left(b\right)}+H^{\left(s\right)}$ are easy [6].

If one uses an finite basis representation (FBR) primitive bend basis and a DVR primitive stretch basis, a matrix element of $\Delta V$ in the product contracted basis is,

$$\begin{array}{ccc}\hfill \langle b^{\prime }{s}^{\prime }|\Delta V(\theta ,r)|bs\rangle & =& \sum _{\genfrac{}{}{0pt}{}{l^{\prime }l}{\alpha }}{D}_{\alpha s^{\prime }}{C}_{l^{\prime }{b}^{\prime }}\langle l^{\prime }|\Delta V(\theta ,r_{\alpha })|l\rangle C_{lb}{D}_{\alpha s}.\hfill \end{array}$$

(22)

This may be re-written

$$\begin{array}{c}\hfill \sum _{\alpha }{F}_{b^{\prime }b,\alpha }{D}_{\alpha s^{\prime }}{D}_{\alpha s},\end{array}$$

(23)

where I have introduced an F matrix [6] defined by,

$$F_{b^{\prime },b,\alpha }=<b^{\prime }|\Delta V(\theta ,r_{\alpha })|b>=\sum _{l^{\prime }l}{C}_{l^{\prime }{b}^{\prime }}{C}_{lb}\langle l^{\prime }|\Delta V(\theta ,r_{\alpha })|l\rangle .$$

(24)

The integral $\langle l^{\prime }|\Delta V(\theta ,r_{\alpha })|l\rangle$ is computed with quadrature. The $F_{b^{\prime }b,\alpha }$ elements are calculated (in parallel) and stored before matrix vector products are evaluated. The $\Delta V$ matrix–vector product

$$u_{b^{\prime }{s}^{\prime }}^{\prime }=\sum _{bs}\langle b^{\prime }{s}^{\prime }|\Delta V|bs\rangle u_{bs},$$

(25)

is done as follows:

$$\begin{array}{ccc}\hfill u_{b\alpha }^{\left(1\right)}& =& \sum _s{D}_{\alpha s}{u}_{bs}\hfill \\ \hfill u_{b^{\prime }\alpha }^{\left(2\right)}& =& \sum _b{F}_{b^{\prime }b\alpha }{u}_{b\alpha }^{\left(1\right)}\hfill \\ \hfill u_{b^{\prime }{s}^{\prime }}^{\prime }& =& \sum _{\alpha }{D}_{\alpha s^{\prime }}{u}_{b^{\prime }\alpha }^{\left(2\right)}.\hfill \end{array}$$

(26)

The idea of reducing the memory cost of contracted-basis calculations by storing a matrix representation of $\Delta V$ was used in [54], where only the bend basis was contracted. The full power of the method is realized only when both stretch and bend bases are contracted [6,59,60,61,62,63,64]. Recently, similar ideas were used for Cl${}^{-}$ -H${}_2$O [65]. Yu has studied several molecules using similar ideas [66,67,68,69]. As in [54], he contracts only the bend part.

5. Using Pruning to Reduce Both Basis and Grid Size

In this section, I present an alternative method for computing vibrational spectra of molecules with more than five atoms. In contrast to the idea of using contracted basis functions, it uses univariate basis functions, but uses only selected products, i.e., the basis of Equation (1) is pruned by removing functions that are deemed unimportant. This makes it possible to obviate the need to store large vectors.

It seems clear that one should discard basis functions that are not necessary. Many authors have implemented basis pruning strategies [16,18,19,53,62,63,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86]. Although pruning has the obvious advantage that it decreases the size of the vectors one must store and the spectral range of the Hamiltonian matrix, if one uses an iterative method, it complicates the evaluation of MVPs. A pruned basis necessarily has less structure than a direct product basis. In this section, I shall discuss how to evaluate MVPs when the pruning strategy retains some product structure.

Pruning is more efficient when used with a FBR and not a DVR [40]. The simplest VBR pruning condition is $n_1+\dots +n_D\le b$. The pruned basis is much smaller than the direct product basis. If $n_c=0,1,\dots ,b$ for $c=1,\dots ,D$ and $b=14$ then the size of the direct product is ∼6 $\times {10}^{11}$, for $D=10$; ∼4 $\times {10}^{17}$, for $D=15$; and ∼3 $\times {10}^{23}$, for $D=20$. On the other hand, if basis functions with $n_1+\dots +n_D>b=14$ are discarded, the basis size increase with D is less than linear: ∼2.0 $\times {10}^6$, for $D=10$, ∼7.7 $\times {10}^7$, for $D=15$, ∼1.4 $\times {10}^9$, for $D=20$. MVPs for the KEO in a pruned basis are straightforward. MVPs for the potential are only straightforward if a direct product quadrature grid is used [83], but if a direct product quadrature is used, one needs to store a potential vector about as large as the direct product vectors one avoids by pruning the basis. The most important advantage of pruning is therefore lost It is possible to find a nondirect product quadrature scheme that uses fewer points and to evaluate potential MVPs by doing sums sequentially. The ideas will be explained with the $n_1+\dots +n_D\le b$ pruning condition, but better pruning conditions will be briefly discussed at the end of this section.

The direct product quadrature is, in a sense, too good because it is so good that many matrix elements with basis functions removed by the pruning are also exact. A nondirect product Smolyak quadrature is better. It has far fewer points but maintains enough structure to allow efficient MVPs. For detail see [7,87,88]. To make a Smolyak quadrature, one needs a family of 1D quadrature rules for each coordinate. Quadratures are labelled by $i_c$, $i_c=1,2,\dots ,i_c^{max}$. The number of points in quadrature rule $i_c$ is $m_c\left(i_c\right)$, where $m_c\left(i_c\right)$ is a non-decreasing function of $i_c$. To evaluate MVPs efficiently, one needs to use points for which all points in rule $i_c-1$ are also in rule $i_c$. Such points are called nested. The standard way to write a Smolyak quadrature is $f(q_1,q_2,q_3,q_4,q_5,q_6)({}^{1}w\left(q_1\right){}^{2}w\left(q_2\right){}^{3}w\left(q_3\right){}^{4}w\left(q_4\right){}^{5}w\left(q_5\right){}^{6}w\left(q_6\right))$ is

$$\begin{array}{c}\hfill \sum _{i_1+\dots +i_6\le H}{C}_{i_1,\dots i_6}^{smol}\sum _{k_1}^{m_1\left(i_1\right)}\sum _{k_2}^{m_2\left(i_2\right)}\sum _{k_3}^{m_3\left(i_3\right)}\sum _{k_4}^{m_4\left(i_4\right)}\sum _{k_5}^{m_5\left(i_5\right)}\sum _{k_6}^{m_6\left(i_6\right)}\\ \hfill {}^{i_1}{w}_{k_1}{}^{i_2}{w}_{k_2}{}^{i_3}{w}_{k_3}{}^{i_4}{w}_{k_4}{}^{i_5}{w}_{k_5}{}^{i_6}{w}_{k_6}f(q_1^{k_1},q_2^{k_2},q_3^{k_3},q_4^{k_4},q_5^{k_5},q_6^{k_6})& ,\end{array}$$

(27)

where $q_c^{k_c}$ is a point in the quadrature labelled by $i_c$, ${}^{i_c}{w}_{k_c}$ is the corresponding weight and the 1D quadratures are designed to approximate

$$\begin{array}{c}\hfill \int d{q}_c{}^{c}w\left(q_c\right)f\left(z_c\left(q_c\right)\right).\end{array}$$

(28)

$C_{i_1,\dots i_6}^{smol}$ are coefficients; see [7]. H is increased until convergence is achieved. The union of the grids for which $i_1+\dots +i_6\le H$ is satisfied is called the Smolyak grid. The size of the Smolyak grid is orders of magnitude smaller than the direct product grid size.

It would be costly to use Equation (27) in the evaluation of MVPs because it would be necessary to evaluate the sum over $i_1+\dots +i_6\le H$ for each MVP. When the 1D quadrature rules are nested, one can [7] replace Equation (27) with

$$\begin{array}{c}={\sum }_{k_1}^{N_1}{\sum }_{k_2}^{N_2}{\sum }_{k_3}^{N_3}{\sum }_{k_4}^{N_4}{\sum }_{k_5}^{N_5}{\sum }_{k_6}^{N_6}w(k_1,k_2,k_3,k_4,k_5,k_6)\\ \times f(q_1^{k_1},q_2^{k_2},q_3^{k_3},q_4^{k_4},q_5^{k_5},q_6^{k_6})& ,\end{array}$$

(29)

where

$$\begin{array}{c}\hfill w(k_1,\dots ,k_6)=\sum _{i_1+\dots i_6\le H}{C}_{i_1,\dots ,i_6}^{smol}{}^{i_1}{w}_{k_1}\dots {}^{i_D}{w}_{k_6},\end{array}$$

(30)

are “super weights” that are pre-computed [89]. $N_c$ is a maximum number of points for coordinate c [7]. $N_c$ depends on $k_{c^{\prime }}$ if $c>c^{\prime }$ and $N_1$ does not depend on $k_1,\dots ,k_D$. Using the super weights, it is possible to evaluate a potential MVP by doing sums sequentially,

$$\begin{array}{c}\hfill u^{\prime }(n_6^{\prime },n_5^{\prime },n_4^{\prime }{n}_3^{\prime },n_2^{\prime },n_1^{\prime })=\sum _{k_1=1}^{N_1}{T}_{n_1^{\prime }{k}_1}\sum _{k_2=1}^{N_2}{T}_{n_2^{\prime }{k}_2}\sum _{k_3=1}^{N_3}{T}_{n_3^{\prime }{k}_3}\\ \hfill \sum _{k_4=1}^{N_4}{T}_{n_4^{\prime }{k}_4}\sum _{k_5=1}^{N_5}{T}_{n_5^{\prime }{k}_5}\sum _{k_6=1}^{N_6}{T}_{n_6^{\prime }{k}_6}\\ \hfill w(k_1,k_2,k_3,k_4,k_5,k_6)V(q_1^{k_1},q_2^{k_2},q_3^{k_3},q_4^{k_4},q_5^{k_5},q_6^{k_6})\\ \hfill \sum _{n_6=0}^{n_6^{\max }}{T}_{n_6{k}_6}\sum _{n_5=0}^{n_5^{\max }}{T}_{n_5{k}_5}\sum _{n_4=0}^{n_4^{\max }}{T}_{n_4{k}_4}\sum _{n_3=0}^{n_3^{\max }}{T}_{n_3{k}_3}\sum _{n_2=0}^{n_2^{\max }}{T}_{n_2{k}_2}\sum _{n_1=0}^{n_1^{\max }}{T}_{n_1{k}_1}\\ \hfill u(n_6,n_5,n_4,n_3,n_2,n_1),\end{array}$$

(31)

where $T_{nk}=h_k^{-1/2}{p}_k\left(z\left(q_k\right)\right)$. $n_c^{max}$ depends on $n_{c^{\prime }}$ if $c<c^{\prime }$.

To use Equation (31), one first sums over $n_1$ to compute an intermediate vector $y{1}_{k_1,n_2,n_3,n_4,n_5,n_6}$ and then sums over $n_2$ to compute an intermediate vector whose components are $y{2}_{k_1,k_2,n_3,n_4,n_5,n_6}$ etc. At each step, the $n_c$ and $k_c$ indices are constrained among themselves. Everything is clearly explained in [90]. The Smolyak grid is a sum of smaller direct product grids and therefore has structure that makes it possible to evaluate MVPs by doing sums sequentially. It is also possible to do sums sequentially for any pruning condition of the form $g^1\left(n_1\right)+\dots +g^D\left(n_D\right)\le b$ [87]. Sometimes, $g^c\left(n_c\right)={\alpha }_c{n}_c$ with ${\alpha }_c=⌊\frac{{\omega }_c}{{\omega }_{lowest}}+0.5⌋$ is a good choice. Often this choice can be improved. In general, basis functions with many non-zero indices for coordinates with large frequencies are unimportant. They can be pushed out of the basis by using $G^c\left(n_c\right)>{\alpha }_c^P{n}_c$. In general, it is important to include in the basis product functions for which there are several non-zero indices for coordinates with small frequencies. Such functions are preferentially included in the basis by using $G^c\left(n_c\right)<{\alpha }_c^P$. These are general guidelines which we have found useful, but they are not specific [87,88].

6. Using Rank Reduction to Avoid Storing Full Dimensional Vectors

Contraction (Section 4) and pruning (Section 5) enable one to avoid storing vectors with as many elements as the direct product basis. For molecules with more than five atoms, this is essential. For example, if $D=12$ and $n=10$ then for a single vector one needs ∼8000 GB of memory to store a vector with $n^D$ components. In this section, I describe another approach for avoiding vectors with $n^D$ components. It $does$ use a direct product basis but exploits advantages of a SOPs PES [91,92,93]. The key idea is that in some cases, the $n^D$ coefficients, used to represent a function, can be computed from a much smaller set of numbers. For example, a product of functions of a single coordinate, ${\varphi }_1\left(q_1\right){\varphi }_2\left(q_2\right)\dots {\varphi }_D\left(q_D\right)$, can be represented as

$$\sum _{i_1=1}^n{f}_{i_1}^{\left(1\right)}{\theta }_{i_1}^1\left(q_1\right)\sum _{i_2=1}^n{f}_{i_2}^{\left(2\right)}{\theta }_{i_2}^2\left(q_2\right)\dots \sum _{i_D=1}^n{f}_{i_D}^{\left(D\right)}{\theta }_{i_D}^D\left(q_D\right)$$

and it is only necessary to store $Dn$ numbers.

We have developed computational methods that solve the Schroedinger equation by projecting into a basis of functions that are sums of products. Although the functions in this basis are SOPs, the basis is not a direct product basis. A single basis function with R terms is determined by only $RDn$ numbers. If D is large, this is much less than $n^D$. The key idea is to use basis functions that are sums of products of optimised factors. It works if the Hamiltonian is itself a SOPs. The SOPs basis functions are represented in a primitive product basis made from 1-D functions ${\theta }_{i_j}^j\left(q_j\right)$ with $i_j=1,\dots ,n_j$ for each coordinate $q_j$. Any function can be expanded in this basis as

$$\mathsf{\Psi }(q_1,\dots ,q_D)\simeq \sum _{i_1=1}^{n_1}\dots \sum _{i_D=1}^{n_D}{F}_{i_1{i}_2\dots i_D}\prod _{j=1}^D{\theta }_{i_j}^j\left(q_j\right).$$

(32)

The goal is to avoid explicitly introducing $F_{i_1{i}_2\dots i_D}$. This is possible if $\mathsf{\Psi }(q_1,\dots ,q_D)$ is a SOPs. In that case,

$$F_{i_1{i}_2\dots i_D}=\sum _{\ell =1}^R\prod _{j=1}^D{f}_{i_j}^{(\ell ,j)}$$

(33)

where $f^{(\ell ,j)}$ is a one-dimensional vector associated with the ℓ-th term and coordinate j. The SOPs format for multidimensional functions is known as the canonical polyadic (CP) decomposition for tensors [94,95,96].

Basis functions are made by applying the Hamiltonian. We have used a shifted block power method [91]. It is imperative that every basis vector be in the form of Equation (33). This is only the case if the Hamiltonian is of the form,

$$H(q_1,\dots ,q_D)=\sum _{k=1}^T\prod _{j=1}^D{h}_{kj}\left(q_j\right),$$

(34)

where $h_{kj}$ is a one-dimensional operator acting in a Hilbert space associated with coordinate $q_j$. PES can be forced into SOPs form by using, for example, potfit [5,97], multigrid potfit [98], or neural network methods [99,100,101].

When $\mathbf{H}$ is applied to a vector $\mathbf{F}$ to obtain a new vector ${\mathbf{F}}^{\prime }$, the number of terms in the vector increases. If there are T terms in $\mathbf{H}$, the rank (number of terms) of ${\mathbf{F}}^{\prime }$ is a factor of T larger than the rank of $\mathbf{F}$. All vectors have the form

$$F_{i_1{i}_2\dots i_D}=\sum _{\ell =1}^R{s}_{\ell }\prod _{j=1}^D{\tilde{f}}_{i_j}^{(\ell ,j)}\mathrm{with}\sum _{i_j}^{n_j}{|{\tilde{f}}_{i_j}^{(\ell ,j)}|}^2=1,$$

(35)

where, for each term (ℓ) and each coordinate (j), ${\tilde{f}}_{i_j}^{(\ell ,j)}$ is a normalized 1-D vector, $s_{\ell }$ is a normalization coefficient, and $n_j$ is the number of basis functions for coordinate j. $\mathbf{H}$ can be applied to $\mathbf{F}$ by evaluating 1-D matrix–vector products,

$${\left(\mathbf{H}F\right)}_{i_1^{\prime }\dots i_D^{\prime }}=\sum _{i_1,i_2,\dots ,i_D}\sum _{k=1}^T\prod _{j^{\prime }=1}^D{\left({\mathbf{h}}_{k{j}^{\prime }}\right)}_{i_{j^{\prime }}^{\prime }{i}_{j^{\prime }}}\sum _{\ell =1}^R\prod _{j=1}^D{s}_{\ell }{\tilde{f}}_{i_j}^{(\ell ,j)}.$$

(36)

$$=\sum _{k=1}^T\sum _{\ell =1}^R\prod _{j=1}^D\sum _{i_j}{\left({\mathbf{h}}_{kj}\right)}_{i_j^{\prime }{i}_j}{s}_{\ell }{\tilde{f}}_{i_j}^{(\ell ,j)},$$

(37)

where ${\left({\mathbf{h}}_{kj}\right)}_{i_j,i_j^{\prime }}=\langle {\theta }_{i_j}^j|h_{kj}|{\theta }_{i_j^{\prime }}^j\rangle$.

To avoid having vectors with many terms, we must reduce the rank. To do this, we replace $F_{i_1{i}_2\dots i_D}^{\mathrm{old}}$

$$\begin{array}{ccc}\hfill F_{i_1{i}_2\dots i_D}^{\mathrm{old}}& =\hfill & \sum _{\ell =1}^{R_{\mathrm{old}}}{s}_{\ell }\prod _{j=1}^D{{}^{\mathrm{old}}\tilde{f}}_{i_j}^{(\ell ,j)}\hfill \\ & ⟹\hfill & F_{i_1{i}_2\dots i_D}^{\mathrm{new}}=\sum _{\ell =1}^{R_{\mathrm{new}}}{s}_{\ell }\prod _{j=1}^D{{}^{\mathrm{new}}\tilde{f}}_{i_j}^{(\ell ,j)},\hfill \end{array}$$

(38)

where $R_{\mathrm{new}}<R_{\mathrm{old}}$ and choose ${}^{\mathrm{new}}{\tilde{f}}_{i_j}^{(\ell ,j)}$ to minimize $\parallel F^{\mathrm{new}}-F^{\mathrm{old}}\parallel$. We use the same $R_{\mathrm{new}}$ for all reductions. An alternating least squares (ALS) algorithm described in [96] is used to carry out the reduction.

In [91], it was demonstrated that these ideas work for a 20D Hamiltonian of coupled oscillators. A rank of only 20 (i.e., 20 terms in each of the basis functions) was sufficient to converge about 40 states of a 20D Hamiltonian. Related ideas were used successfully for molecules with as many as 10 atoms [91,92,93].

7. Conclusions

Iterative eigensolvers make it possible to calculate vibrational spectra without storing a Hamiltonian matrix. The Lanczos algorithm, filter diagonalization, and a re-started Lanczos or Arnoldi method available as ARPACK [102] are common iterative eigensolvers.

It is easiest to use an iterative eigensolver when the Hamiltonian is a SOP. A Taylor series potential is a SOPs. In many cases, the vibrational KEO is a SOPs. In normal coordinates, it is only a SOPs if one expands elements of the effective moment of inertia tensor [103] or sets them to zero (approximates the KEO). If the PES is not a SOPs, it can be massaged into SOPs form [101,104]. It is harder to use an iterative eigensolver when the Hamiltonian is not an SOPs. In this case, quadrature (or collocation) is used and the Hamiltonian matrix is usually not sparse. In a product basis, the cost of Hamiltonian MVPs scales as $n^{D+1}$ [40,48,49]. This favourable scaling is obtained by exploiting structure. Product basis/product grid methods are methods of first resort for molecules with four or five atoms. They also work extremely well for Van der Waals molecules when only the intermonomer coordinates are treated explicitly [8,11,12,105,106].

In this review article, I describe three methods that obviate the need to store vectors with as many components as the product basis. The first method uses basis functions that are eigenfunctions of a Hamiltonian obtained by setting a subset of the coordinates equal to reference values. In Section 4, a method is described for evaluating MVPs in a product contracted basis. It does not require storing vectors as large as the primitive product basis set. The key idea is to store an intermediate matrix, called the F matrix. The second method uses basis functions that are products of univariate functions. Functions deemed unimportant are removed from a direct-product basis by imposing a pruning condition. The pruning condition is chosen so that the pruned basis has structure. When used with a general PES, it is necessary to use the pruned basis in conjunction with a nondirect product quadrature grid that satisfies two requirements. It must have fewer points than the product quadrature and it must have sufficient structure to be able to evaluate MVPs by doing sums sequentially. A Smolyak quadrature satisfies both requirements. However, to use it, the quadrature must be written not as a sum over quadrature levels, as is usually the case, but as a sum over points (i.e., not Equation (27) but Equation (30)). The third method uses SOPs basis functions. It works only if the Hamiltonian is a SOPs. The basis functions are determined by reducing the rank of the vectors obtained from MVPs. With rank reduction methods, it is possible to compute vibrational spectra for molecules with more than a dozen atoms [107].