Modelling interactions of cationic dimers in He droplets: microsolvation trends in He_nK₂⁺ clusters

Abstract

We report the results of a detailed theoretical investigation of small K₂⁺-doped He clusters. The structural characteristics and stabilities of such cations are determined from ab initio electronic structure calculations at the MRCI+Q level of theory. The underlying interactions show a multireference character and such effects are analyzed. The interaction potentials are constructed employing an interpolation technique within the inverse problem theory method, while the nuclear quantum effects are computed for the trimers, their spatial arrangements are discussed, and information was extracted on the orientational anisotropy of the forces. We found that energetically the most stable conformer corresponds to linear arrangements that are taking place under large amplitude vibrations, with high zero-point energy. We have further looked into the behavior of higher-order species with various He atoms surrounding the cationic dopant. By using a sum of potentials approach and an evolutionary programming method, we analyzed the structural stability of clusters with up to six He atoms in comparison with interactions energies obtained from MRCI+Q quantum chemistry computations. Structures containing He_n motifs that characterize pure rare gas clusters, appear for the larger K₂⁺-doped He clusters, showing selective growth during the microsolvation process of the alkali-dimer cation surrounded by He atoms. Such results indicate the existence of local solvation microstructures in these aggregates, where the cationic impurity could get trapped for a short time, contributing to the slow ionic mobility observed experimentally in ultra-cold He-droplets.

---

Research in the field of cluster science continues to grow aiming to establish the connection with nanoscale science, as clusters could serve as building blocks for new materials with designed/tailored properties.¹ In some cases the system evolves from finite-size toward bulk properties smoothly allowing extrapolation from few-body to many-body behavior, while sometimes significant changes in the nanoscopic properties occur, enabling different applications.²

Among such clusters, the investigation of nano-sized He droplets has received a great deal of attention during the last few years, being a challenge for both theoreticians and experimentalists.^3–22 Nowadays, He nanodroplets are used as ultracold homogeneous matrices for spectroscopic studies, providing a unique environment for the isolation of otherwise hardly accessible molecules, allowing detailed study of a wide variety of molecular dopants, neutral or charged.^{3,4,8,15,23–26}

From experimental and theoretical evidence it turns out that alkali metal atoms and dimer dopants have a series of interesting properties due to the unusual bonding behavior of such systems.^14,27–34 They both are weakly bound, attached to the surface of the He droplet usually forming a dimple-like deformation on it, and atoms typically desorb upon excitation, or when more alkali atoms are attached, then they surf on the surface eventually forming cold molecules, e.g. high-spin alkali dimers, through collisions.^35–39 When moving to the charged alkali dopants, the interaction between the cationic impurity and the surrounding He atoms is much stronger, leading to the formation of solvation shells, the so-called snowball structure.^{28–30,40,41} Such highly compressed He regions owing to electrostriction, indicate the existence of compact localized microstructures (solid-like shells) surrounded by less localized He atomic regimes (liquid-like shells), and have been widely used to interpret the low mobility of ions experimentally observed.^42,43

Concerning investigations involving charged dopants, information on the underlying intermolecular, markedly orientational ionic forces, and on the microscopic structures and energetics of alkali-cation-He complexes is essential. In this vein, several theoretical studies have been reported on alkali cations in He clusters,^{27,31,40,41,44,45} while there is also a number of them exploring the lighter alkali dimer cationic complexes,^10,46–48 indicating critical dependence on the minor details of the interaction potential, and computational difficulties for carrying out modern quantum chemistry computations in weakly interacting heavy-atom-containing molecular systems. Moreover, recently high-resolution mass spectra of the lighter Na₂⁺ and K₂⁺ alkali dimers in He_n clusters have been obtained by electron ionization of doped He droplets.³⁰ It has been found³⁰ that the size dependence of the experimental ion yield of the clusters with up to 20 He atoms exhibits distinct anomalies that have not been observed previously, indicating that the number of atoms predicted in the first solvation subshell is overestimated by theory.³⁰ Such anomalies have been found to be more pronounced in the Na₂⁺He_n with n = 2 and 6, while in the K₂⁺He_n clusters minor anomalies are observed at n = 2, 5, 9 and 12.³⁰

Given the importance of the underlying interactions in the microsolvation processes of such charged impurities in He-clusters, their accurate description requires detailed knowledge of the electronic properties. The advances in the computational capabilities, both software and hardware, have made very precise calculations of such finite-size complexes containing an increasing number of He atoms affordable. However, the intriguing interplay between weak intermolecular interactions and the role of quantum effects on their binding present a real challenge for the theorists. Thus, the present work is directed at studying the structures and energies of different finite-size He_n–K₂⁺ clusters. We intend to study the forces between the potassium-dimer cation and He atom by high level ab initio multi-reference calculations for constructing reliable potential energy surfaces able to describe high-order complexes. We build up a global analytical ab initio-based model potential for the trimer, based on the reproducing kernel Hilbert-space (RKHS) and the inverse problem theory,⁴⁹ which was employed to calculate the vibrational quantum states and zero-point effects. Furthermore, we adopt a simplified sum of potential approach, up to three-body terms, to represent K₂⁺-doped clusters containing an increasing number of He atoms, and their structuring was investigated by optimal potential configurations through an evolutionary programming (EP) algorithm.^50,51

This article is organized as follows: the next section describes the computational details and methods implemented, as well as the results obtained. It contains four subsections including: the electronic structure calculations of the trimer interaction energies, the construction of its potential energy surface by a general interpolation scheme, the trimer's bound quantum states, the sum of potential representation, structures and energetics of higher order K₂⁺-doped He clusters, respectively. Finally, the last section summarizes some concluding remarks.

1 Computational methods and results

1.1 Ab initio electronic structure calculations: multireference character effects

All ab initio calculations have been performed with the Molpro program.⁵² The potential energy of a He atom interacting with the K₂⁺(²Σ⁺_g) ion are computed using the Jacobi coordinates (r, R, θ), where r and R are the vectors joining the two atoms in the K₂⁺ dimer, and the center of mass of the K₂⁺ with the He atom, respectively, with θ being the angle between them. The interaction energies using the supermolecular approach are given by: , where E_He–K2⁺ is the total energy of the complex, and E^BSSE is the basis set superposition error (BSSE) corrected total energies of the corresponding He and K₂⁺ monomers. The BSSE correction was included for all molecular configurations using the counterpoise method.⁵³

In our calculations, the r bond length of the K₂⁺ dimer ion was kept fixed at its experimental equilibrium value of r_e = 4.4 Å,⁵⁴ while the intermolecular R distances ranged from R = 4.5 to 20 Å, and the angle θ varied between 0 and 90° by steps of 10°. Scalar relativistic effects were accounted for by using relativistic and quasi-relativistic effective small-core potentials (ECPs), such as ECP10MDF(11s,11p,5d,3f)⁵⁵ and ECP10MWB(7s,6p),⁵⁶ respectively, for the K atoms. Core–core and core–valence correlation effects were entirely neglected in the small-core calculations. Also, for the K atoms we also used all electrons basis sets, such as the 6-31G(3df)⁵⁷ and def2-QZVP,⁵⁸ while for the He atom we used the large augmented correlation consistent basis sets aug-cc-pV5Z (AV5Z) and aug-cc-pV6Z (AV6Z).^59,60 Interaction energies were computed at three different levels of theory: second-order Møller–Plesset perturbation theory (MP2), couple-cluster single and double excitations with perturbative triples (CCSD(T)), and internally contracted multireference configuration interaction (MRCI+Q). In the MRCI calculations, the orbitals were optimized in the state-averaged complete active space self-consisted-field (CASSCF) calculations, and the CASSCF wave functions were used as reference functions. The Davidson correction (+Q) was applied to the MRCI energies to account approximately for higher excitations, and to reduce size consistency error. All grid configurations and their corresponding energies were organized and analyzed using DENEB software utilities.⁶¹

In Table 1 we list the calculated interaction energies using the indicated basis sets for He and K atoms at the linear and T-shaped geometries of the HeK₂⁺ complex. We also report energies at the complete basis set (CBS) limit, using for each basis set for the K atoms, the AV5Z and AV6Z basis sets for He atoms to extrapolate the energies to the (approximate) CBS value. We performed extrapolation of the correlation energies utilizing the two-point single inverse power function first introduced by Schwartz⁶², with X = 5 and 6. As it can be seen, by using the relativistic ECPs for the K atoms, the computed energies are found to be in accord with those from all electron calculations for both geometries of the system at each level of calculation. In turn, we could notice that the MP2, CCSD(T) and MRCI+Q interaction energies for the T-shaped geometry show some relatively small differences between them, while for the linear geometry the MRCI+Q energy values are considerably much lower than those of MP2 and CCSD(T) for all basis sets used.

Table 1 Interaction energies (in cm⁻¹) at the MP2, CCSD(T) and MRCI+Q levels of theory using the indicated basis sets for linear and T-shaped configurations of the He-K₂⁺ complex

K/He atom	MP2	CCSD(T)	MRCI+Q
ECP(Basis set)	θ = 0/90°
ECP10MWB(7s,6p)/AV5Z	−24.84/−3.67	−26.86/−4.40	−33.40/−4.27
ECP10MWB(7s,6p)/AV6Z	−25.38/−3.71	−27.42/−4.46	−34.30/−4.33
ECP10MWB(7s,6p)/CBS^56	−26.12/−3.76	−28.18/−4.54	−35.25/−5.21
ECP10MDF(11s,11p,5d,3f)/AV5Z	−25.37/−4.54	−27.64/−5.82	−35.56/−5.60
ECP10MDF(11s,11p,5d,3f)/AV6Z	−25.49/−4.58	−27.36/−5.95	−36.33/−5.81
ECP10MDF(11s,11p,5d,3f)/CBS^56	−25.65/−4.63	−28.02/−6.12	−37.38/−6.09
All e^−(6-31G(3df))/AV5Z	−25.05/−4.15	−26.84/−5.21	−34.07/−5.03
All e^−(6-31G(3df))/AV6Z	−25.23/−4.20	−27.00/−5.30	−34.57/−5.11
All e^−(6-31G(3df))/CBS^56	−25.47/−4.26	−27.26/−5.42	−35.25/−5.21
All e^−(def2-QZVP)/AV5Z	−25.45/−4.34	−27.14/−5.53	−35.56/−5.33
All e^−(def2-QZVP)/AV6Z	−25.63/−4.39	−27.33/−5.61	−35.99/−5.40
All e^−(def2-QZVP)/CBS^56	−25.87/−4.45	−27.59/−5.71	−36.58/−5.49

---

Thus, we performed further single-reference MP2 and CCSD(T) and multi-reference MRCI+Q calculations at various angular orientations of He–K₂⁺(²Σ⁺_g) by minimizing the energy at each R distance using the all-electron 6-31G(3df) and def2-QZVP basis sets, as well as the ECP10MDF(11s,11p,5d,3f) ones for K atoms and the AV6Z for the He atom. Fig. 1 shows the minimum energy path as a function of θ angle between the linear and T-shaped potential wells for such calculations. One can see that the energy values from each ab initio method slightly depend on the type of basis set used here for the K atoms, with the ECP10MDF(11s,11p,5d,3f) data being lower by less than 0.4 cm⁻¹ from the energies obtained with the def2-QZVP basis set. Also, we found small differences, of less than 0.2 cm⁻¹, between the CCSD(T) and MRCI+Q energies (see inset plot in Fig. 1) for angular orientations with θ ≥ 50°, while in contrast, for linear and bent orientations with θ < 50° the interaction energies highly depend on the level of theory employed. In particular, at these configurations the single-reference methods estimate energies of near 10 cm⁻¹ higher than those of the MRCI+Q one. In a first step, we monitored the T₁, D₁ and T₁/D₁-diagnostics,^63,64 as proposed for open shell complexes, from the single reference UCCSD(T) coupled cluster calculations. These diagnostics have been widely used to test the validity of single-reference methods. Such criteria have been useful for qualitative analysis, while quantitative relationships suggested for the possible significant multireference character have been difficult to interpret. In general, the criteria T₁ < 0.02, D₁ < 0.05 and the T₁/D₁ < 0.4 have been proposed for CCSD calculations. We also checked the expectation value of 〈S²−S_z²−S_z〉 of spin contamination in the UCSSD wavefunction, as well as the weight C₀² of the leading configuration in the CASSCF wavefunction⁶⁵ as given by the MRCI computations, with C₀² values smaller than 0.90, indicating the significant multireference character of the system.

Fig. 1 Interaction energies obtained from the indicated MP2, CCSD(T) and MRCI+Q ab initio methods and basis sets along the minimum energy path as a function of angle θ for the HeK₂⁺ (see text).

Such analysis provides similar trends and in Fig. S1 (see ESI†) we present the values of each of the above most widely used CC diagnostics as a function of the R distance, for the selected values of angle θ (see left panel), and as a function of θ with R = R_e(θ) (see right panel). We found that T₁ spans from 0.3 to 0.001, D₁ from 0.5 to 0.002, while T₁/D₁ from 0.9 to 0.6 as θ goes from 0 to 90° and R distances are between 4.5 to 6.5 Å. One can see that both T₁ and D₁ criteria, indicate single-reference character of the system around the corresponding equilibrium configurations, while their ratio, used to examine the homogeneity of the electronic structure, suggests the need for a multireference electron correlation procedure, especially for θ < 60° (see right panels in Fig. S1, ESI†), as well as all diagnostics shown in the left panels of Fig. S1 (ESI†) for selected R and θ configurations. On the basis of these comparisons, and taking into account the MRCI+Q results obtained for the K₂⁺ spectroscopic constants (Table S1 in ESI†),^77–82 we decided to perform the final calculations at the MRCI level of theory using the small-core ECP10MDF pseudo-potentials with the large (11s,11p,5d,3f) basis set for the potassium atom and the aug-cc-pV6Z for the He atom.

1.2 Representation of the potential energy surface: orientational anisotropy effects

For representing the potential energy surface of He_n–K₂⁺ we employed a general interpolation scheme based on the RKHS method proposed by Ho and Rabitz.⁴⁹ Following the RKHS scheme⁴⁹ the interaction energy is written as,where y = cos θ, N_R and N_θ are the number of ab initio calculated points in the R and θ coordinates, respectively. The expressions of the one-dimensional reproducing kernel functions, q^n,m₁ and q₂ for the distance-like, R, and angle-like, θ, variables, respectively⁴⁹ are given by,where, x_> and x_< are the largest and the smallest of the x and x′ values, respectively. The n and m superscripts refer to the order of smoothness of the function, with values n = 2 and m = 3, and at its asymptotic behavior R⁻⁴ is the leading dispersion interaction between the He atom and the K₂⁺ diatom. B is the beta function, ₂F₁ is the Gaussian hyper-geometric function,⁴⁹ and P_l are the Legendre polynomials with l = 0, 2, 4, 6, 8, 10, 12, 14, 16 and 18.

The v_ij coefficients were obtained by solving eqn (1), with V(R_i,θ_j,r_e) being the calculated ab initio MRCI+Q interaction energies at the corresponding (R_i,θ_j,r_e) grid point. Fig. S2 (see the left panel, ESI†) shows the RKHS potential curves together with the ab initio MRCI+Q interaction energies along the R coordinate for each θ angle from 0 to 90°. One can see the angular anisotropy of the PES, with the potential curve for θ = 0 (linear configurations) to present the deepest minimum, while the interaction for other angular orientations becomes less attractive with those of the T-shaped to count for −5.846 cm⁻¹. Further, in Fig. S2 (see the right panel, ESI†) the long range part of the interaction, as obtained from the MRCI+Q data for linear and T-shaped configurations compared very well with the expected asymptotic behavior, with α_He being the polarizability of the He atom given in ref. 66, indicating the quality of the electronic quantum calculations at these potential regions.

In Fig. 2 a 2D contour plot of the RKHS potential surface is displayed. The K₂⁺ bond length is kept fixed at its equilibrium r_e = 4.4 Å, while the He atom is moving in the (θ, R)-plane. The potential has two symmetric minima at linear geometries with R_e = 5.62 Å, and a well-depth value of 36.355 cm⁻¹, while the barrier between them is 30.509 cm⁻¹ higher in energy, and corresponds to a T-shaped configuration with the He atom at R = 5.82 Å from the K₂⁺ molecular ion.

Fig. 2 2D contour plots of the RKHS PES of the He–K₂⁺(²Σ⁺_g) in the (θ, R)-plane for the fixed K₂⁺ bond length at 4.4 Å. The equipotential curves correspond to energies of −36, −35, −32, −27, −19, −13, −9, −7, −6, −5, −2, and −1.5 cm⁻¹.

Table 2 list the interaction energies of the He–K₂⁺(²Σ⁺_g) complex along the minimum energy path for different angular orientations of the He atom computed by the present MRCI+Q, RKHS PES and previously reported studies.^46,47 One can see that RKHS PES describes the He–K₂⁺ interaction in excellent agreement with the MRCI+Q data, with an overall average difference of 0.04 cm⁻¹. Also comparison with the previously reported values^46,47 shows some significant differences along all angular orientations.

Table 2 Comparison of the MRCI+Q interaction energies with the RKHS potential, and previously reported data along the minimum energy path for the He–K₂⁺ complex. δE is the energy difference between the MRCI+Q interaction energy and those obtained from the RKHS potential. Energies in cm⁻¹, angles in deg., and distances in Å

(θ, R)	ΔEMRCI+Q	V(R,θ,re)/δ[scr E, script letter E]	Literature
a From ref. 46. b From ref. 47.
(0, 5.62)	−36.344	−36.355/0.011	—
(0, 5.5)	−35.465	−35.490/0.025	−38^a
(0, 5.4)	−32.353	−32.420/0.068	−70^b
(10, 5.62)	−35.490	−35.484/−0.006	—
(10, —)	—	—	−35.92^a
(20, 5.41)	−32.728	−32.731/0.003	—
(20, —)	—	—	−31.48^a
(30, 5.31)	−27.295	−27.314/0.019	—
(30, —)	—	—	−22.96^a
(40, 5.41)	−19.823	−19.808/−0.015	—
(40, —)	—	—	—14.81^a
(50, 5.51)	−13.685	−13.701/0.016	−
(50, −)	—	—	−8.88^a
(60, 5.62)	−9.800	−9.779/−0.021	—
(60, —)	—	—	−5.92^a
(70, 5.72)	−7.529	−7.558/0.029	—
(70, —)	—	—	−4.4^a
(80, 5.82)	−6.344	−6.308/−0.036	—
(80, —)	—	—	−4.07^a
(90, 5.82)	−5.808	−5.846/0.038	—
(90, 6.1)	−5.633	−5.679/0.046	−4^a
(90, 5.53)	−5.215	−5.273/0.058	−6^b

---

1.3 Bound state calculations: quantum effects

We performed variational quantum nuclear calculations employing the present RKHS PES for the He–K₂⁺(²Σ⁺_g) complex, and the full spectrum of its bound eigenstates was obtained. As the spin-rotation constant of the open shell K₂⁺(²Σ⁺_g) diatom is two orders of magnitude smaller than its rotational constant B, we neglected this coupling, and its effect on the rotational levels of a ²Σ molecule.^67,68 Thus, the corresponding Hamiltonian for the complex is the same as that for a closed-shell ¹Σ diatom, and in Jacobi coordinates it is written as,with μ₁ and μ₂ being the reduced masses of the He–K₂⁺ and K₂⁺, respectively. The ĵ and l̂ are the angular momentum operators associated with the vectors r and R, respectively, leading to a total angular momentum Ĵ = l̂ + ĵ, and V(R, θ, r) is the RKHS potential given in eqn (1).

We are interested in the vibrational bound eigenstates of the whole system, so we carried out variational discrete variable representation (DVR) calculations, following the procedure described in our previous studies.^69–71 The Hamiltonian is represented as a product of radial, f_n(R), and angular, Θ^(JMp)_jΩ, basis functions. For the R coordinate we used a DVR basis set based on the particle in a box eigenfunctions,⁶⁹ while the angular basis functions are eigenfunctions of the parity (p), with M being the projection of J on the space-fixed z-axis, and Ω its projection on the body-fixed z-axis, which is chosen here along the R vector. Specifically, for the present J = 0 calculations we used 100 DVR radial grid points over the range from 4.5 to 20 Å, and orthonormalized Legendre polynomials {P_j(cosθ)} in θ coordinate with up to j = 40, and 41 values for even (p = +1) and odd (p = −1) parity symmetry, respectively, achieving a convergence of 0.005 cm⁻¹ in the calculations. The rotational constant is taken equal to 0.0445 cm⁻¹ for the K₂⁺, and the standard atomic masses of He and K are 4.0026 and 39.0983 μ, respectively. By diagonalizing the corresponding Hamiltonian matrix we obtained the eigenenergies and eigenfunctions.

In Table 3 we list the energies of all 14 vibrational bound states (even(e) and odd(o)) of the HeK₂⁺ using the present MRCI+Q/RKHS surface, and compare them with the previously reported theoretical data⁴⁶ available. As can be seen, the ground and the first excited vibrational states are double degenerated with even and odd parity states having the same energy, while the states with energy values higher than the T-shaped potential barrier show splitting between the even and odd symmetries that increases as the energy increases. By comparing with the previous data available for an MP2/6-31G(3df) surface,⁴⁶ one can see that the present PES supports in total 14 bound states, three more states than the MP2 surface, with their energy spacings being quite different. The MRCI+Q/RKHS surface predicts a shallower potential well with a more stable HeK₂⁺ conformer than the MP2 one. The computed binding energies are −18.67 and −17.47 cm⁻¹ for the MRCI+Q/RKHS and MP2 PESs, respectively, with the corresponding zero-point energies being 17.67 and 20.53 cm⁻¹, accounting for 48.6 and 54% of their potential well-depths.

Table 3 Vibrational energies (in cm⁻¹) for the even/odd parity bound states of the He–K₂⁺ complex (see text) and their comparison with values from ref. 46

No. state (n)	This work	From ref. 46
0	−18.67/−18.67	−17.47/−17.47
1	−11.00/−11.00	−8.57/−8.57
2	−4.98/−4.96	−3.57/−3.57
3	−4.39/−4.34	−2.48/−2.33
4	−2.07/−1–13	−0.73/−0.26
5	−0.94/−0.47	−0.27/—
6	−0.19/−0.16	—

---

In Fig. 3 and 4 we display the radial and angular probability distributions, respectively, for all even parity vibrational states of the He–K₂⁺ cluster. One can clearly see that ground and first excited vibrational states are localized in the linear potential well, while higher lying states (above the T-shaped barrier) are delocalized, with He atoms undergoing large amplitude radial and bending motions. The corresponding probability distributions for the n = 2–6 states are moving to larger R values showing an increasing number of nodes in both R and θ coordinates, presenting a more complicated nodal pattern gradually. Both radial and angular probability distributions of these states, especially of the n = 4–6 ones, show clear differences between those of even and odd parity symmetry (see Fig. S3 in ESI†), with the even state wavefunctions up to n = 4 describing the bent He–K₂⁺ complexes. We should also note that the even and odd n = 5 and 6 states are localized in the angular coordinate, while their radial distributions are broad especially for the even n = 6 state.

Fig. 3 Radial distributions of all even parity bound states of HeK₂⁺ using the present MRCI+Q/RKHS surface.

Fig. 4 Angular distributions of all even parity bound states of HeK₂⁺ using the MRCI+Q/RKHS surface.

1.4 Energetics and structures of higher-order clusters: Microsolvation effects

Following our previous work on larger He-doped clusters^18,72–75 we represent the analytical PES of higher order He_nK₂⁺ using the sum-of-potentials approach based on the sum of the MRCI+Q/RKHS three-body parameterized HeK₂⁺ interactions plus the He–He one,where the corresponding terms are the present MRCI+Q/RKHS potential of the HeK₂⁺ complex. and is the potential function for He₂ given in ref. 76. are the vectors connecting the center of mass of K₂⁺ with each He atom, while θ_i are the angles between each vector and the K₂⁺ axis (see Fig. 5). In turn, we follow the (#T1/T2,#L1/#L2,#B1/B2) notation to label the number of He atoms in T-shaped (T1 for θ = 90° and T2 for θ = 270°), linear (L1 for θ = 0° and L2 for θ = 180°), and for bent configurations, B1 for 90 < θ < 270 and B2 for 270 < θ < 90 for any ϕ value, respectively. In this way one could classify the different structures of higher-order He_nK₂⁺ clusters.

Fig. 5 Coordinate system used for the He_nK₂⁺ clusters.

In order to validate the sum-of-potential approach (see eqn (4)) adapted here, we carried out MRCI+Q/ECP10MDF/AV6Z calculations for the tetratomic He₂K₂⁺ cluster. In Fig. 6 we display the potential energy values obtained from the analytical expression of eqn (4) together with the MRCI+Q computed interaction energies for the He₂K₂⁺ cluster along its R₁ and R₂ distances at selected orientantional configurations, such as those with He atoms in linear and/or T-shaped geometries. In particular, we chose to check interaction energies along three different cluster orientations, the (0/0, 1/1, 0/0), (1/1, 0/0, 0/0) and (1/0, 1/0, 0/0), corresponding to the global and low-energy local minima on He₂K₂⁺ PES. In this way, one could validate the proposed analytical potential approach, as the energies obtained from the sum-of-products expression of eqn (4) are in excellent agreement with the corresponding ab initio MRCI+Q values, using the same basis sets for K and He atoms as in the trimer case for all configurations considered here for the tetratomic He₂K₂⁺ cluster.

Fig. 6 Potential curves obtained for the He₂K₂⁺ cluster using the sum-of-potentials approach of eqn (4), (solid lines), together with the calculated MRCI+Q/AV6Z interaction energies (solid-circle lines) as a function of (a) R₁ for the (0/0, 1/1, 0/0) structure, (b) R₁ =R₂ for the (1/1, 0/0, 0/0) configuration, and (c and d) R₁ and R₂ for the (1/0, 1/0, 0/0) configuration, respectively. The arrows in each configuration indicate which He atoms are moving along the R₁/R₂ coordinates.

In turn, we carried out a systematic search to identify optimal energy structures, such as global and local minima, on the He_nK₂⁺ PESs of higher-order clusters. An evolutionary programming algorithm⁵⁰ is employed to find the lowest energy structures and study the structural optimization of He_nK₂⁺ clusters, with n up to 6 He atoms. Such algorithms have been proven to be an attractive and efficient way to solve numerical optimization problems in multi-dimensional space, as gradients of the potential are not needed, while the search space-narrowing is automatically realized by the self-adaptation of mutation.^18,50,51 Briefly, we start by generating an initial population of [scr M, script letter M] = 10 individuals (for each He_nK₂⁺ cluster in our case). Each individual is characterized by a pair of real vectors (χ_i,η_i)_{i=1−[scr M, script letter M])}, containing the Cartesian coordinates, χ_i, of all cluster atoms and their standard deviation, η_i, for Gaussian mutation (strategy parameter) that controls the evolution of the dispersion of the population in time. The initial coordinates χ_i with η_i = 1 are chosen randomly in the interval (0, Δ), with Δ being a displacement factor that increases the resolution power. Each parent set (χ_i,η_i) evolves to generate a new population by mutation,⁵⁰ while for each individual of the joint parent-child group (2[scr M, script letter M] individuals) q (tournament size equals 10) opponents are randomly chosen from the 2[scr M, script letter M] − 1 individuals, to compare each other, and the individual in each encounter with the lower potential energy wins. The best [scr M, script letter M] individuals then become parents for the next generation and so on. Convergence is achieved when the potential energy difference between two consecutive generations is below a threshold value of 10⁻³ cm⁻¹.

Thus, for each He_nK₂⁺ cluster we localized the optimized structures on its PES (see eqn (4)), and by carrying out additional frequency analysis of the Hessian eigenvalues we could characterize the nature of the stationary PES points. We were able to identify several minima on each cluster's PES, with those lying close in energy to the global minimum being important as they could be accessible by including quantum zero-point effects.

In turn, in Table S2 (see the ESI†) we list the computed energies of various low-energy structures for the He_nK₂⁺ clusters with n up to 6 He atoms, obtained from the EP method⁵⁰ using the sum-of-potential approach (see eqn (4)), together with their comparison with previously reported values available in the literature.^10,46 One can see that the difference in the energy values varies, depending on the size cluster and the optimal configuration, with most energies of the global minima showing differences of a few cm⁻¹, while for the remaining optimal structures we also found differences in their ordering, especially for the cluster with n = 4 and 5 He atoms.

In Fig. 7 we display the configurations of the selected optimal structures for each cluster. For the He₂K₂⁺ we show six optimal structures with the (0/0, 0/0, 2/0) structure being the global minimum on the PES at an energy of −76.5 cm⁻¹, while the linear (0/0, 1/1, 0/0) and (0/0, 1/0, 1/0) ones are found to be very close in energy, within 3.8 and 5.4 cm⁻¹, respectively. Also, all T-shaped configurations are found to be relatively less stable than the others, with the (1/1, 0/0, 0/0) located at an energy of just −11.64 cm⁻¹ (see Fig. 7). For the He₃K₂⁺ cluster five optimal structures are shown with the global minimum (0/0, 0/0, 3/0) structure found at an energy of −123.9 cm⁻¹, while the next structures are higher at relatively close energies of −112.8, −106.1, −85.4 and −78.6 cm⁻¹ corresponding to bent, linear and T-shaped configurations. In turn, the He₄K₂⁺ optimal structures are found well separated in energy with the global minimum at −166.1 cm⁻¹ for the (0/0, 0/0, 4/0) configuration, while the two lower in energy (0/0, 0/0, 5/0)/(0/0, 0/0, 6/0) and (0/0, 0/0, 3/2)/(0/0, 0/0, 3/3) configurations for the He₅K₂⁺/He₆K₂⁺ clusters are close in energy within just ∼3.9 cm⁻¹. By analysing the cluster structures as the number of He increases, one can figure out that the He atoms prefer to locate themselves at one end of the K₂⁺ forming He_n-motifs, such as triangles, rhombs, trapezoids and rings, or symmetrically at its ends forming equilateral He₃ in the case of He_5,6K₂⁺ systems. The stronger He-K₂⁺ interaction at linear configurations forces the sequential He atoms to get attached at the two ends of the dimeric cation. Such selective growth of the small He_nK₂⁺ clusters around the two sides of the cation is due to the competition between the weak He–K₂⁺ and He–He interactions, giving rise to localized He-networks. The formation of such individual compact motifs during the microsolvation process could cause the trapping of the dopant for a short time, and such features could be related to the slow mobility of ions in He-droplets. However, we should point out that all structures discussed here correspond to classical optimal structures on the PES of each cluster, and quantum zero-point energy effects are expected to cause significant changes in their stability, especially for those that are close in energy.

Fig. 7 Schematic representation of optimal low-lying structures and their energetics obtained using the sum-of-potential He_n=1−6K₂⁺ PESs.

Finally, in Fig. 8 we report the classical binding, D_e, and evaporation energies, D_n, given as the difference between the total energies E_n−1 and E_n of adjacent size clusters (with n > 1) in their most stable (global minima, see inset plots) configurations, D_n = −E_n + E_n−1, for each He_nK₂⁺ cluster, as a function of the He atoms number n (cluster size), and compare them with experimental values of He_nK₂⁺ ion yields, as determined recently by high-resolution mass spectroscopy.³⁰ One can see that the D_e values drop down rapidly, with a rate of increase in the binding of 2.1D_e¹, 3.4D_e¹, 4.5D_e¹, 5.6D_e¹ and 6.9D_e¹ with D_e¹ = 36.36 cm⁻¹ for the n = 1 cluster. Above n = 4 the rate of energy acquisition is higher, with cluster structures marked by the formation of the He₃ motifs at the ends of the K₂⁺ cation. Distinct minor anomalies in the ion yield with the increase in the size n of the He_nK₂⁺ clusters have been observed,³⁰ and are most likely caused by anomalies in evaporation energies. As we can see the patterns of experimental ion yield values and theoretical evaporation energies are similar for such small size clusters, especially for the n = 3 and 4 cases, where the energy separation between the global potential minimum and low-lying local minima is large (more than 10 cm⁻¹). This suggests that such classical minima structures can be classified as bound quantum states once zero-point corrections are included.

Fig. 8 Binding, D_e¹ and evaporation, D_n, energies for the global potential minima (see inset plots) as a function of the cluster size n. Experimental ion yield values³⁰ are also shown for comparison reasons.

2 Summary and conclusions

The present work is focused on the study of intermolecular interactions between the K₂⁺ dimeric cations and one or few He atoms. Such small size species are theoretically tractable by high-level ab initio computations, and reliable interaction forces could be obtained. Thus, as a first step the He–K₂⁺ PES has been constructed from MRCI+Q/ECP10MDD/AV6Z data using the RKHS interpolation scheme. The strength of intermolecular interaction depends of the charge and size of the impurity, and in contrast to lighter alkali-cation dimers, we found that the He–K₂⁺ interaction shows a high multireference character, especially at configurations with θ < 60°. The ZPE and all excited vibrational states of the trimer were computed from nuclear quantum calculations. We found a high ZPE value that counts up to 48.6% of its potential well-depth, although quantum effects do not alter the stability of the classical potential minima. The most stable conformer is linear, while higher vibrational bound states exhibit large amplitude motions, and the states lying at energies above the T-shaped barrier give rise to double structures.

The sum-of-potentials approach, based on the ab initio three-body MRCI+Q HeK₂ and two-body He–He terms, has been checked and found to provide an accurate and realistic representation of higher-order He_nK₂⁺ clusters. Small size He_nK₂⁺ clusters can be used as model microsolutions, and accurate knowledge of their potential can provide a very useful initial information on structures and solvation shells. The most stable potential structures for each cluster were obtained via an evolutionary programming algorithm energy optimizations. New insights into their energetics highlight the competition between weak He–K₂⁺ and He–He intermolecular interactions, which control the stabilization of low-lying energy conformers. We clearly see that during the initial cluster growth with n > 2 formation of He n motifs, such as triangles, rhombs or ring-like structures are more favored, indicating a selective growth of these clusters during the microsolvation process. The existence of such localized compact structures should also be expected to influence short-time solute–solvent dynamics, and related with slow mobility of ions observed in He-droplets. Furthermore, by comparing binding and evaporation energies with experimental ion yield values as a function of size n, we identify some minor anomalies for small He_nK₂⁺ clusters. Such classical description of the energetics and spatial arrangements of the clustering process should be extended to explicitly include thermal and nuclear quantum effects. Although thermal effects at typically low-T in He droplets, are not expected to alter the classical solvation behavior of K₂⁺ in He clusters, quantum effects may affect the microscopic behavior of such doped clusters, especially in cases where the global and local potential minima are close in energy. The anisotropy of the underlying PES is reflected in the solvation structure, and such information could be used for describing the transition from a molecular aggregate to a dissolved molecule.

The computational cost for nuclear full quantum treatments is becoming really high on moving to larger dimensional systems. However, benchmark studies for small species could serve for testing other approaches to deal with higher dimensional systems with an increasing number of particles. Currently, path-integral molecular dynamics (PIMD) simulations are in progress for such He n-nanodroplet model systems for studying microsolvation quantum and thermal effects for modeling short-time solute–solvent dynamics.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors thank Prof. Pablo Villarreal and Prof. Hamid Berriche for useful discussions. The authors thank the “Centro de Calculo del IFF/SGAI-CSIC” and the CESGA-Supercomputing centre for allocation of computer time. This work has been supported by the MINECO grant No. FIS2017-83157-P and Communidad de Madrid grant Ref: IND2017/AMB-7696, and COST Action CA18212(MD-GAS). N. A. acknowledges financial support from the Ministry of Higher Education and Scientific Research of Tunisia.

References

1. P. Jena and Q. Sun, Chem. Rev., 2018, 118, 5755–5870.
2. G. W. F. Drake, Springer Handbook of Atomic, Molecular, and Optical Physics, Springer-Verlag, New York, New York, 2006.
3. J. P. Toennies, A. F. Vilesov and K. B. Whaley, Phys. Today, 2001, 54, 31–37.
4. J. P. Toennies and A. F. Vilesov, Angew. Chem., Int. Ed., 2004, 43, 2622–2648.
5. J. Tang and A. R. W. McKellar, J. Chem. Phys., 2003, 119, 5467–5477.
6. Y. Xu, W. Jäger, J. Tang and A. R. W. McKellar, Phys. Rev. Lett., 2003, 91, 163401.
7. C. Di Paola, F. A. Gianturco, D. López-Durán, M. P. de Lara-Castells, G. Delgado-Barrio, P. Villarreal and J. Jellinek, ChemPhysChem, 2005, 6, 1348–1356.
8. M. Y. Choi, G. E. Douberly, T. M. Falconer, W. K. Lewis, C. M. Lindsay, J. M. Merritt, P. L. Stiles and R. E. Miller, Int. Rev. Phys. Chem., 2006, 25, 15–75.
9. M. P. de Lara-Castells, R. Prosmiti, G. Delgado-Barrio, D. López-Durán, P. Villarreal, F. A. Gianturco and J. Jellinek, Phys. Rev. A: At., Mol., Opt. Phys., 2006, 74, 053201.
10. F. Marinetti, L. Uranga-Piña, E. Coccia, D. López-Durán, E. Bodo and F. A. Gianturco, J. Phys. Chem. A, 2007, 111, 12289–12294.
11. R. E. Zillich and K. B. Whaley, J. Phys. Chem. A, 2007, 111, 7489–7498.
12. D. Bonhommeau, M. Lewerenz and N. Halberstadt, J. Chem. Phys., 2008, 128, 054302.
13. J. A. Ramilowski, D. Mikosz, A. A. Farrelly, J. L. C. Fajín and B. Fernández, J. Phys. Chem. A, 2007, 111, 12275–12288.
14. R. Prosmiti, G. Delgado-Barrio, P. Villarreal, E. Yurtsever, E. Coccia and F. A. Gianturco, J. Phys. Chem. A, 2009, 113, 14718–14729.
15. C. Callegari and W. E. Ernst, Handbook of High-Resolution Spectroscopy, Hoboken, NJ, 2011, pp. 1551–1594.
16. G. Guillon, A. Zanchet, M. Leino, A. Viel and R. E. Zillich, J. Phys. Chem. A, 2011, 115, 6918–6926.
17. Y. Ajili, K. Hammami, N. E. Jaidane, M. Lanza, Y. N. Kalugina, F. Lique and M. Hochlaf, Phys. Chem. Chem. Phys., 2013, 15, 10062–10070.
18. R. Pérez de Tudela, P. Barragán, A. Valdés and R. Prosmiti, J. Phys. Chem. A, 2014, 118, 6492–6500.
19. P. Villarreal, R. Rodríguez-Cantano, T. González-Lezana, R. Prosmiti, G. Delgado-Barrio and F. A. Gianturco, J. Phys. Chem. A, 2015, 119, 11574–11582.
20. R. Rodríguez-Cantano, T. González-Lezana and P. Villarreal, Int. Rev. Phys. Chem., 2016, 35, 37–68.
21. F. Coppens, J. von Vangerow, M. Barranco, N. Halberstadt, F. Stienkemeier, M. Pi and M. Mudrich, Phys. Chem. Chem. Phys., 2018, 20, 9309–9320.
22. D. Verma, R. M. P. Tanyag, S. M. O. O'Connell and A. F. Vilesov, Adv. Phys.: X, 2019, 4, 1553569.
23. K. Nauta and R. E. Miller, Science, 2000, 287, 293–295.
24. S. Moroni and S. Baroni, Comput. Phys. Commun., 2005, 169, 404–407.
25. F. Stienkemeier and K. K. Lehmann, J. Phys. B: At., Mol. Opt. Phys., 2006, 39, R127–R166.
26. M. Gatchell, P. Martini, F. Laimer, M. Goulart, F. Calvo and P. Scheier, Faraday Discuss., 2019, 217, 276–289.
27. C. Di Paola, F. Sebastianelli, E. Bodo, F. A. Gianturco and M. Yurtsever, J. Chem. Theory Comput., 2005, 1, 1045–1054.
28. S. Müller, M. Mudrich and F. Stienkemeier, J. Chem. Phys., 2009, 131, 044319.
29. M. Theisen, F. Lackner and W. E. Ernst, Phys. Chem. Chem. Phys., 2010, 12, 14861–14863.
30. L. An der Lan, P. Bartl, C. Leidlmair, R. Jochum, S. Denifl, O. Echt and P. Scheier, Chem. – Eur. J., 2012, 18, 4411–4418.
31. M. Rastogi, C. Leidlmair, L. An der Lan, J. Ortiz de Zárate, R. Pérez de Tudela, M. Bartolomei, M. I. Hernández, J. Campos-Martínez, T. González-Lezana, J. Hernández-Rojas, J. Bretón, P. Scheier and M. Gatchell, Phys. Chem. Chem. Phys., 2018, 20, 25569–25576.
32. A. Castillo-García, T. González-Lezana, G. Delgado-Barrio and P. Villarreal, Eur. Phys. J. D, 2018, 72, 102.
33. M. Blancafort-Jorquera, A. Vilá and M. González, Phys. Chem. Chem. Phys., 2018, 29737–29753.
34. P. Stipanović and L. V. Markić, J. Phys. B: At., Mol. Opt. Phys., 2018, 51, 155101.
35. G. Auböck, J. Nagl, C. Callegari and W. E. Ernst, J. Phys. Chem. A, 2007, 111, 7404–7410.
36. E. Bodo, F. A. Gianturco, E. Yurtsever and M. Yurtsever, Mol. Phys., 2005, 103, 3223–3231.
37. G. Auböck, M. Aymar, O. Dulieu and W. E. Ernst, J. Chem. Phys., 2010, 132, 054304.
38. W. C. S. Roman Krems and B. Friedrich, Cold Molecules: Theory, Experiment, Applications, CRC Press, Boca Raton, 2009.
39. R. Rodríguez-Cantano, T. González-Lezana, R. Prosmiti, G. Delgado-Barrio, P. Villarreal and J. Jellinek, J. Chem. Phys., 2015, 142, 164304.
40. E. Coccia, E. Bodo, F. Marinetti, F. A. Gianturco, E. Yildrim, M. Yurtsever and E. Yurtsever, J. Chem. Phys., 2007, 126, 124319.
41. S. Paolini, F. Ancilotto and F. Toigo, J. Chem. Phys., 2007, 126, 124317.
42. W. I. Glaberson and W. W. Johnson, J. Low Temp. Phys., 1975, 20, 313–338.
43. P. Moroshkin, V. Lebedev and A. Weis, Phys. Rev. Lett., 2009, 102, 115301.
44. M. Rossi, M. Verona, D. E. Galli and L. Reatto, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 212510.
45. D. E. Galli, D. M. Ceperley and L. Reatto, J. Phys. Chem. A, 2011, 115, 7300–7309.
46. E. Bodo, E. Yurtsever, M. Yurtsever and F. A. Gianturco, J. Chem. Phys., 2006, 124(7), 74320.
47. C. Ghanmi, H. J. Al Qarni, O. Al-Hagan and H. Berriche, AIP Conf. Proc., 2018, 1976, 020026.
48. N. Alharzali, H. Berriche, P. Villarreal and R. Prosmiti, J. Phys. Chem. A, 2019, 123, 7814–7821.
49. T. Ho and H. Rabitz, J. Chem. Phys., 1996, 104, 2584–2597.
50. M. Iwamatsu, Comput. Phys. Commun., 2001, 142, 214–218.
51. D. Fogel, Evolutionary Computation: Toward a New Philosophy of Machine Intelligence, IEEE Press-Wiley, Hoboken, New Jersey, 2006.
52. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby and M. M. Schütz, et al., MOLPRO, Version 2012.1, A Package of Ab Initio Programs, see http://www.molpro.net.
53. S. F. Boys and F. Bernardi, Mol. Phys., 1970, 19, 553–566.
54. A. Henriet, J. Phys. B: At. Mol. Phys., 1985, 18, 3085–3103.
55. I. S. Lim, P. Schwerdtfeger, B. Metz and H. Stoll, J. Chem. Phys., 2005, 122, 104103.
56. T. Leininger, A. Nicklass, W. K”uchle, H. Stoll, M. Dolg and A. Bergner, Chem. Phys. Lett., 1996, 255, 274–280.
57. K. Schuchardt, B. Didier, T. Elsethagen, L. Sun, V. Gurumoorthi, J. Chase, J. Li and T. Windus, J. Chem. Inf. Model., 2007, 47, 1045–1052.
58. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
59. D. E. Woon and T. H. Dunning, J. Chem. Phys., 1994, 100, 2975–2988.
60. D. Feller, J. Comput. Chem., 1996, 17, 1571–1586.
61. DENEB 1.30 beta: the Nanotechnology Software by Atelgraphics, 2020, https://www.atelgraphics.com.
62. C. Schwartz, Phys. Rev., 1962, 126, 1015–1019.
63. T. J. Lee and P. R. Taylor, Int. J. Quantum Chem., 1989, 36, 199–207.
64. T. J. Lee, Chem. Phys. Lett., 2003, 372, 362–367.
65. W. Jiang, N. J. DeYonker and A. K. Wilson, J. Chem. Theory Comput., 2012, 8, 460–468.
66. R. Prosmiti, C. Cunha, P. Villarreal and G. Delgado-Barrio, J. Chem. Phys., 2003, 119, 4216–4222.
67. L. Veseth, J. Phys. B: At. Mol. Phys., 1971, 4, 20–27.
68. M. Dubernet, D. Flower and J. M. Hutson, J. Chem. Phys., 1991, 94, 7602–7618.
69. R. Prosmiti, C. Cunha, P. Villarreal and G. Delgado-Barrio, J. Chem. Phys., 2002, 117, 7017–7023.
70. L. Garcia-Gutierrez, L. Delgado-Tellez, A. Valdés, R. Prosmiti, P. Villarreal and G. Delgado-Barrio, J. Phys. Chem. A, 2009, 113, 5754–5762.
71. L. Delgado-Tellez, A. Valdés, R. Prosmiti, P. Villarreal and G. Delgado-Barrio, J. Chem. Phys., 2011, 134, 214304.
72. A. Valdés, R. Prosmiti, P. Villarreal and G. Delgado-Barrio, J. Chem. Phys., 2005, 122, 044305.
73. A. Valdés, R. Prosmiti, P. Villarreal and G. Delgado-Barrio, J. Chem. Phys., 2006, 125, 014313.
74. Á. Valdés, P. Barragán, R. Pérez de Tudela, J. S. Medina, L. Delgado-Tellez and R. Prosmiti, Chem. Phys., 2012, 399, 39–45.
75. A. Valdés and R. Prosmiti, J. Phys. Chem. A, 2013, 117, 7217–7223.
76. R. A. Aziz and M. J. Slaman, J. Chem. Phys., 1991, 94, 8047–8053.
77. P. Fuentealba, H. Preuss, H. Stoll and L. Von Szentpály, Chem. Phys. Lett., 1982, 89, 418–422.
78. S. Magnier and M. Aubert-Frécon, J. Quant. Spectrosc. Radiat. Transfer, 2003, 78, 217–225.
79. A. Jraij, A. R. Allouche, S. Magnier and M. Aubert-Frecon, Can. J. Phys., 2008, 86, 1409–1415.
80. S. Leutwyler, M. Hofmann, H.-P. Harri and E. Schumacher, Chem. Phys. Lett., 1981, 77, 257–260.
81. M. Broyer, J. Chevaleyre, G. Delacretaz, S. Martin and L. Wöste, Chem. Phys. Lett., 1983, 99, 206–212.
82. W. C. Stwalley and H. Wang, J. Mol. Spectrosc., 1999, 195, 194–228.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d0cp05406b

---