Towards breaking the curse of dimensionality in (ro)vibrational computations of molecular systems with multiple large-amplitude motions

TL;DR

This paper introduces a new computational approach combining a numerical kinetic-energy operator and a Smolyak grid method to efficiently solve the (ro)vibrational Schrödinger equation for complex polyatomic molecules with multiple large-amplitude motions.

Contribution

It presents a novel black-box variational method that reduces computational complexity for vibrational calculations of floppy, polyatomic molecules using combined advanced techniques.

Findings

Successfully applied to CH₄·Ar complex in full 12D vibrational space.

Achieved significant reduction in grid size while maintaining accuracy.

Demonstrated feasibility of system-adapted coordinates and iterative eigensolvers.

Abstract

Methodological progress is reported in the challenging direction of a black-box-type variational solution of the (ro)vibrational Schr\"odinger equation applicable to floppy, polyatomic systems with multiple large-amplitude motions. This progress is achieved through the combination of (i) the numerical kinetic-energy operator (KEO) approach of [E. M\'atyus, G. Czak\'o, and A. G. Cs\'asz\'ar, J. Chem. Phys. 130, 134112 (2009)] and (ii) the Smolyak non-product grid method of [G. Avila and T. Carrington, Jr., J. Chem. Phys. 131, 174103 (2009)]. The numerical representation of the KEO makes it possible to choose internal coordinates and a body-fixed frame best suited for the molecular system. The Smolyak scheme reduces the size of the direct-product grid representation by orders of magnitude, while retaining some of the useful features of it. As a result, multi-dimensional (ro)vibrational states are computed with system-adapted coordinates, a compact basis- and grid-representation, and an iterative eigensolver. Details of the methodological developments and the first numerical applications are presented for the CH$_4\cdot$Ar complex treated in full (12D) vibrational dimensionality.

Tables

Table 1. Table 1: Convergence tests for the bound-state vibrational energies of CH ⋅ 4 {}_{4}\cdot Ar (3D) using spherical coordinates, ( R , cos ⁡ θ , ϕ ) 𝑅 𝜃 italic-ϕ (R,\cos\theta,\phi) . The vibrational energies, ν ~ ~ 𝜈 \tilde{\nu} in cm -1 and referenced to the ZPVE, were computed with GENIUSH-DVR and the PES of Ref.  Geleijns  et al. ( 2001 , 2002 ) . The vibrational states are labelled with the (approximate) j 𝑗 j methane angular momentum quantum number, the n R subscript 𝑛 𝑅 n_{R} radial excitation index, and the Γ Γ \Gamma T d subscript 𝑇 d T_{\rm d} (M) irrep label. The R 𝑅 R and ϕ italic-ϕ \phi degrees of freedom are described using generalized Laguerre basis functions ( ℒ n α = 2 subscript superscript ℒ 𝛼 2 𝑛 \mathcal{L}^{\alpha=2}_{n} ) scaled to [ 2.64 , 30 ] 2.64 30 [2.64,30]  Å, and Fourier functions defined over the [ 0 , 2 ​ π ) 0 2 𝜋 [0,2\pi) interval, respectively. Legendre or sincot-Legendre DVR is used for cos ⁡ θ 𝜃 \cos\theta . The number of basis functions and grid points is given for each set as ( K R max , K c max , K ϕ max ) superscript subscript 𝐾 𝑅 max superscript subscript 𝐾 𝑐 max superscript subscript 𝐾 italic-ϕ max (K_{R}^{\rm max},K_{c}^{\rm max},K_{\phi}^{\rm max}) . The test sets, i = A 3 ​ D , B 3 ​ D , 𝑖 subscript 𝐴 3 D subscript 𝐵 3 D i=A_{\rm 3D},B_{\rm 3D}, and C 3 ​ D subscript 𝐶 3 D C_{\rm 3D} are compared with the “final”, benchmark values of F 3 ​ D subscript 𝐹 3 D F_{\rm 3D} , Δ ​ ν ~ i = ν ~ i − ν ~ F 3 ​ D Δ subscript ~ 𝜈 𝑖 subscript ~ 𝜈 𝑖 subscript ~ 𝜈 subscript 𝐹 3 D \Delta\tilde{\nu}_{i}=\tilde{\nu}_{i}-\tilde{\nu}_{F_{\rm 3D}} .

				$\cos \theta$: Legendre-DVR		$\cos \theta$: sincot-Legendre
				$A_{3D}:$ (111,111,31)		$B_{3D}:$ (111,21,17)		$C_{3D}:$ (111,31,31)		$F_{3D}:$ (151,31,31)
$n$	$j$	$n_R$	$\Gamma$	${\tilde{\nu }}_{A_{3D}}$	$\Delta {\tilde{\nu }}_{A_{3D}}$	${\tilde{\nu }}_{B_{3D}}$	$\Delta {\tilde{\nu }}_{B_{3D}}$	${\tilde{\nu }}_{C_{3D}}$	$\Delta {\tilde{\nu }}_{C_{3D}}$	${\tilde{\nu }}_{F_{3D}}$
ZPVE	$0$	$0$	$A_1$	51.200	0.000	51.200	0.000	51.200	0.000	51.200
1	$1$	$0$	$F_2$	9.107	$-$0.002	9.109	0.000	9.109	0.000	9.109
2	$1$	$0$	$F_2$	9.107	$-$0.002	9.109	0.000	9.109	0.000	9.109
3	$1$	$0$	$F_2$	9.109	0.000	9.109	0.000	9.109	0.000	9.109
4	$0$	$1$	$A_1$	29.188	0.000	29.188	0.000	29.188	0.000	29.188
5	$2$	$0$	$F_2$	31.384	$-$0.004	31.388	0.000	31.388	0.000	31.388
6	$2$	$0$	$F_2$	31.384	$-$0.004	31.388	0.000	31.388	0.000	31.388
7	$2$	$0$	$F_2$	31.388	0.000	31.388	0.000	31.388	0.000	31.388
8	$2$	$0$	$E$	31.942	0.000	31.942	0.000	31.942	0.000	31.942
9	$2$	$0$	$E$	31.942	0.000	31.942	0.000	31.942	0.000	31.942
10	$1$	$1$	$F_2$	44.570	$-$0.004	44.573	0.000	44.573	0.000	44.573
11	$1$	$1$	$F_2$	44.570	$-$0.004	44.573	0.000	44.573	0.000	44.573
12	$1$	$1$	$F_2$	44.573	0.000	44.573	0.000	44.573	0.000	44.573
13	$0$	$2$	$A_1$	53.036	0.000	53.036	0.000	53.036	0.000	53.036
14	$2$	$1$	$F_2$	56.228	$-$0.004	56.232	0.000	56.232	0.000	56.232
15	$2$	$1$	$F_2$	56.228	$-$0.004	56.232	0.000	56.232	0.000	56.232
16	$2$	$1$	$F_2$	56.232	0.000	56.232	0.000	56.232	0.000	56.232
17	$2$	$1$	$E$	64.046	0.000	64.046	0.000	64.046	0.000	64.046
18	$2$	$1$	$E$	64.046	0.000	64.046	0.000	64.046	0.000	64.046
19	$3$	$0$	$F_2$	65.825	$-$0.013	65.837	0.000	65.837	0.000	65.837
20	$3$	$0$	$F_2$	65.825	$-$0.013	65.837	0.000	65.837	0.000	65.837
21	$3$	$0$	$F_2$	65.837	0.000	65.837	0.000	65.837	0.000	65.837
22	$1$	$2$	$F_1$	66.066	$-$0.004	66.070	0.000	66.070	0.000	66.070
23	$1$	$2$	$F_1$	66.066	$-$0.004	66.070	0.000	66.070	0.000	66.070
24	$1$	$2$	$F_1$	66.070	0.000	66.070	0.000	66.070	0.000	66.070
25	$0$	$3$	$A_1$	70.313	0.000	70.313	0.000	70.313	0.000	70.313
26	$3$	$0$	$A_1$	73.497	0.000	73.497	0.000	73.497	0.000	73.497
27	$2$	$2$	$F_2$	75.340	$-$0.007	75.347	0.000	75.347	0.000	75.347
28	$2$	$2$	$F_2$	75.340	$-$0.007	75.347	0.000	75.347	0.000	75.347
29	$2$	$2$	$F_2$	75.347	0.000	75.347	0.000	75.347	0.000	75.347
30	$1$	$3$	$F_2$	80.280	$-$0.003	80.283	0.000	80.283	0.000	80.283
31	$1$	$3$	$F_2$	80.280	$-$0.003	80.283	0.000	80.283	0.000	80.283
32	$1$	$3$	$F_2$	80.283	0.000	80.283	0.000	80.283	0.000	80.283
33	$0$	$4$	$A_1$	83.085	0.000	83.085	0.000	83.085	0.000	83.085
34	$1$	$4$	$F_2$	88.186	$-$0.003	88.189	0.000	88.189	0.000	88.189
35	$1$	$4$	$F_2$	88.186	$-$0.003	88.189	0.000	88.189	0.000	88.189
36	$1$	$4$	$F_2$	88.189	0.000	88.189	0.000	88.189	0.000	88.189
37	$2$	$4$	$E$	88.826	0.000	88.826	0.000	88.826	0.000	88.826
38	$2$	$4$	$E$	88.826	0.000	88.826	0.000	88.826	0.000	88.826
39	$0$	$5$	$A_1$	89.427	0.000	89.427	0.000	89.427	0.000	89.427

Table 2. Table 2: Deviation of the vibrational energies, cm − 1 superscript cm 1 \text{cm}^{-1} , of the CH 4 molecule obtained with GENIUSH-Smolyak with a pruned basis and grid, from the tightly converged results of Ref.  Wang and Carrington Jr ( 2014 ) , with increasing the b 𝑏 b and H 𝐻 H values in the basis and the grid pruning conditions, Eqs. ( 41 ) and ( 44 ), respectively. In general, H = b + D + 2 ≥ b + D 𝐻 𝑏 𝐷 2 𝑏 𝐷 H=b+D+2\geq b+D was found to be sufficient to converge the results (note that D = 9 𝐷 9 D=9 for isolated methane). The corresponding number of Smolyak points, N Smol subscript 𝑁 Smol N_{\rm Smol} , is also shown.

$n$	Deviation from Ref. Wang and Carrington Jr (2014)							Ref. Wang and Carrington Jr (2014)
$b$:	2	3	4	5	6	7	8
$H_{\text{9D}}$:	13	14	15	16	17	18	19
$N_{\text{Smol}}$:	3 481	11 833	35 929	97 561	241 201	556 707	1 202 691
ZPV	41.18	2.51	0.66	0.57	0.07	0.02	0.02	9651.29
1	47.81	44.46	3.03	0.81	0.65	0.09	0.03	10961.76
2	47.81	44.45	3.03	0.81	0.65	0.09	0.03	10961.76
3	47.82	44.45	3.03	0.81	0.65	0.09	0.03	10961.76
4	45.91	42.46	2.97	0.75	0.61	0.08	0.03	11184.76
5	45.93	42.47	2.97	0.75	0.61	0.08	0.03	11184.76
6	81.12	57.35	46.15	4.79	1.23	0.74	0.14	12238.29
7	75.65	55.48	48.12	4.18	1.14	0.77	0.13	12265.12
8	75.73	55.48	48.11	4.17	1.14	0.77	0.13	12265.13
9	76.42	55.48	48.11	4.17	1.13	0.77	0.13	12265.13
10	65.82	53.24	47.43	3.49	1.00	0.72	0.11	12275.73
11	65.85	53.24	47.43	3.49	1.00	0.72	0.11	12275.74
12	78.12	53.43	43.26	3.87	1.03	0.66	0.11	12481.49
13	78.77	53.48	43.26	3.88	1.03	0.66	0.11	12481.49
14	78.77	53.48	43.27	3.87	1.03	0.66	0.11	12481.49
15	72.37	51.38	45.63	3.57	0.93	0.69	0.10	12497.25
16	72.37	51.38	45.65	3.56	0.94	0.69	0.10	12497.25
17	73.08	51.44	45.65	3.57	0.94	0.69	0.10	12497.26
18	83.68	72.76	15.82	2.88	1.76	0.47	0.12	12568.47
19	86.04	74.55	16.18	2.89	1.79	0.48	0.12	12670.73
20	86.07	74.56	16.18	2.89	1.79	0.48	0.12	12670.73
21	86.07	74.56	16.18	2.89	1.79	0.48	0.12	12670.73

Table 3. Table 3: Intramolecular (methane, “Met”) basis set and grid choices used in the 12D CH ⋅ 4 {}_{4}\cdot Ar vibrational computations with the basis and grid pruning conditions n q 1 + … + n q 9 ≤ b subscript 𝑛 subscript 𝑞 1 … subscript 𝑛 subscript 𝑞 9 𝑏 n_{q_{1}}+\ldots+n_{q_{9}}\leq b and i q 1 + … + i q 9 ≤ H subscript 𝑖 subscript 𝑞 1 … subscript 𝑖 subscript 𝑞 9 𝐻 i_{q_{1}}+\ldots+i_{q_{9}}\leq H , respectively. H = b + D + 2 ≥ b + D 𝐻 𝑏 𝐷 2 𝑏 𝐷 H=b+D+2\geq b+D was found to be sufficient to converge the results ( D = 12 𝐷 12 D=12 ). The number of basis functions, N 𝑁 N , and grid points, K 𝐾 K , are also given for the methane (“Met”, “Smol”), for the intermolecular (“Inter”), and for the full (12D) computations.

	Intramolecular (CH4, 9D)					Intermolecular (3D)			CH${}_4{}^{}\cdot$Ar (12D)
Label	$b$	$H$	$N_{\mathrm{Met}}$	$K_{\mathrm{Smol}}/{10}^3$		$N_{\mathrm{Inter}}/{10}^3$	$K_{\mathrm{Inter}}/{10}^3$		$N_{12D}/{10}^5$	$K_{12D}/{10}^7$
$A$	0	14	1	$0.163$		4.28	6.30		0.0464	0.113
$B$	1	15	10	$0.871$		4.28	6.30		0.464	0.604
$C$	2	16	55	$3.48$		4.28	6.30		2.55	2.41
$D$	3	17	220	$11.8$		4.28	6.30		10.2	8.20

Table 4. Table 4: Vibrational bound-state energies, ν ~ ~ 𝜈 \tilde{\nu} in cm − 1 superscript cm 1 \text{cm}^{-1} , referenced to the ZPVE of CH ⋅ 4 {}_{4}\cdot Ar computed in full (12D) vibrational dimensionality using the GENIUSH program extended with the Smolyak algorithm in the present work. The potential energy was approximated with the sum of the molecule-atom interaction PES of Ref.  Geleijns  et al. ( 2001 , 2002 ) and the isolated methane PES of Ref.  Wang and Carrington Jr ( 2014 ) . The A , B , C , 𝐴 𝐵 𝐶 A,B,C, and D 𝐷 D basis and grid representations defined in Table  3 correspond to an increasing b = 0 , 1 , 2 , 𝑏 0 1 2 b=0,1,2, and 3 value in the methane basis functions’ pruning condition, Eq. ( 41 ). Convergence of the results can be estimated based on the deviation from computation D 𝐷 D , Δ ​ ν ~ i = ν ~ i − ν ~ D Δ subscript ~ 𝜈 𝑖 subscript ~ 𝜈 𝑖 subscript ~ 𝜈 𝐷 \Delta\tilde{\nu}_{i}=\tilde{\nu}_{i}-\tilde{\nu}_{D} (see also Table  LABEL:ch4res ). For comparison, the benchmark 3D computation with a fixed, effective ( v = 0 𝑣 0 v=0 ) methane geometry ( F 3 ​ D subscript 𝐹 3 D F_{\rm 3D} ) is also shown (taken from Table  LABEL:1table ).

	12D							3D
Label	$A$ ($b=0$)		$B$ ($b=1$)		$C$ ($b=2$)		$D$ ($b=3$)	(Table LABEL:1table)
	${\tilde{\nu }}_A$	$\Delta {\tilde{\nu }}_A$	${\tilde{\nu }}_B$	$\Delta {\tilde{\nu }}_B$	${\tilde{\nu }}_C$	$\Delta {\tilde{\nu }}_C$	${\tilde{\nu }}_D$	${\tilde{\nu }}_{F_{3D}}$
ZPV	9695.262	132.242	9604.164	41.144	9600.706	37.686	9563.019	51.200
1	9.398	0.285	9.139	0.026	9.112	$-$0.002	9.113	9.109
2	9.398	0.285	9.139	0.026	9.112	$-$0.002	9.113	9.109
3	9.398	0.285	9.139	0.026	9.112	$-$0.002	9.113	9.109
4	29.275	0.086	29.197	0.008	29.189	$-$0.000	29.189	29.188
5	31.970	0.575	31.447	0.052	31.392	$-$0.003	31.395	31.388
6	31.970	0.575	31.447	0.052	31.392	$-$0.003	31.395	31.388
7	31.970	0.574	31.447	0.052	31.392	$-$0.003	31.395	31.388
8	32.687	0.736	32.016	0.066	31.946	$-$0.004	31.950	31.942
9	32.687	0.736	32.017	0.066	31.946	$-$0.004	31.950	31.942
10	45.042	0.463	44.620	0.041	44.576	$-$0.003	44.579	44.573
11	45.042	0.463	44.620	0.041	44.577	$-$0.003	44.579	44.573
12	45.042	0.463	44.620	0.041	44.577	$-$0.003	44.579	44.573
13	53.156	0.119	53.048	0.011	53.036	$-$0.001	53.037	53.036
14	57.039	0.799	56.313	0.073	56.236	$-$0.005	56.240	56.232
15	57.039	0.799	56.313	0.073	56.236	$-$0.005	56.240	56.232
16	57.039	0.799	56.313	0.073	56.236	$-$0.005	56.240	56.232
17	64.807	0.753	64.122	0.068	64.050	$-$0.004	64.055	64.046
18	64.807	0.753	64.122	0.068	64.050	$-$0.004	64.055	64.046
19	66.819	0.970	65.989	0.141	65.839	$-$0.009	65.848	65.837
20	66.819	0.971	65.989	0.141	65.839	$-$0.009	65.848	65.837
21	66.819	0.970	65.989	0.141	65.839	$-$0.009	65.848	65.837
22	67.414	1.337	66.143	0.066	66.072	$-$0.004	66.076	66.070
23	67.414	1.337	66.143	0.067	66.072	$-$0.004	66.076	66.070
24	67.414	1.337	66.143	0.067	66.072	$-$0.004	66.077	66.070
25	70.705	0.388	70.360	0.043	70.314	$-$0.003	70.317	70.313
26	74.623	1.118	73.597	0.092	73.499	$-$0.006	73.505	73.497
27	76.276	0.920	75.438	0.081	75.351	$-$0.005	75.356	75.347
28	76.276	0.920	75.438	0.081	75.351	$-$0.005	75.356	75.347
29	76.276	0.920	75.442	0.086	75.351	$-$0.005	75.356	75.347
30	80.808	0.517	80.338	0.047	80.287	$-$0.003	80.290	80.283
31	80.808	0.517	80.338	0.047	80.288	$-$0.003	80.291	80.283
32	80.808	0.517	80.337	0.047	80.288	$-$0.003	80.291	80.283
33	83.156	0.067	83.093	0.004	83.088	$-$0.000	83.088	83.085
34	88.844	0.647	88.255	0.057	88.194	$-$0.004	88.197	88.189
35	88.844	0.647	88.254	0.056	88.194	$-$0.004	88.198	88.189
36	88.844	0.647	88.254	0.056	88.194	$-$0.004	88.198	88.189
37	89.588	0.753	88.903	0.068	88.830	$-$0.004	88.835	88.826
38	89.588	0.753	88.902	0.067	88.830	$-$0.004	88.835	88.826
39	89.505	0.017	89.431	$-$0.057	89.488	0.000	89.488	89.427

Equations

$${\hat{H}}_{\mathrm{vib}}{\Psi }_i({\xi }_1,\dots ,{\xi }_D)=E_i{\Psi }_i({\xi }_1,\dots ,{\xi }_D),$$

$${\Psi }_i({\xi }_1,\dots ,{\xi }_D)=\sum _{n=0}^{N_{\mathrm{bas}}}{C}_n^i{\psi }_n({\xi }_1,\dots ,{\xi }_D),$$

$${\psi }_n({\xi }_1,\dots ,{\xi }_D)=\prod _{\chi =1}^D{\psi }_{n_{\chi }}^{(\chi )}({\xi }_{\chi })$$

$${\Psi }_i({\xi }_1,\dots ,{\xi }_D)=\sum _{n_1=0}^{N_1^{\max }}\dots \sum _{n_D=0}^{N_D^{\max }}{C}_{n_1,\dots ,n_D}^i\prod _{\chi =1}^D{\psi }_{n_{\chi }}^{(\chi )}({\xi }_{\chi }),$$

$${\Psi }_i({\xi }_1,\dots ,{\xi }_D)=\sum _{f(n_1,\dots ,n_D)\le b}{C}_{n_1,\dots ,n_D}^i\prod _{\chi =1}^D{\psi }_{n_{\chi }}^{(\chi )}({\xi }_{\chi }),$$

$$\hat{V}(\mathbf{\xi })=\sum _{m_1=0}^{M_1^{\max }}\dots \sum _{m_D=0}^{M_D^{\max }}{A}_{m_1,\dots ,m_D}\prod _{\chi =1}^D{\mathcal{V}}_{m_{\chi }}^{(\chi )}({\xi }_{\chi }).$$

$$⟨{\psi }_{n_{\chi }^{\prime }}^{(\chi )}\mid {\mathcal{V}}_{m_{\chi }}^{(\chi )}\mid {\psi }_{n_{\chi }}^{(\chi )}⟩$$

$$\approx \sum _{k_{\chi }=1}^{K_{\mathrm{quad}}}{w}_{\chi ,k_{\chi }}{\psi }_{n_{\chi }^{\prime }}^{(\chi )}({\xi }_{\chi ,k_{\chi }}){\mathcal{V}}_{m_{\chi }}^{(\chi )}({\xi }_{\chi ,k_{\chi }}){\psi }_{n_{\chi }}^{(\chi )}({\xi }_{\chi ,k_{\chi }}),$$

$$V({\xi }_1,\dots ,{\xi }_D)$$

$$+\sum _{i=1}^D\sum _{j=i}^D{V}^{i,j}({\xi }_i,{\xi }_j)$$

$$+\sum _{i=1}^D\sum _{j=i}^D\sum _{k=j}^D{V}^{i,j,k}({\xi }_i,{\xi }_j,{\xi }_k)$$

$$+\sum _{i=1}^D\sum _{j=i}^D\sum _{k=j}^D\sum _{l=k}^D{V}^{i,j,k,l}({\xi }_i,{\xi }_j,{\xi }_k,{\xi }_l)$$

$$+\sum _{i=1}^D\sum _{j=i}^D\sum _{k=j}^D\sum _{l=k}^D\sum _{m=l}^D{V}^{i,j,k,l,m}({\xi }_i,{\xi }_j,{\xi }_k,{\xi }_l,{\xi }_m).$$

$$\int \dots \int \prod _{{\chi }^{\prime }=1}^D{\psi }_{n_{{\chi }^{\prime }}^{\prime }}^{({\chi }^{\prime })}({\xi }_{{\chi }^{\prime }})V({\xi }_1,\dots ,{\xi }_D)\prod _{\chi =1}^D{\psi }_{n_{\chi }}^{(\chi )}({\xi }_{\chi })\text{d}{\xi }_1\dots \text{d}{\xi }_D$$

$$\approx \sum _{K=1}^{K_{\mathrm{quad}}}{w}_K\prod _{{\chi }^{\prime }=1}^D{\psi }_{n_{{\chi }^{\prime }}^{\prime }}^{({\chi }^{\prime })}({\xi }_{{\chi }^{\prime },k_{{\chi }^{\prime }}})V({\xi }_{1,k_1},\dots ,{\xi }_{D,k_D})\prod _{\chi =1}^D{\psi }_{n_{\chi }}^{(\chi )}({\xi }_{\chi ,k_{\chi }}),$$

$$\int \dots \int F({\xi }_1,\dots ,{\xi }_D)\text{d}{\xi }_1\dots \text{d}{\xi }_D\approx \sum _{k_1=1}^{K_1^{\max }}\dots \sum _{k_D=1}^{K_D^{\max }}{w}_{1,k_1}\dots w_{D,k_D}F({\xi }_{1,k_1},\dots ,{\xi }_{D,k_D}),$$

$$\int w(\xi )F(\xi )\text{d}\xi =\sum _{k=1}^{K_{\mathrm{quad}}^{\mathrm{G}}}{w}_k^{\mathrm{G}}F({\xi }_k^{\mathrm{G}}),\mathrm{for}\mathrm{any}F(\xi )=\sum _{n=0}^{2(K_{\mathrm{quad}}^{\mathrm{G}}-1)}{A}_n{\xi }^n,$$

$$I({\xi }_1,\dots ,{\xi }_D)=\sum _{f(n_1,\dots ,n_D)\le b}{\mathcal{I}}_{n_1,\dots ,n_D}\prod _{\chi =1}^D{\Phi }_{n_{\chi }}^{(\chi )}({\xi }_{\chi }).$$

$$\prod _{\chi =1}^D{\Phi }_{n_{\chi }}^{(\chi )}({\xi }_{\chi })\text{with}f(n_1,\dots ,n_D)\le b$$

$$\prod _{\chi =1}^D{\Phi }_{n_{\chi }}^{(\chi )}({\xi }_{\chi })\text{with}f(n_1,\dots ,n_D)>b.$$

$$\int \dots \int F({\xi }_1,\dots ,{\xi }_D)\text{d}{\xi }_1\dots \text{d}{\xi }_D\approx \sum _{m=1}^{K_{\mathrm{quad}}}{w}_{m}F({\xi }_{1,m_1},\dots ,{\xi }_{D,m_D})$$

$$\text{with}{K}_{\mathrm{quad}}\ll \prod _{\chi =1}^D{K}_{\chi }^{\max }$$

$${\int }_0^{\pi }\text{d}\theta {\int }_0^{2\pi }\text{d}\varphi Y_{l,m}(\theta ,\varphi )Y_{l^{\prime },m^{\prime }}(\theta ,\varphi )\sin (\theta )$$

$$\text{with}l,l^{\prime }\le l^{\max }\text{and}m,m^{\prime }\le l^{\max },$$

$$N_{\mathrm{Leb}}\sim \frac{4}{3}(l^{\max }+1{)}^2$$

$${\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{H}_l(q_1)H_{l^{\prime }}(q_1)H_m(q_2)H_{m^{\prime }}(q_2)H_n(q_3)H_{n^{\prime }}(q_3){\text{e}}^{-q_1^2-q_2^2-q_3^2}\text{d}{q}_1\text{d}{q}_2\text{d}{q}_3$$

$$\text{with the restrictions }l+n+m\le 4\text{ and }{l}^{\prime }+n^{\prime }+m^{\prime }\le 5,$$

$${\hat{T}}_{\mathrm{Pod}}^v=-\frac{1}{2}\sum _{i=1}^D\sum _{j=1}^D{\tilde{g}}^{-1/4}\frac{\partial }{\partial {\xi }_i}{G}_{i,j}{\tilde{g}}^{1/2}\frac{\partial }{\partial {\xi }_j}{\tilde{g}}^{-1/4}$$

$${\hat{T}}_{\mathrm{rearr}}^{\mathrm{v}}$$

$$\text{with}U$$

$$=\frac{1}{32}\sum _{kl=1}^D\left[\frac{G_{kl}}{{\tilde{G}}^2}\frac{\partial \tilde{G}}{\partial {\xi }_k}\frac{\partial \tilde{G}}{\partial {\xi }_l}-4\frac{\partial }{\partial {\xi }_k}\left(\frac{G_{kl}}{\tilde{G}}\frac{\partial \tilde{G}}{\partial {\xi }_l}\right)\right],$$

$${\hat{T}}_{\mathrm{frearr}}^{\mathrm{v}}$$

$$\text{with}{B}_j$$

$$\hat{I}\approx \sum _{n=0}^{N_{\text{aux}}}\mid n⟩⟨n\mid .$$

$$\mid n_1,n_2,n_3⟩,0\le n_1\le N_1^{\max },0\le n_2\le N_2^{\max },0\le n_3\le N_3^{\max },$$

$$v^1={\tilde{g}}^{-1/4}{v}^0$$

Full text

I Introduction

Molecular systems with many vibrational degrees of freedom, including multiple fluxional modes have been challenging for nuclear motion theory (also known as quantum dynamics) for decades. These systems are difficult to handle because

1. they require a curvilinear coordinate representation, for which we might not have an analytic kinetic energy operator (KEO) readily available;
2. their wave functions are spread over multiple wells of the potential energy surface (PES); and 3) assume the evaluation of high-dimensional (sometimes singular) integrals due to the multiple, coupled (curvilinear) internal degrees of freedom.

There are important, high-dimensional molecular systems with multiple, large-amplitude motions. For example, molecular complexes belong to this class. Molecular complexes are prototypes for molecular interactions and they can be probed in high-resolution spectroscopy experiments. Weakly-bound complexes have a shallow PES valley, so they exhibit only a few, low-energy transitions between bound states, but they usually have a rich predissociation spectrum which can be probed in overtone spectroscopy experiments.

The theory of molecular complexes has been restricted to the explicit quantum mechanical description of the inter-monomer modes, while the monomers were held fixed, described with some rigid, effective structure Huang et al. (2008); van der Avoird et al. (2010). An explicit consideration of monomer-flexibility effects Tennyson and Sutcliffe (1983); Zhang et al. (1995) has come to the focus only in recent years Leforestier et al. (2012); Wang and Carrington, Jr. (2017, 2018). This is not surprising: adding the monomer degrees of freedom to the quantum dynamics treatment rapidly increases the vibrational dimensionality, while in molecular complexes, monomer flexibility effects are usually small, so they can be averaged upon a first look at the system. At the same time, the flexibility of monomers, through the kinetic and the potential energy couplings, plays a central role in the energy transfer between the inter- and the intra-molecular degrees of freedom during the (ro)vibrational and collision dynamics.

Motivation for the present work is provided by these ideas, but we hope that the methodological developments described in this article will become useful for solving the (ro)vibrational Schrödinger equation of (high-dimensional, floppy) molecular systems, in general.

II Curse of dimensionality in vibrational computations

We focus in the present work on the variational solution of the Schrödinger equation including the (ro)vibrational Hamiltonian of $D$ vibrational degrees of freedom, $\boldsymbol{\xi}=(\xi_{1},\ldots,\xi_{D})$,

[TABLE]

where the vibrational wave function is as a linear combination of $N_{\rm bas}$ orthogonal basis functions

[TABLE]

and the expansion coefficients $C^{i}_{n}$ are obtained as the elements of the eigenvectors of the Hamiltonian matrix. The Hamiltonian matrix elements are computed with some appropriate (multi-dimensional) integration scheme. If the basis set is well chosen in this finite basis representation (FBR) scheme, the lowest eigenvalues of the Hamiltonian matrix converge to the exact energies by increasing $N_{\rm bas}$. The most common way to build the multi-dimensional basis functions is to use a direct-product ansatz

[TABLE]

constructed from the $\psi^{(\chi)}_{n_{\chi}}(\xi_{\chi})$ orthogonal basis functions.

II.1 Curse of dimensionality due to the multi-dimensional vibrational basis

By adopting a direct-product basis set, the vibrational wave functions are represented as a linear combination

[TABLE]

in which the number of terms (multi-dimensional basis functions) scales exponentially with the vibrational dimensionality, $N_{\rm bas}=\prod_{\chi=1}^{D}(N_{\chi}^{\rm max}+1)$. For low-dimensional systems this is not a problem, but many challenging systems are high dimensional. Our example system, CH${}_{4}\cdot$Ar has twelve vibrational degrees of freedom. For a 12-dimensional (12D) problem, if we pick 10 basis functions per coordinate (a reasonable starting point if the coordinates are equally coupled), the number of product basis functions will be $10^{12}$. In this representation, we would need to store a vector with $10^{12}$ elements to represent a single vibrational state, which would require ca. 7.3 TB of memory in double precision arithmetics. For this reason, beyond ca. 9 vibrational dimensions, it is necessary to develop and use methods which attenuate the curse of dimensionality in the basis set.

There are different strategies for breaking the exponential growth of the vibrational basis. The first option is to improve the quality of the basis functions in order to decrease the number of functions per coordinate, at least for a subset of the coordinates. The second option is to find a way to identify and discard the basis functions from the direct-product basis set, which have little effect on the accuracy of the computed eigenvalues. The first alternative is efficiently realized by the multi-configuration time-dependent Hartree (MCTDH) method Meyer et al. (2009); Beck et al. (2000), the canonical polyadic (CP) approach Leclerc and Carrington (2014); Thomas and Carrington (2017); Bowman et al. (2003) or in a contracted basis representation obtained by solving reduced-dimensionality eigenproblems Bacic and Light (1989); Henderson and Tennyson (1990); Mladenovic (2002). The second alternative is achieved by finding physically motivated restrictions on the basis set indices. These restrictions can be as simple as the selection of an appropriate multi-polyad Avila and Carrington, Jr. (2009, 2011a), $f(n_{1},\ldots,n_{D})\leq b$, for which the wave-function expansion reads as

[TABLE]

This basis-pruning strategy will be used later in this work. More elaborate basis-pruning restrictions are used, for example, in the MULTIMODE program Halverson and Poirier (2015); Brown and Carrington (2016).

II.2 Curse of dimensionality due to multi-dimensional integrals

Reducing the number of the multi-dimensional basis functions solves only half of the problem. In (ro)vibrational computations, multi-dimensional integrals must be evaluated to construct the Hamiltonian matrix.

There are two common ways to cope with the integrals problem. The first option is to expand the Hamiltonian in a Sum-of-Products form (SOP). For example, the potential energy in a SOP form is

[TABLE]

Using the SOP form, multi-dimensional integrals are obtained as the sum of products of 1-dimensional (1D) integrals,

[TABLE]

which is evaluated with a 1D numerical quadrature using the $w_{\chi,k_{\chi}}$ and $\xi_{\chi,k_{\chi}}$ quadrature weights and points, respectively, defined for the $\xi_{\chi}$ coordinate (in this work, we account for the Jacobian in the wave function). The integrals converge to their exact value upon the increase of the number of quadrature points, $K_{\rm quad}$. The SOP form is useful when a small number of terms is sufficient in Eq. (6) to represent the Hamiltonian. This form is usually employed in MCTDH and in the CP method Meyer et al. (2009); Beck et al. (2000); Leclerc and Carrington (2014); Thomas and Carrington (2017); Bowman et al. (2003). There are methods which can find an excellent ‘basis set’ for the SOP representation of the Hamiltonian Jackle and Meyer (1996, 1996). If the SOP representation, however, requires an excessive number of function evaluations over a multi-dimensional grid of the vibrational coordinates, the exponential scale up with the dimension is re-introduced. This feature is related to the fact that a SOP representation of the Hamiltonian can be as expensive as the representation of the multi-dimensional wave function. In any case, an effective way for attenuating this type of curse of dimensionality was proposed in Ref. Jackle and Meyer (1996).

As an alternative to a sum-of-product representation of the Hamiltonian, one can approximate it with a truncated multi-mode expansion of $n$th-order terms Halverson and Poirier (2015); Brown and Carrington (2016); Bowman et al. (2003); Ziegler and Rauhut (2019). For example, a five-mode expansion of the the potential energy is

[TABLE]

This expansion is exact if $n=D$, but under certain circumstances (also depending on the $\xi_{1},\ldots\xi_{D}$ coordinates) it is very well converged with $n<D$. Using this approximation, the integrals are evaluated using a $D=1,2,\ldots,n$ dimensional direct-product Gauss quadrature, and thereby, the curse of dimensionality is attenuated.

If we want to use the Hamiltonian directly, without approximating or expanding it, we have to tackle the direct evaluation of multi-dimensional integrals by multi-dimensional quadrature. In this case, the integral of the potential energy over a multi-dimensional basis set is evaluated as

[TABLE]

where $w_{k}$ is the multi-dimensional quadrature weight for the $\xi_{\chi,k_{\chi}}$ points, we used the condensed summation index $K\leftrightarrow(k_{1},k_{2},\ldots,k_{D})$. The integral approaches its exact value as the $K_{\rm quad}$ number of points is increased. The most common multi-dimensional quadrature is the multi-dimensional direct-product quadrature

[TABLE]

where $w_{\chi,k_{\chi}}$ and $q_{\chi,k_{\chi}}$ $(\chi=1,\ldots,D)$ are the 1D quadrature weights and points for the $\chi$th coordinate. 1D quadrature rules are most often Gauss (G) quadrature rules, which integrate exactly

[TABLE]

and $d=2(K^{\rm G}_{\rm quad}-1)$ is called the (1D) accuracy of the Gauss quadrature.

A multi-dimensional direct-product quadrature integration suffers from a similar curse of dimensionality problem as a multi-dimensional direct-product basis set: the number of quadrature points, $K_{\rm quad}=\prod_{\chi=1}^{D}K_{\chi}^{\rm max}$, increases exponentially with the vibrational dimensionality. To continue the 12D example from the previous section in which we had 10 basis functions per coordinate, we choose 13 quadrature points per coordinate (a reasonable value) to evaluate the integrals. Then, the number of points in a direct-product grid is $2.33\times 10^{13}$. Storage of this many double-precision numbers would require 170 TB.

As it was explained earlier, the curse of dimensionality in the basis set can be attenuated by identifying and discarding the product basis functions, which are not necessary for the desired precision of the vibrational states. Then, we may think about attenuating the curse of dimensionality in the quadrature grid by using grids which have a non-product structure. In general terms, the application of non-product quadrature grids can be justified, if the integrand $I(\xi_{1},\ldots,\xi_{D})$ is smooth, i.e., it can be expanded with respect to a pruned, product basis set:

[TABLE]

For smooth functions, it makes sense to distinguish between necessary product basis functions:

[TABLE]

and non-necessary product basis functions:

[TABLE]

The total number of necessary and non-necessary product basis functions scales exponentially with the dimension and this is the reason why the total number of product quadrature grid points, which integrate the overlap of all these functions exactly, also scales exponentially with the dimension. If we need to integrate accurately only the necessary product basis functions, the number of which does not grow exponentially with the dimensionality, it is possible to find a multi-dimensional quadrature, which integrates exactly only the necessary basis functions and which does not grow exponentially with the dimension. In such an approach, the curse of dimensionality in the integration grid can be attenuated, i.e.,

[TABLE]

during the course of the evaluation of the Hamiltonian terms (without approximating by some expansion). Optimal non-product quadratures exist for special cases, two of them are explained in the following paragraphs.

II.2.1 An optimal, two-dimensional, non-product quadrature

The most popular non-product quadrature grid is probably the Lebedev quadrature designed to integrate spherical harmonics Lebedev and Laikov (1999). Lebedev grids are used in density functional theory Murray et al. (1993) and they have been used also in rovibrational computations Wang and Carrington (2003). In particular, if we want to obtain the exact value of all integrals, related to the overlap of the spherical harmonics functions, by numerical integration

[TABLE]

we would need to use a total number of $2(l^{\rm max}+1)^{2}$ grid points in the two-dimensional direct-product grid composed of Gauss–Legendre quadrature points for the $\theta$ and Gauss–Chebyshev (first kind) quadrature points for the $\phi$ coordinate. Note that in the expansion of the $Y_{l,m}(\theta,\phi)Y_{l^{\prime},m^{\prime}}(\theta,\phi)\sin(\theta)$ integrand in terms of the product-basis functions, one has to comply with the two restrictions, $m\leq l$ and $m^{\prime}\leq l^{\prime}$. By taking into account these two restrictions, a (smaller) non-product quadrature grid, called Lebedev grid, can be constructed for the numerical integration which includes only

[TABLE]

points, instead of the $2(l^{\rm max}+1)^{2}$ points of the 2D direct-product grid. For example for $l^{\text{max}}=5$, there are a total number of 36 spherical harmonics functions. The calculate exactly the overlap of these functions, we would need $2(5+1)^{2}=72$ points in the 2D direct-product grid, whereas it is sufficient to use 50 ($\sim 4/3\cdot(5+1)^{2}=48$) Lebedev points Lebedev and Laikov (1999). Note that there is not any general formula for the Lebedev quadrature, but the weights and points are tabulated for several two-dimensional maximum accuracy values.

II.2.2 An optimal, three-dimensional, non-product quadrature

Our next example is about the calculation of the exact value of the overlap integrals in a numerical integration scheme, for products of harmonic oscillator functions,

[TABLE]

where $H_{n}$ is the $n$th Hermite polynomial. The smallest, 3D Gauss–Hermite direct-product grid, which recovers the exact value for all these integrals contains $5^{3}=125$ points. By explicitly considering the restrictions in Eq. (19), we may realize that there are only 35 different product functions in the integrand. The smallest non-product grid (for a maximum multi-dimensional accuracy of 9), which recovers the exact value of the integrals for the possible integrands consists of only 77 points Stroud (1971). We note that the corresponding Smolyak grid consists of 93 points, which is less than the direct-product grid, but more than the optimal non-product grid.

In spite of the fact that the optimal multi-dimensional, non-product quadratures use the smallest number of points, they have some handicaps. First, the construction of optimal, non-product quadratures may be cumbersome. There are only a limited number of cases for which the optimal multi-dimensional quadrature is tabulated in the literature (in practice, limited to $D=2$ or 3 for the available cases) Stroud (1971): the points and weights are available only for certain types of polynomials and for limited values of a maximum multi-dimensional accuracy. Second, optimal non-product quadratures lack any structure, which is a serious disadvantage in rovibrational applications Wang and Carrington (2003). If a non-product grid has some structure (reminiscent of a direct-product grid), then it can be used to compute sums over the 1D quadrature points sequentially, which is an important algorithmic element in efficient variational vibrational approaches.

II.2.3 The Smolyak scheme for non-product grids with a structure

There is a simple way to construct non-product quadrature grids, first proposed by the Russian mathematician Sergey A. Smolyak. The Smolyak grid may be slightly larger than the optimal non-product grid but it retains some useful features of direct-product grids. The Smolyak scheme was first adopted for solving the (ro)vibrational Schrödinger equation by Avila and Carrington in 2009 Avila and Carrington, Jr. (2009, 2011a) who exploited that the Smolyak grid is built from a sequence of quadrature rules, and its special structure makes it possible to compute the potential and kinetic energy matrix-vector products by doing sums sequentially. It is possible to combine the Smolyak algorithm with optimal non-product grids of Stroud Stroud (1971), i.e., non-product Smolyak quadrature grids of high-dimensional systems can be constructed from sequences of Stroud-kind non-product quadratures (if the desired Stroud quadrature is available). Although the Stroud–Smolyak grids have less structure, they require fewer points than Smolyak quadratures built from 1D quadrature rules. This direction has been pursued in (ro)vibrational computations by Lauvergnat since 2014 Lauvergnat and Nauts (2014).

III Definition of the (ro)vibrational Hamiltonian in GENIUSH

The GENIUSH protocol, as it was proposed in 2009 Mátyus et al. (2009), aimed for the development of a universal and exact procedure for the (near-)variational solution of the (ro)vibrational Schrödinger equation. Its central part is the numerical construction of the kinetic energy terms over a grid—thereby, the burdensome derivation and implementation of the kinetic energy operator for various molecular and coordinate choices was eliminated. The GENIUSH program was developed using the discrete variable representation (DVR) Light and Carrington Jr. (2007), and it suffered from the curse of dimensionality (Section II). The present work aims for the elimination of this bottleneck, both in respect of the basis and the grid representations, using the ideas first described by Avila and Carrington in 2009 Avila and Carrington, Jr. (2009).

III.1 Numerical representation of the kinetic-energy operator

The GENIUSH program determines the KEO coefficients numerically, over a grid, from the user’s definition of the vibrational coordinates, $\xi_{i}\ (i=1,2,\ldots,D)$ (and body-fixed frame definition, which is relevant for rovibrational computations). Arbitrary coordinates and frames can be defined by writing down the Cartesian coordinates (in the body-fixed frame) in terms of the vibrational coordinates, $\xi_{i}\ (i=1,2,\ldots,D)$. From this coordinate conversion subroutine (written by the user if not yet available in the code), the program numerically evaluates the mass-weighted metric tensor, $\boldsymbol{g}\in\mathbb{R}^{(D+3)\times(D+3)}$, from the vibrational and the rotational $\boldsymbol{t}$ vectors over the coordinate grid. The vibrational $\boldsymbol{t}$ vectors are obtained by two-sided finite differences, for which a step size of $10^{-5}-10^{-7}$ atomic units has been used. In principle, the numerical but exact differentiation scheme of Yachmenev and Yurchenko Yachmenev and Yurchenko (2015) (using chain rule sequences and the derivatives of ‘all’ possible elementary functions and thereby extending Ref. Lauvergnat and Nauts (2002)) could also be used to eliminate the numerical differentiation step.

The $\boldsymbol{G}$ matrix is calculated by inverting $\boldsymbol{g}$, $\boldsymbol{G}=\boldsymbol{g}^{-1}\in\mathbb{R}^{(D+3)\times(D+3)}$, over the grid points of the vibrational coordinates. In this notation the last three rows and columns of $\boldsymbol{g}$ and $\boldsymbol{G}$ correspond to the rotational coordinates. The vibrational kinetic-energy operator has usually been written in the Podolsky form Podolsky (1928)

[TABLE]

with $\tilde{g}=\text{det}\boldsymbol{g}$, because it requires the calculation of only first coordinate derivatives. The volume element for this Hamiltonian Podolsky (1928); Mátyus et al. (2009); Mátyus (2018), and for all its rearranged variants, Eqs. (21), (24), (60) appearing later in this article, is $\text{d}V=\prod_{i=1}^{D}\text{d}\xi_{i}$. Ref. Mátyus et al. (2009) also used a general but “rearranged” form of the (ro)vibrational Hamiltonian

[TABLE]

which can be further rearranged to

[TABLE]

This last form was used by Lauvergnat and Nauts in their numerical KEO approach Lauvergnat and Nauts (2002). Eqs. (21)–(23) and (24)–(25) require third-order derivatives of the coordinates, which are obtained in GENIUSH by using quadruple precision arithmetic to ensure numerical stability for the finite differences. All functions appearing next to the differential operators in Eqs. (20)–(25) have been available from the original implementation Mátyus et al. (2009), so we were able to change between different KEO representations, which has turned out to be necessary for this work (vide infra).

As a first step for implementing the Smolyak algorithm, we had to replace the original DVR implementation with FBR, because we wanted to discard functions from the direct product using simple, physical arguments, e.g., to restrict the basis to a certain (multi) polyad, Eq. (5).

It is important to notice that the application of the Podolsky form, Eq. (20), assumes the insertion of multiple (truncated) resolutions of identities in the basis during the construction of the KEO representation. In our earlier DVR applications, this did not cause any problem, but since we are aiming for a compact FBR, an accurate representation of the Podolsky form could be ensured only if an auxiliary basis set was introduced to converge the completeness relation

[TABLE]

For example, in a 3D FBR computation with a basis set

[TABLE]

the matrix-vector products

[TABLE]

would have to be expanded with respect to a larger, basis

[TABLE]

where $m$ is determined by the coordinate-dependence of the $\tilde{g}^{-1/4}$ and $G_{i,j}\tilde{g}^{1/2}$ multi-dimensional functions. For the example of the H2O molecule, $m=4$ was found to be sufficient to compute the first fifty vibrational states. So, in this 3D problem, the use of an auxiliary basis set introduces only a modest increase in the computational cost. For a 12D problem, however, an $m=4$ choice would increase the basis space by two orders of magnitude!

For this reason, we will use (the rearranged and) the fully rearranged form of the KEO, Eqs. (21)–(24), which did not require the introduction of any additional (auxiliary) functions in an FBR computation. Further details concerning the matrix representation of the KEO, including a pragmatic ‘treatment’ of the KEO singularities, ubiquitous in floppy systems, will be explained in Section IV.4.

III.1.1 Definition of the coordinates for the example of CH${}_{4}\cdot$Ar

The vibrational dynamics of the CH${}_{4}\cdot$Ar complex was described using the $\xi_{1}=R\in[0,+\infty)$, $\xi_{2}=\theta\in[0,\pi]$, $\xi_{3}=\phi\in[0,2\pi)$ spherical coordinates, and the nine dimensionless normal coordinates of the isolated CH4 molecule, $\xi_{3+i}=q_{i}\in(-\infty,+\infty)$ ($i=1,\ldots,9$). At the reference structure (necessary to define the normal coordinates), the methane was oriented in the most symmetric fashion in the Cartesian space with the C atom is at the origin (this orientation also ensured that the KEO singularity is not at the equilibrium structure of the complex):

[TABLE]

and $r=r_{\rm CH}^{\rm eq}=2.052\ 410\ 803$ bohr was the equilibrium C–H distance corresponding to the PES of Ref. Wang and Carrington Jr (2014). The GENIUSH program evaluates functions appearing in the KEO from a coordinate conversion routine in which the instantaneous (body-fixed) Cartesian coordinates must be specified in terms of the internal coordinates. The Cartesian positions of the carbon and the hydrogen atoms were calculated from the $q_{1},\ldots,q_{9}$ normal coordinate values as

[TABLE]

where $\epsilon=x,y,z$ and $i=1,2,\ldots 5$, and the Cartesian coordinates of the Ar atom, $\epsilon_{6}$ $(\epsilon=x,y,z)$, were measured from the center of mass of the methane moiety and were obtained as

[TABLE]

In the last step of the calculation of the Cartesian coordinates, the center of mass of the complex was shifted to the origin. The orientetation of the body-fixed frame corresponding to the coordinates just described corresponds to the orientation of the frame used to define the methane’s normal coordinates. A more sophisticated choice of the body-fixed frame can be useful to make rovibrational computations efficient. In the present work however, we focus on the computation of the vibrational states. We used the atomic masses Coursey et al. (2015) $m(\text{H})=1.007\,825\,032\,23$ u, $m(\text{C})=12$ u, and $m(\text{Ar})=39.962\,383\,123\,7$ u throughout this work.

III.1.2 Potential energy surface

Due to the lack of any full-dimensional methane-argon potential energy surface, we used the sum of the 3D intermolecular potential energy surface of Ref. Geleijns et al. (2001, 2002) and the 9D methane PES from Wang and Carrington Wang and Carrington Jr (2014). This setup allows us to study the kinetic coupling of this weakly bound complex. Should a full-dimensional PES become available, the computations can be adapted to it.

IV Implementation of the Smolyak scheme in GENIUSH

IV.1 Pruning the basis functions

For the example of the CH${}_{4}\cdot$Ar complex described with the $(R,\theta,\phi,q_{1},\ldots,q_{9})$ vibrational coordinates defined in Section III.1.1, we chose the following 1D basis functions: $\mathcal{L}^{(\alpha)}$ generalized Laguerre basis functions (with $\alpha=2$) or tridiagonal Morse basis functions for $R$; Legendre basis functions (and variants of them) or Jacobi associated basis functions for $\theta$; Fourier functions, composed of $\cos(n_{\phi}\phi),\sin(n_{\phi}\phi)$, for $\phi$; and harmonic oscillator functions for the $q_{1},\ldots,q_{9}$ methane normal coordinates. As a result, the direct-product expansion of the vibrational wave function can be written as

[TABLE]

This direct-product basis representation, for the typical values of $N_{R}^{\rm max}>10$, $N_{\theta}^{\rm max}>20$, $N_{\phi}^{\rm max}>15$, and $b>2$, would include $>5.90\cdot 10^{7}$ functions. To reduce the basis set size, we prune the basis representation of the methane fragment

[TABLE]

by replacing the 0 and $b$ lower and upper summation limits of each normal coordinate with the basis-pruning condition

[TABLE]

which we call ‘standard’ pruning. This condition is a natural choice for normal coordinates and harmonic oscillator basis functions, which provide a good ‘zeroth-order’ model. This standard pruning, equivalent to choosing a big polyad of states, allows us to discard basis functions, for which the coupling between the intramolecular basis functions (through the full Hamiltonian) is small and for which the zeroth-order energies are very different. The larger the $b$ value in Eq. (35), the more accurate (higher excited) vibrational states of methane are obtained. (If we focused on the computation of highly excited methane states, it would be better to use a more sophisticated pruning condition.) For the intermolecular basis set $\psi_{R,n_{R}}(R)\psi_{\theta,n_{\theta}}(\theta)\psi_{\phi,n_{\phi}}(\phi)$ we do not introduce any pruning, because the selected functions are not close to any zeroth-order approximate basis set for this system, so we cannot discard any of the functions based on simple arguments. Nevertheless, standard pruning of the methane part already reduces the basis set substantially. The storage of one vector in the direct-product basis set with 10 basis functions per coordinate would require ca. 8 TB of memory, while using standard pruning, Eq. (35), reduces this value to 0.39 GB.

IV.2 Pruning the grid with the Smolyak scheme

The GENIUSH program computes the values of the $G_{i,j},B_{i},$ and $U$ multi-dimensional functions of the KEO at multi-dimensional points of the vibrational coordinates. Since we do not use any interpolation procedure to fit $G_{i,j},B_{i},$ and $V+U$ to special analytic functions, a multi-dimensional quadrature grid is necessary to obtain the integrals.

It is straightforward to design non-product quadrature grids for the evaluation of the multi-dimensional integrals of the Hamiltonian operator with the standard basis-pruning condition, Eq. (35). For the example of the CH${}_{4}\cdot$Ar complex (see Sec. III.1.1 for the coordinate definition and Sec. IV.1 for the basis set and the pruning condition), the 12D Smolyak integration operator of order $H$ Avila and Carrington, Jr. (2009, 2011a), is

[TABLE]

and the general grid-pruning condition is

[TABLE]

The $i_{\chi}$th incremental operator is defined as

[TABLE]

with $\hat{Q}^{0}_{\chi}=0$. The action of the $i_{\chi}$th operator, $\hat{Q}^{i_{\chi}}_{\chi}$, on an $f$ function is its (numerical, quadrature) integral:

[TABLE]

corresponding to the $w^{i_{\chi}}_{\chi,m}$ quadrature weights and $q^{i_{\chi}}_{\chi,m}$ quadrature points, $m=1,\ldots,m_{\chi}(i_{\chi})$ within the $i_{\chi}$th grid.

We also note that Eq. (36) can be written as a linear combination of the 1D integration operators (instead of using the incremental operators) as

[TABLE]

which allows us to better understand the structure of the Smolyak grid. The Smolyak quadrature grid is a linear combination of product quadrature grids with different 1D accuracies, while it has a smaller number of points than the product grid $\otimes\prod_{i_{\chi}=1}^{12}\hat{Q}^{i^{\rm max}_{\chi}}_{\chi}=\hat{Q}^{i^{\rm max}_{R}}_{R}\otimes\hat{Q}^{i^{\rm max}_{\theta}}_{\theta}\otimes\hat{Q}^{i^{\rm max}_{\phi}}_{\phi}\otimes\hat{Q}^{i^{\rm max}_{q_{1}}}_{q_{1}}\otimes\ldots\otimes\hat{Q}^{i^{\rm max}_{q_{9}}}_{q_{9}}$, where $i^{\rm max}_{\chi}=H-\sum_{\chi^{\prime}\neq\chi}g^{\chi^{\prime}}(1)$ is determined by the smallest value of the pruning function for the other coordinates, Eq. (44). If a product basis function, $f^{(R)}_{n_{R}}(R)f^{(\theta)}_{n_{\theta}}(\theta)f^{(\phi)}_{n_{\phi}}(\phi)f^{(q_{1})}_{n_{q_{1}}}(q_{1})\ldots f^{(q_{9})}_{n_{q_{9}}}(q_{9})$, can be integrated exactly by the product quadrature grid $\hat{Q}^{i_{R}}_{R}\otimes\hat{Q}^{i_{\theta}}_{\theta}\otimes\hat{Q}^{i_{\phi}}_{\phi}\otimes\hat{Q}^{i_{q_{1}}}_{q_{1}}\otimes\ldots\otimes\hat{Q}^{i_{q_{9}}}_{q_{9}}$, that product basis function is also exactly integrated by the Smolyak quadrature grid $\hat{Q}(12,H)$, because it comprises this smaller product grid Cools et al. (1999).

To ensure accurate integration, we have to tune three factors: a) the pruning function, $g^{\chi}(i_{\chi})$ (which must be a monotonic increasing function); b) the value of $H$ (the larger, the better); and b) the number of grid points, $m_{\chi}(i_{\chi})$, in the 1D grids determined by the smallest possible Smolyak grid which integrates accurately the Hamiltonian for a selected, pruned, multi-dimensional basis set.

For the case of CH${}_{4}\cdot$Ar, the basis-set pruning condition was (Section IV.1)

[TABLE]

i.e., the intermolecular basis was retained in its product form and pruning was introduced for the methane fragment. The corresponding non-product grid includes the intermolecular grid in its product form, and a pruned intramolecular grid implemented using following grid-pruning functions:

[TABLE]

The corresponding $\hat{Q}^{i_{\chi}}_{\chi}$ integration operators are chosen as

[TABLE]

This choice of the integration operators allowed us to use the 12D Smolyak operator for the special case when the first three degrees of freedom are described with a direct-product grid. $\hat{Q}^{M_{R}^{\rm max}}_{R}$, $\hat{Q}^{M_{\theta}^{\rm max}}_{\theta}$, and $\hat{Q}^{M_{\phi}^{\rm max}}_{\phi}$ label the integration operators corresponding to the spherical degrees of freedom, and each of them is constructed with a Gauss quadrature rule with $M^{\rm max}_{\chi}$ points and $d_{\chi}=2M^{\rm max}_{\chi}-1$ maximum accuracy.

The $\hat{Q}^{m_{\chi}(i_{\chi})}_{\chi}\ (\chi=q_{1},\ldots,q_{9})$ operators, corresponding to the normal coordinates, are constructed using a nested Hermite quadrature with a maximum degree of $d_{\chi}(i_{\chi})\geq 2i_{\chi}-1$, and $d_{\chi}(i_{\chi})=1$, 5, 5, 7, 15, 15, 15, 15, 17, 29, 29, 29, 31, 33, 61, 61$,\ldots$ for the $i_{\chi}=1,2,3\ldots$ sequence of Eq. (48) (also note that the same quadrature is used for each dimensionless normal coordinate). Nesting means that all quadrature points of the quadrature rule $\hat{Q}^{j}_{\chi}$ also appear in the quadrature rule $\hat{Q}^{j+1}_{\chi}$. It is important that we need to have nested grids to be able to use a Smolyak quadrature efficiently. For this reason, we always use the smallest grid which is nested, e.g., for $i_{\chi}=2$ we use a three-point quadrature, $m_{\chi}(2)=3$, in Eq. (48), because there is not any nested, two-point Hermite quadrature. Nested Hermite grids are listed in tables, see for example Ref. spa .

In this paragraph, we compare the orders of magnitudes for a direct-product and for a Smolyak grid just defined for the example of CH${}_{4}\cdot$Ar. Let us assume, that we have a direct-product basis set with $0\leq n_{R},n_{\theta},n_{\phi}\leq 9$ functions for the spherical degrees of freedom, and $0\leq n_{q_{1}}\ldots+n_{q_{9}}\leq b=3$ for the methane’s degrees of freedom. The smallest 12D product Gauss grid which gives correctly the overlap integrals for this basis set includes $10^{3}\cdot 4^{9}=2.62\cdot 10^{8}$ points. To integrate the overlap for this basis set exactly, we need to choose $H=15$ for the 12D Smolyak grid, which includes $10^{3}\times 871=8.71\cdot 10^{5}$ points, almost three orders of magnitude less than the 12D direct-product Gauss grid. Certainly, an even more significant reduction in the grid size (in comparison with a direct-product grid) can be achieved, if a larger number of degrees of freedom is included in the pruning Avila and Carrington, Jr. (2011b).

The smallest necessary value of $H$ can be calculated from the basis-pruning condition and the value of $b$ as follows. To compute exactly an overlap integral of $2b$ polynomial degree, it is necessary to have a maximum degree of $2i_{\chi}-1\geq 2b$, i.e., $i_{\chi}\geq b+1/2$. Then, using the grid-pruning condition, Eq. (48) and the fact that $i_{\chi}\geq 1$, we must have $H\geq b+D$, which makes $H\geq 3+12=15$ for a 12D problem with $b=3$. In the numerical applications, we choose an $H$ value slightly larger than this minimal necessary value: $H=b+D+2$ was usually found to be sufficient to converge the results for the example computations (Section V).

It is important that the Smolyak algorithm uses nested sequences of quadrature rules. Nesting ensures that the non-product grid has a special structure. By exploiting this structure, a multi-dimensional integral of a multi-variable function, $F(x_{1},\ldots,x_{D})$, can be re-written as

[TABLE]

where the structure of the Smolyak grid appears in the second equation through the $k_{i}^{\rm max}$ upper summation indexes Avila and Carrington, Jr. (2009, 2011a): $k_{1}^{\rm max}$ depends on $H$; $k_{2}^{\rm max}$ depends on $H$ and $k_{1}^{\rm max}$; $k_{3}^{\rm max}$ depends on $H$ and $k_{1}^{\rm max},k_{2}^{\rm max}$, etc. It is important to notice that the multi-dimensional integral, Eq. (49), can be written in a sequential sums form (second equation in Eq. (49)) only for structured grids, otherwise only the first, computationally more demanding, form is applicable.

IV.3 An efficient matrix-vector product algorithm for computing eigenvalues and eigenvectors with an iterative eigensolver

We develop a method to compute (ro)vibrational states of polyatomic molecules with multiple large-amplitude motions. Probably, the most common way to tackle (ro)vibrational problems is to compute the Hamiltonian matrix elements in FBR, and then diagonalize the Hamiltonian matrix following the pioneering work of Whitehead and Handy Whitehead and Handy (1975). For polyatomic molecules and complexes, the size of the basis set, even if we use a pruned, product basis, may be larger than 100 000 ($10^{5}$), and a corresponding non-product quadrature grid would consists of more than 10 000 000 ($10^{7}$) points. Unless the Hamiltonian matrix is very sparse and the system has a high permutation-inversion symmetry, the ‘traditional’ route of using a direct eigensolver is not feasible for time and memory reasons.

Using iterative eigensolvers is a practical alternative Brunet et al. (1988); Bramley and Carrington, Jr. (1994), which allows us to compute eigenvalues and eigenvectors without storing or even explicitly computing the Hamiltonian matrix elements. The key algorithmic element in relation with iterative eigensolvers, is the efficient multiplication of an input vector with the Hamiltonian matrix.

In this section, we develop an efficient matrix-vector products algorithm in relation with the numerical KEO approach (Section III) and the Smolyak scheme (Section IV.2). The multiplication is made efficient by exploiting the structure of the pruned basis set and the structure of the non-product Smolyak grid. Multiplication with the potential energy matrix is carried out as

[TABLE]

with the condensed indexing of the basis labels, grid labels, and multi-dimensional grid points:

[TABLE]

respectively. $T_{n}(x_{k})$ is the value of the basis function with index $n$ at point $x_{k}$ and $W^{{{\rm S}}}$ collects the multi-dimensional quadrature weights. In the Fortran implementation we use two condensed indexes for the intermediate vectors, labelled with $K_{1\ldots\Delta}\leftrightarrow(k_{1},\ldots,k_{\Delta})$ ‘partial’ grid and the corresponding $N_{\Delta+1,\ldots,D}\leftrightarrow(n_{\Delta+1},\ldots,n_{D})$ ‘partial’ basis index. The operations are performed in parallel using the OpenMP protocol. The $n^{\rm max}_{\chi}$ and $k^{\rm max}_{\chi}$ values for each coordinate, i.e., the structure of the basis and the grid, are determined from the basis and the grid pruning conditions.

For our present numerical example, CH${}_{4}\cdot$Ar, $D=12$ and $\boldsymbol{\xi}=(R,c,\phi,q_{1},\ldots,q_{9})$ (henceforth, we use the short labelling $c=\cos\theta$). According to the basis pruning condition, Eq. (41), the upper summation indexes for the basis labels are

[TABLE]

The grid pruning condition in Eq. (44) determines the structure of the quadrature indexes according to

[TABLE]

where $S(k)$ is the index of the smallest quadrature rule in the nested sequence of Hermite quadratures that contains $k$ points. For the Hermite sequence used in the present work, the $S(k)$ values are obtained from Eq. (48):

[TABLE]

If the FBR method is used for the intermolecular coordinates $R$, $c$, or $\phi$, we use more grid points than basis functions $K^{\rm max}>N^{\rm max}$ in order to get exact integrals (typically, $K^{\rm max}-N^{\rm max}\approx 5$ was sufficient to achieve convergence). If the DVR scheme is used (due to the reasons explained in Section IV.4), then we have the same number of points and functions, so $K^{\rm max}=N^{\rm max}+1$.

IV.4 Singularity concerns and a hybrid DVR-FBR solution

We have numerically identified that KEO we use to describe the CH${}_{4}\cdot$Ar complex has singularities along the $\theta$ spherical angle (also related to $\phi$). These singularities appear at $\theta=0$ and $\theta=\pi$ ($c=\cos\theta=\pm 1$), and they represent a considerable challenge for a non-analytic KEO representation, especially because two kinds of singularities appear:

[TABLE]

This singular property can be discerned from numerical tests with the numerical KEO coefficients and by calculating matrix elements for the $G_{i,j}\partial^{2}/\partial\mathcal{R}_{i}\partial\mathcal{R}_{j}$ terms using an associated Jacobi basis set, $J_{n}^{\alpha,\beta}(c)$, for example.

An obvious way to avoid these types of singular integrals for analytic KEOs would be to use the 2D spherical harmonics functions for $\theta$ and $\phi$. This option is the way to go for tailor-made approaches, but it would destroy the simplicity and generality of a universal (ro)vibrational approach we are developing, especially if there are several groups of spherical coordinates $(\theta_{i},\phi_{i})$ $i=1,2,\ldots$ in the system Wang and Carrington (2003). In particular, the application of spherical harmonics would require the development of special matrix-vector product routines for each $i=1,2,\ldots$ values.

Another possibility would be to use Jacobi associated functions, $J_{n}^{\alpha,\beta}(c)$, with $\alpha$ and $\beta$ close to zero. We could follow this alternative, if an analytic KEO and analytic KEO integrals were available.

Since we develop a universal method for numerical KEO representations, we need to find a multi-dimensional quadrature which allows us to evaluate all the different kinds of integrals appearing in the KEO without knowing its exact, analytic form, but knowing only the characteristic singular behavior, Eq. (58). Let us use Jacobi associated functions with $\alpha=\beta=0.001$ for $c$. Then, we have to find a quadrature rule which integrates exactly the following types of integrals simultaneously

[TABLE]

with $0\leq n,n^{\prime}\leq N$. Gauss-quadrature rules exist for each integral in Eq. (59) separately, but there is not any single Gauss quadrature that integrates exactly all three types of integrals, whereas in the numerical KEO, it is not possible to separate different terms of different singular behavior (which we know again from numerical test calculations). Then, the next logical step is to find a (non-Gauss) quadrature rule of $M$ points that gives exactly all the integrals in Eq. (59) at the same time. We determined such a quadrature using a two-step procedure. First, we optimized the quadrature points with a simplex algorithm and calculated the quadrature weights by solving an overdetermined set of equations; this set of points and weights was refined by optimizing both the quadrature points and weights with the simplex algorithm. Unfortunately, this (non-Gauss) quadrature includes a large number of points $(K\gg N)$ (three times as many as a single Gauss-quadrature rule) and some of the points come extremely close to the singular points at $c=-1$ and $c=1$. Since GENIUSH calculates the $G_{i,j}$ elements through finite differences, the finite step size will place limitations on increasing the number of quadrature points. Due to the large number of points and their accumulation near the singular values, we cannot accept this special quadrature as a practical solution for the problem, but we will use this (non-Gauss) quadrature rule to check the practical ideas we explain in the following paragraphs.

Since we do not have analytic integral expressions, and it is not possible to find any compact (Gauss) numerical integration scheme which ensures exact integration, let us consider approximate integrals (which become accurate at the limit of a large number of points). First of all, non-exact integration, due to the singularities, manifests itself in a non-symmetric matrix representation of the KEO in Eq. (24). Then, instead of aiming for exact integrals (with a compact grid), let’s aim for a symmetric matrix representation at the first place. Construction of a symmetric matrix representation is straightforward by using Legendre-DVR (or the variants of it discussed below) and the inherently more symmetric general KEO in Eq. (21) for $c=\cos(\theta)$. We will ensure a symmetric representation in the same way, as in the original DVR-based GENIUSH implementation Mátyus et al. (2009) (see also Ref. Wei and Carrington, Jr. (1994) concerning the Legendre polynomials) and in its applications to floppy systems Fábri et al. (2013); Mátyus et al. (2014); Fábri et al. (2014); Sarka and Császár (2016); Sarka et al. (2016, 2017), which did not suffer from the present singularity problems but which did suffer from the curse of dimensionality. So, we handle the singular coordinate $c$ as we would do it in GENIUSH-DVR, for the rest of the coordinates, we use FBR.

So, instead of using the fully rearranged KEO, Eq. (24), for which we obtain a non-symmetric matrix representation due to inexact integration (the off-diagonal elements with different basis indexes of $c$ fail to be equal unless they are exactly integrated), we re-write the KEO for the $c$ coordinate into the more symmetric form

[TABLE]

Using this KEO and a hybrid DVR(c)-FBR representation, the Hamiltonian matrix is real, symmetric by construction and the matrix elements for functions with the same $c$ index are the same as the ones we get using the fully rearranged KEO, Eq. (24). We have carried out an additional test for this hybrid DVR-FBR approach. First, we performed a fully FBR computation with the fully rearranged KEO, Eq. (24), using a Jacobi associated basis set for $c$ with $\alpha=0.01$, $\beta=0.01$ and a (non-Gauss) quadrature developed to calculate accurately the integrals of Eq. (59). This non-Gauss quadrature included a very large number of points for $c$, so we could afford only a small basis and grid for the other degrees of freedom. We repeated the computation using the same, small basis set for the non-$c$ coordinates and DVR with the symmetric KEO, Eq. (60), for $c$. The two computations resulted in the same eigenvalues, which provides a numerical test for our practical DVR-FBR approach (of course, the eigenvalues obtained in this way were different from the converged values due to the smallness of the non-$c$ basis set). So, in this sense, using DVR(c)-FBR and the KEO in Eq. (60) has the correct “limiting” (convergence) behavior, while it ensures a symmetric matrix representation by construction.

IV.5 Matrix-vector products in the hybrid DVR-FBR

The matrix-vectors products in the hybrid DVR-FBR are carried out similarly to Eq. (50). In what follows we list the necessary changes in comparison with the fully FBR PES multiplication, Eq. (50) to accommodate the hybrid FBR-DVR representation for the KEO of Eq. (60). We also note that in the hybrid DVR-FBR scheme, the $W^{{{\rm S}}}_{K_{1\ldots D}}$ Smolyak weights were obtained using a quadrature rule for the $c$ coordinate with weights equal to one.

The matrix-vector product for the potential (and the $U$ pseudo-potential) term is carried out as in Eq. (50), but for the $c$ coordinate, we make the following replacements:

[TABLE] 2. 2.

The matrix-vector product for the $\frac{\partial}{\partial c}G_{c,c}\frac{\partial}{\partial c}$ term is calculated as in Eq. (50), but for the $c$ coordinate, we make the following replacements:

[TABLE]

with

[TABLE]

where $\Theta_{n_{c}}(c)$ is the $n_{c}$th (cot-, sincot-)Legendre-DVR function with $K_{c}^{\rm max}=N_{c}^{\rm max}+1$ quadrature points (vide infra). 3. 3.

The matrix-vector product for the $\frac{\partial}{\partial c}G_{c,\mathcal{R}}\frac{\partial}{\partial\mathcal{R}}$ term, where $\mathcal{R}$ is not the $c$ coordinate, is calculated as in Eq. (50) with the following replacements:

[TABLE] 4. 4.

The matrix-vector product for the $\frac{\partial}{\partial\mathcal{R}}G_{\mathcal{R},c}\frac{\partial}{\partial c}$ term, where $\mathcal{R}$ is not the $c$ coordinate, is calculated as in Eq. (50) with the following replacements:

[TABLE] 5. 5.

The matrix-vector product for the $G_{\mathcal{R},\mathcal{R}}\frac{\partial^{2}}{\partial\mathcal{R}^{2}}$ term, where $\mathcal{R}$ is not the $c$ coordinate, is calculated as in Eq. (50) with the following replacements:

[TABLE] 6. 6.

The matrix-vector product for the $B_{\mathcal{R}}\frac{\partial}{\partial\mathcal{R}}$ term, where $\mathcal{R}$ is not the $c$ coordinate, is calculated as in Eq. (50) with the following changes

[TABLE] 7. 7.

The matrix-vector product for the $G_{\mathcal{R}_{s},\mathcal{R}_{t}}\frac{\partial^{2}}{\partial\mathcal{R}^{2}}$ term, where $\mathcal{R}_{s}$ and $\mathcal{R}_{t}$ are not the $c$ coordinate, is calculated as in Eq. (50) with the replacements

[TABLE]

IV.6 Analysis and improvements for the intermolecular representation

To test the convergence properties, and to determine the optimal basis set and grid sizes for our example system, CH${}_{4}\cdot$Ar, we performed reduced-dimensionality computations. Intermolecular (3D) computations were performed with a fixed methane structure corresponding to the effective rotational constant, $B_{v=0}=5.246\,981\,98\ \text{cm}^{-1}$ (and an effective C–H distance of $\langle R(\text{CH})\rangle_{v=0}=1.107\,117\,44$ bohr) obtained with the ground-state vibrational wavefunction of CH4 with pruning condition $b=8$ (see Section IV.7) and using the isolated methane’s PES Wang and Carrington Jr (2014).

IV.6.1 Intermolecular angular representation: Legendre, cot-Legendre and sincot-Legendre DVRs

Since regions near the singularities, Eq. (58), are dynamically relevant for the CH${}_{4}\cdot$Ar complex, using Legendre DVR for the coordinate $c=\cos\theta$ is an inefficient choice: more than 120 points are needed to converge all vibrational bound states of CH${}_{4}\cdot$Ar (3D) within 0.01 cm*-1*.

In 2010, Schiffel and Manthe Schiffel and Manthe (2010) proposed more efficient alternatives to Legendre DVR to be used for the type of singularities we have to tackle. First of all, the quadrature is improved by selecting the quadrature points, different from the Legendre points, as the inverse cotangent of the eigenvalues ($w_{i}$) of the following matrix

[TABLE]

where $L_{n}(c)$ is the $n$th normalized Legendre function. These integrals are calculated exactly using the Gauss–Chebyshev quadrature with a sufficiently large number of points. Using the eigenvectors, $\boldsymbol{A}$, of $\boldsymbol{P}$ the cot-Legendre DVR basis functions are defined as

[TABLE]

and the first derivative matrix, $\boldsymbol{M}$, for the cot-Legendre DVR functions is

[TABLE]

In our test calculations, it was sufficient to use 50 cot-Legendre DVR points to converge all bound states of the CH${}_{4}\cdot$Ar in 3D (within 0.01 cm*-1*) (see also Table LABEL:1table).

Schiffel and Manthe Schiffel and Manthe (2010) continued and proposed further improvements by extending the basis set. They have noticed that some eigenfunctions of the KEO in spherical coordinates have a $\sin(\theta)$ ‘component’ close to the singularities, so they extended the Legendre basis set with sine functions. Their new basis set included $L_{n}(c)$, $n=0,\ldots,N_{c}^{\rm max}-s$ and $\sin(\theta),\ldots,\sin(s\theta)$, where $s=2$ was sufficient (and stable without any over-completeness problems, which would occur for larger $s$ values) in most applications. A corresponding DVR basis set, called ‘sincot-Legendre DVR basis’, is obtained in the following procedure:

Orthogonal basis functions are created from the set $\{L_{n}(x),(n=0,\ldots,N_{c}^{\rm max}-2),\sin\theta,\sin 2\theta\}$ by diagonalizing the corresponding overlap matrix $S^{\rm sin-cos}$. The orthogonal basis functions, $L_{n}^{\rm sin-cos}(c)$, are calculated using the eigenvectors of the overlap matrix. 2. 2.

A $\boldsymbol{P}^{\rm sin-cos}$ matrix is introduced with the elements

[TABLE]

The DVR points are the inverse cotangent of the $w_{i}$ eigenvalues of $\boldsymbol{P}^{\rm sin-cos}$. The sincot-Legendre DVR basis functions are obtained from the eigenvectors of the $\boldsymbol{P}^{\rm sin-cos}$ matrix, collected in $\boldsymbol{A}^{\rm sin-cos}$, as

[TABLE] 3. 3.

The first derivative matrix, $\boldsymbol{M}^{\rm sin-cos}$, for sincot-Legendre DVR is

[TABLE]

The integrals for the $\boldsymbol{S}^{\rm sin-cos}$, $\boldsymbol{P}^{\rm sin-cos}$, and $\boldsymbol{M}^{\rm sin-cos}$ matrices can be calculated analytically using elementary properties of trigonometric functions and they were tabulated in Ref. Schiffel and Manthe (2010).

We used the sincot-Legendre DVR points and the corresponding first derivative matrix elements (as an alternative to Legendre DVR) in the matrix-vector multiplication procedure described in Section IV.4. Our 3D test computations show that it is sufficient to use 21 sincot-Legendre DVR points for coordinate $c$ to converge all the bound states within 0.01 cm*-1* for CH${}_{4}\cdot$Ar, which is a significant reduction compared to the original Legendre DVR which required more than 120 points. The performance of a few different representations for the $c$ coordinate is compared in Table LABEL:1table. In all computations, we used the $\mathcal{L}^{\alpha}_{n}$ generalized Laguerre polynomials (with $\alpha=2$) for $R$, scaled to the $[2.64,30]$ Å interval, and Fourier functions for $\phi$. The number of points used for the $R,\cos\theta,$ and $\phi$ degrees of freedom in the three test sets of the table is

• •

$A_{\rm 3D}$: $(K^{\rm max}_{R},K^{\rm max}_{c},K^{\rm max}_{\phi})=(81,101(\text{L}),101)$ using Legendre (L) DVR for $c$

• •

$B_{\rm 3D}$: $(K^{\rm max}_{R},K^{\rm max}_{c},K^{\rm max}_{\phi})=(61,21(\text{SCL}),17)$ using sincot-Legendre (SCL) DVR for $c$

• •

$C_{\rm 3D}$: $(K^{\rm max}_{R},K^{\rm max}_{c},K^{\rm max}_{\phi})=(61,31(\text{SCL}),31)$ using sincot-Legendre (SCL) DVR for $c$

It is important to observe in Table LABEL:1table that the vibrational states are not perfectly converged even with a very large number (more than 100) of Legendre DVR points. On the contrary, almost perfect results are obtained with less than 30 sincot-Legendre DVR points. Another important observation (relevant for the 12D applications in Section V) is that we can use fewer Fourier basis functions for $\phi$, than (sincot-Legendre) functions for $\theta$ to converge the 3D vibrational energies.

IV.6.2 Intermolecular radial representation: Laguerre and Morse-tridiagonal basis sets

If we choose the $\mathcal{L}_{n}^{(\alpha)}$ generalized Laguerre basis functions (with $\alpha=2$) for the $R$ radial coordinate, we have to use a large number, more than 30, basis functions to converge the vibrational bound states. Since in the present work we focus on the computation of bound states, it is better to use tridiagonal Morse basis set Wei and Carrington, Jr. (1992); Johnson and Reinhardt (1986); Tennyson and Sutcliffe (1982). The parameters of Morse function functions were determined according to the equations in Ref. Wei and Carrington, Jr. (1992) with $D=143.49$ cm*-1*, $\alpha=0.65$, and $\gamma=0.00033$. These parameters were adjusted to obtain 13 functions that recover the exact vibrational energies for the bound states of the radial Hamiltonian

[TABLE]

where $\mu_{{\rm CH}_{4},{\rm Ar}}$ is the reduced mass of methane and argon, and $\theta^{\rm eq}$ and $\phi^{\rm eq}$ are the equilibrium values of the 3D PES. Since the CH${}_{4}\cdot$Ar complex is a very isotropic system, the parameters and the radial basis set determined in this way should be useful over the entire range of the $\theta\in[0,\pi]$ and $\phi\in[0,2\pi)$ coordinates.

IV.7 Analysis of the intramolecular representation: vibrational states of CH4

The vibrational basis set used to describe the intramolecular vibrational dynamics, i.e., vibrations of the methane molecule, was constructed from the harmonic oscillator basis set with the standard pruning condition, $\sum_{k=1}^{9}n_{q_{k}}\leq b$ in Eq. (41), and the Smolyak quadrature with $\sum_{k=1}^{9}g^{q_{k}}(i_{q_{k}})\leq H$ in Eq. (44). Table LABEL:ch4res shows the convergence of the lowest vibrational states by increasing $b$ and $H$.

As to the 12D computation of CH${}_{4}\cdot$Ar, the bound states correspond to the zero-point vibrational state (ZPV) of CH4, we focused on the lowest-energy states of CH4. Of course, more accurate results for the isolated methane molecule can be obtained by increasing the size of the Smolyak grid, which is perfectly feasible for a 9D computation.

In a minimalistic setup (to be transferred for the 12D computations), we chose a representation which allowed us to converge the fundamental vibrational energies within 1 $\text{cm}^{-1}$. In this representation the 9D Smolyak grid includes more than 100 000 points, which is approximately an order of magnitude larger than what is necessary for a meaningful representation of the zero-point vibration.

V Full-dimensional (12D) results for methane-argon

All bound-state vibrational energies were computed for the CH${}_{4}\cdot$Ar complex in full (12D) vibrational dimensionality (Table LABEL:12D). The basis and the grid representations are selected based on the convergence tests carried out for the inter- and intra-molecular representations (Sections IV.6 and IV.7). Concerning the intermolecular representation, it is composed of Morse-tridiagonal basis functions with $N^{\rm max}_{R}=12$, sincot-Legendre-DVR basis functions with $N^{\rm max}_{c}=20$, and Fourier functions with $N^{\rm max}_{\phi}=17$. The number of quadrature points was $K_{R}^{\rm max}=15$, $K_{c}^{\rm max}=21$, and $K_{\phi}^{\rm max}=20$. As to the methane fragment, we used four different intramolecular representations, with $b=0,1,2,$ and 3 values, which allowed us to check the convergence of the ZPVE and the vibrational energies in the full-dimensional treatment.

Table 3 gives an overview of the orders of magnitudes of the basis and the grid representations employed in the final 12D computations. The largest computation (set $D$ in the table) includes 82 002 690 (8.20$\cdot 10^{7}$) quadrature points and 1 021 020 ($1.02\cdot 10^{6}$) basis functions. The numerical KEO terms, Eq. (60), and the PES are stored as double precision reals (in Fortran) at every grid point, which amounts to a $(12\cdot 13/2+12+1)\cdot 8.20\cdot 10^{7}=60$ GB memory usage. The dimensionality of the Lanczos vectors are determined by the number of basis functions, so one Lanczos vector occupies a negligible amount of 8 MB of memory. To multiply a trial vector with the Hamiltonian matrix took ca. 230 seconds on 51 processor cores, and we had to perform ca. 10 000 matrix-vector multiplications to obtain the 40 states reported in Table LABEL:12D using an in-house Lanczos implementation (it might be possible to reduce the number of matrix-vector products with a Lanczos and a pre-conditioning algorithm optimized for the present system).

Based on the isolated-methane test computations (Table LABEL:ch4res) the error in the ZPVE for $b=2$ and $3$ is 41 and 2.5 $\text{cm}^{-1}$, respectively. The vibrational energies of the complex (referenced to the ZPVE) change less than 0.01 $\text{cm}^{-1}$ by increasing the $b$ value from 2 to 3, hence we may accept them as converged for $b=3$. The ZPVE of the complex is probably accurate within a few $\text{cm}^{-1}$ with $b=3$ similarly to the case of the isolated methane (Table LABEL:ch4res). We only note that a full 12D computation with $b=4$ would also be feasible with the current implementation, but it would only change the ZPVE, since the vibrational energies were converged already with $b=3$.

We also show the $b=0$ results, which correspond to a single harmonic oscillator function for methane (the product of the zeroth harmonic oscillator basis functions for $q_{1},\ldots,q_{9}$). Since the present model includes only kinetic coupling (the PES coupling is also probably very small), the deviation of $\tilde{\nu}_{A}$ ($b=0$) and $\tilde{\nu}_{A}$ ($b=3$) is due to the structural differences of methane: the effective structure for the $b=0$ ground-state harmonic oscillator basis function is the equilibrium structure, whereas $b>0$ accounts for structural distortions due to anharmonicity effects. This change is related to the common wisdom in reduced-dimensionality computations of weakly-bound complexes that it is better to use effective (vibrationally averaged) monomer structures than equilibrium monomer structures Jeziorska et al. (2000). In agreement with this prescription, the 3D computation (column $F_{\rm 3D}$ in Table 3) performed with an effective methane structure corresponding to the (isolated) ground-state vibration very well reproduces the 12D result (remember that only kinetic coupling is included in the present computation, due to the lack of a 12D fully coupled PES).

VI Summary, conclusions, and outlook

The numerical kinetic-energy operator (KEO) approach as implemented in the GENIUSH program Mátyus et al. (2009) has been extended with the Smolyak algorithm Avila and Carrington, Jr. (2009, 2011a), which opens a promising route towards variational (ro)vibrational computations for polyatomic systems with multiple large-amplitude motions.

A direct, variational solution of the (ro)vibrational Schrödinger equation of polyatomic systems (without imposing constraints on the coordinates) is difficult due to the high vibrational dimensionality, which generates an exponential growth in the direct-product basis used to represent the wave functions, and an exponential growth in the direct-product grid necessary to calculate integrals of multi-dimensional operator terms in the Hamiltonian.

If coordinates well-suited for the motions in the system and good zeroth-order basis functions can be found for each coordinate, it is not necessary to use a direct-product basis, but a much smaller, ‘pruned’ basis can be defined, the size of which does not scale exponentially with the number of vibrational degrees of freedom. If it is possible to prune a direct-product basis, it is also possible to find a pruned product grid to calculate integrals. The Smolyak scheme of Avila and Carrington Avila and Carrington, Jr. (2009, 2011a) makes it possible to define non-product (Smolyak) grids, which are orders of magnitude smaller than a direct-product grid but which retain some of the practical features of a direct-product grid. Most importantly, Smolyak grids can be efficiently used in computing matrix-vector products and efficient matrix-vector products allow us to compute eigenvalues and eigenfunctions with an iterative (Lanczos) eigensolver without storing or even explicitly computing the Hamiltonian matrix elements.

In the present work, the combination of these ideas with the numerical KEO approach of GENIUSH were elaborated and explained for all stages of the vibrational computation of the floppy CH${}_{4}\cdot$Ar complex treated in full vibrational dimensionality. Due to the highly fluxional nature of this system, regions of the curvilinear coordinate domains above which the KEO has singularities are dynamically important.

In a fully finite basis representation (FBR) treatment of the numerical KEO, the Hamiltonian matrix fails to be Hermitian due to inaccurate integration of the singularities in general coordinates. As a practical way to avoid these singularity problems in FBR, we proposed to use (efficient) DVRs and an inherently symmetric form of the general KEO for the singular coordinate(s), which ensures a symmetric matrix representation by construction and correct limiting (convergence) behavior at the same time. In practice, this hybrid DVR-FBR treatment allows us to converge all bound vibrational states of CH${}_{4}\cdot$Ar.

In general, this hybrid DVR-FBR approach makes it possible to continue using

1. numerical KEOs; and
2. a general and simple starting product basis sets and grids (both pruned according to physically motivated restrictions) for systems with multiple large-amplitude motions; and ultimately, to (further) develop a universal, black-box-type (ro)vibrational procedure practical for polyatomic systems. Extension of the algorithm for $J>0$ rotational quantum number is straightforward, limitations might be set by the memory requirements and the computational time.

We can foresee future possible improvements of the present procedure to (at least partially) eliminate the current bottlenecks in terms of memory usage (storage of the numerical KEO terms over the grid, see for example Ref. Nauts and Lauvergnat (2018)) and perhaps also in terms of the computational cost. Furthermore, the present developments, in particular the fact that the Smolyak grid is several orders of magnitude smaller than the direct product grid, can be combined with the basis-set contraction idea Bacic and Light (1989); Yurchenko et al. (2009); Wang and Carrington, Jr. (2018). With these or other developments, it will become possible to directly access the predissociation spectral range corresponding to the molecule’s fundamental (and lowest overtone) vibrations in weakly or more strongly bound complexes of the size of CH${}_{4}\cdot$Ar, i.e., with $D=12$ or perhaps beyond this value. In general, a careful choice of the coordinate set, the basis, and the grid representation are required to make full use of the ideas combined, developed, and described in the present work. We hope that these ideas will find applications, also beyond the realm of molecular complexes, among high-dimensional molecular systems with multiple large-amplitude motions.

**Supplementary Material

**Definition of the normal coordinates used for the methane fragment is provided in the Supplementary Material.

**Acknowledgment

**Financial support of the Swiss National Science Foundation through a PROMYS Grant (no. IZ11Z0_166525) is gratefully acknowledged. We thank PRACE for a Preparatory Access Grant during 2017–2018, and NIIFI for allocating computer time at The Hungarian Computing Infrastructure (Miskolc node). Xiao-Gang Wang and Tucker Carrington helped us to get through a difficult stage of the development work by sharing their well-tested sincot-Legendre-DVR implementation (codvr.f90) with us Wang and Carrington, Jr. (2012). We are indebted to our colleagues and co-authors, in particular Tucker Carrington, Attila Császár, and Gábor Czakó, with whom we worked together on Refs. Avila and Carrington, Jr. (2009) and Mátyus et al. (2009), and also the colleagues who later contributed to the further developments and successful applications of the initial methods over the past 10 years.