# Attachment of Hydrogen Molecules to Atomic Ions (Na+, Cl−): Examination of an Adiabatic Separation of the H2 Rotational Motion

Esther García-Arroyo, José Campos-Martínez, Massimiliano Bartolomei, Marta I. Hernández, Fernando Pirani, Nadine Halberstadt

<figure><img src="image_1.png" /></figure>

► To cite this version:

Esther García-Arroyo, José Campos-Martínez, Massimiliano Bartolomei, Marta I. Hernández, Fernando Pirani, et al.. Attachment of Hydrogen Molecules to Atomic Ions (Na+, Cl−): Examination of an Adiabatic Separation of the H2 Rotational Motion. ChemPhysChem, 2023, 24 (23), pp.e202300424. 10.1002/cphc.202300424. hal-04300081

HAL Id: hal-04300081

https://hal.science/hal-04300081v1

Submitted on 26 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

Accepted Article

**Title:** Attachment of Hydrogen Molecules to Atomic Ions (Na+, Cl-):
Examination of an Adiabatic Separation of the H2 Rotational
Motion

**Authors:** Esther García-Arroyo, José Campos-Martínez, Massimiliano Bartolomei, Marta I. Hernández, Fernando Pirani, and Nadine Halberstadt

This manuscript has been accepted after peer review and appears as an
Accepted Article online prior to editing, proofing, and formal publication
of the final Version of Record (VoR). The VoR will be published online
in Early View as soon as possible and may be different to this Accepted
Article as a result of editing. Readers should obtain the VoR from the
journal website shown below when it is published to ensure accuracy of
information. The authors are responsible for the content of this Accepted
Article.

**To be cited as:** *ChemPhysChem* **2023**, e202300424

**Link to VoR:** https://doi.org/10.1002/cphc.202300424
Attachment of Hydrogen Molecules to Atomic
Ions (Na⁺, Cl⁻): Examination of an Adiabatic
Separation of the H₂ Rotational Motion

Esther García-Arroyo,[a] José Campos-Martínez,[a] Massimiliano Bartolomei,[a]
Marta I. Hernández,[*a] Fernando Pirani,[b] Nadine Halberstadt[c]

*This paper is dedicated to our colleague and friend Prof. Pablo Villarreal on his 70th birthday.*

Interactions between molecular hydrogen and ions are of in-
terest in cluster science, astrochemistry and hydrogen stor-
age. In dynamical simulations, H₂ molecules are usually
modelled as point particles, an approximation that can fail
for anisotropic interactions. Here, we apply an adiabatic
separation of the H₂ rotational motion to build effective
pseudoatom-ion potentials and in turn study the proper-
ties of (H₂)ₙNa⁺/Cl⁻ clusters. These interaction potentials
are based on high-level *ab initio* calculations and Improved
Lennard-Jones parametrizations, while the subsequent dy-
namics has been performed by quantum Monte Carlo cal-
culations. By comparisons with simulations explicitly de-
scribing the molecular rotations, it is concluded that the
present adiabatic model is very adequate. Interestingly, we
find differences in the cluster stabilities and coordination
shells depending on the spin isomer considered (*para-* or
*ortho*-H₂), especially for the anionic clusters.

**Introduction**

Molecular hydrogen (H₂) is important in several fields such as astrochemistry [1,2], cluster science [3-5], solid phases at high pressures [6,7] and interactions of the gas with porous materials [8,9]. The last topic is of much interest for the development of safe and efficient hydrogen storage platforms, recognized as a priority in hydrogen energy research [10].

An understanding of these processes can be gained from
atomistic simulations of the motion of the atomic nuclei in-
tervening in the processes and subjected to intermolecular
potentials (or force fields). Due to their light masses, the
dynamics of the H₂ nuclei must be studied within quantum

mechanics. Indeed, it is known that nuclear quantum effects
(NQE) (zero-point energy, tunneling, quantized hindered ro-
tor states, etc.) can be significant in these systems even at
room temperature [7]. One important NQE is due to the
indistinguishability of the protons (fermions) forming the
molecule, which imposes antisymmetry in the total nuclear
wavefunction [11,12]. This leads to hydrogen being a mix-
ture of two different species or spin isomers: *para*-hydrogen
(pH₂), with a singlet nuclear spin and a rotational quantum
number *j*=*even*, and *ortho*-hydrogen (oH₂), with a triplet
nuclear spin and *j*=*odd*. Since the *ortho-para* spin conver-
sion is very slow in the gas phase, both species coexist in
the gas at the proportion given by the nuclear statistics:
25% of pH₂ and 75% of oH₂, each variety showing distinct
properties [13].

For processes involving many H₂ molecules (as in cluster
or adsorption studies) it is quite appropriate to apply quan-
tum Monte Carlo methods, such as Diffusion Monte Carlo
(DMC) and Path Integral Monte Carlo (PIMC). Usually,
these molecules are described as point particles or pseu-
doatoms (PsAt) by using effective potentials obtained from
averaging the interactions over the molecule orientations,
here denoted as Av-PsAt model. This is often a good ap-
proximation for the ground rotational state of pH₂, j=0,
provided that the interactions with other partners are nearly
isotropic. However, the Av-PsAt model can fail for aniso-
tropic interactions and/or for the description of higher ro-
tational states, as the ground rotational state of oH₂, j=1,
which imposes quantized projections of the angular momen-
tum operator. This is the case of H₂-ion interactions, whose
potential energy surfaces (PESs) are anisotropic because of
dominant electrostatic and induction forces [14], hence mak-
ing H₂ behave as a hindered rotor when attached to the
ion [11]. For these situations, the H₂ rotational degrees-of-
freedom should be explicitly included in the PIMC [15] or
DMC [16] quantum simulations. The calculations, however,
become cumbersome, particularly for oH₂, because it is not
straightforward to describe the nodes of the excited state
rotational function [17].

In this work we apply an alternative pseudoatom model
to represent H₂-ion interactions. The approach is based on
an adiabatic separation of the H₂ rotational motion from
the H₂-ion radial mode, a treatment originally proposed by
Holmgren et al[¹⁸] and later applied to study various H₂-
molecule complexes[¹⁹⁻²³]. For each value of the H₂-ion dis-
tance, *R*, the quantum states of the Hamiltonian restricted
to the rotational motion are obtained, leading to an effective
potential only depending on *R* so that the H₂ molecules be-
come described as point particles. The treatment is applied

[a] E. García-Arroyo, Dr. J. Campos-Martínez, Dr. M. Bartolomei, Dr. M. I. Hernández*
Instituto de Física Fundamental, Consejo Superior de Investigaciones Científicas (IFF-CSIC), Serrano 123, 28006 Madrid, Spain
E-mail: marta@iff.csic.es

[b] Prof. F. Pirani
Dipartimento di Chimica, Biologia e Biotecnologie, Università di Perugia, via Elce di Sotto 8, 06123, Perugia, Italy

[c] Prof. N. Halberstadt
Laboratoire Collisions, Agrégats, Réactivité (LCAR), Université de Toulouse, CNRS, 31062 Toulouse, France
to study the attachment of a number n of H₂ molecules to atomic ions -i.e. ion-doped hydrogen clusters- which is carried out using the DMC method. This adiabatic PsAt (Ad-PsAt) model is assessed by means of comparisons with accurate DMC calculations that explicitly describe H₂ molecules as rigid rotors (RigRot). Interestingly, our approach allows us to describe the properties of the clusters formed either by para or by ortho hydrogen isomers, by building effective potentials correlating with the j = 0 or j = 1 rotational states, respectively. One aspect of interest for hydrogen storage applications is the determination of the ion coordination number, i.e., the number of H₂ molecules that are directly attached to the ion, forming the first shell of the cluster. Since H₂-ion binding energies are generally larger for oH₂ than for pH₂ [24,25], it can be expected that more ortho molecules will be attached to the ion.

We have selected two ions for our study, Na⁺ and Cl⁻, which interact differently with hydrogen. Indeed, the equilibrium geometry of H₂Cl⁻ is linear while that of H₂Na⁺ is T-shaped, due to the different sign of the ion charge. Also, H₂Na⁺ binding energies are larger, but the difference between ortho and para binding energies is more significant in H₂Cl⁻. These factors are expected to affect the properties of the clusters. (H₂)ₙNa⁺ clusters have been recently investigated by some of us and other coworkers in a joint experimental-computational work [26]. Cluster populations, measured by a mass spectrometer after formation inside a helium nanodroplet, showed peaks revealing specially stable cluster sizes, in agreement with the calculations. The work is here extended with the application of the Ad-PsAt model to both p/oH₂ isomers. Regarding chlorine clusters, recording of the vibrational spectra of both oH₂Cl⁻ and pH₂Cl⁻ dimers was achieved in recent experiments [27], which confirmed theoretical predictions about strong differences in the spectral features of both spin isomers [25]. While we are not aware of experiments on larger clusters, some works have addressed their electronic structure [28,29], showing the capacity of the anion for the attachment of hydrogen molecules. In this work we present a new ab initio-based PES for H₂-Cl⁻ and report the properties of (H₂)ₙCl⁻ clusters as obtained within all the models previously mentioned.

The paper is organized as follows. Next section describes the theoretical models and the related interaction potentials used. This is followed by a Results and Discussion section, and finally, a Conclusion. Details about the pairwise interactions and the DMC calculations are provided in the Supporting Information (SI).

## Theoretical Models and Interaction Potentials

The total interaction potential or PES here adopted for a (H₂)ₙM cluster is written as

$$V_{\text{TOT}} = \sum_{i=1}^{n} V_{\text{H}_2(i)\text{-M}} + \sum_{i<k}^{n} V_{\text{H}_2(i)\text{-H}_2(k)} \quad (1)$$

where V<sub>H₂(i)M</sub> and V<sub>H₂(i)H₂(k)</sub> indicate pairwise potentials between the i-th molecule and the ion and between a couple of H₂ molecules i and k, respectively. In this section we present various models of these pairwise interactions for clusters formed with M=Na⁺ and Cl⁻ dopants.

Since most of the discussion below is focused on the M-H₂ interactions, for which we apply a new adiabatic approach,

here we briefly report the models adopted for the H₂-H₂ interaction. On the one hand, within the following RigRot model, we have used an analytical potential depending on the intermolecular distance as well as the relative orientations of H₂-H₂. On the other hand, a simpler analytical form, only depending on the intermolecular distance and obtained by averaging over all the relative orientations, is used in combination with the two pseudoatom models (Ad-PsAt and Ad-PsAt) considered for the H₂-ion interaction and reported in detail below. Explicit expressions of these potentials and values of the related parameters are given in the SI. It should be possible to develop an adiabatic separation of the angular degrees of freedom of the H₂-H₂ interaction as well and, in this way, to achieve a more realistic representation of the interactions between para or ortho molecules. In this work, however, we wish to keep the description of these interactions as simple as possible and leave more refined approaches for future investigations.

## Rigid Rotor (RigRot) Model and Interaction Potentials

In the RigRot model, the position of the center of mass as well as the orientation of the H₂ molecules are specified, while the intramolecular distance of H₂ is kept fixed to a value consistent with the vibrational ground state. This is the most complete model presented here so it is taken as a reference to test the subsequent PsAt approximations.

In this approach, V<sub>H₂-M</sub> terms depend on the distance R between M and the center of mass of H₂ as well as on the angle θ formed between the H₂ axis and the vector joining M and the H₂ center of mass (Jacobi coordinates). The H₂-M interaction is dominated by electrostatic, induction and van der Waals contributions and it is analytically represented as a sum of an atom-bond Improved Lennard-Jones function [30,31] and the electrostatic energy between the ion charge and the H₂ quadrupole moment. These parameters have been optimized from comparisons with high level ab initio calculations, which, in the case of H₂-Na⁺, have been reported in Ref. [24]. Details on new ab initio calculations and optimized parameters for the H₂-Cl⁻ interaction are reported in the SI. It is worth noting that the new H₂-Cl⁻ interaction potential is in good agreement with previous ab initio calculations of the complex [32,33]. An extended atom-bond Improved Lennard-Jones formulation plus a electrostatic contribution have been used to represent the H₂-H₂ interaction [34]. For the clusters doped with Na⁺, a further three-body (3B) interaction term has been added, consisting on the interaction between the dipoles that the ion induces on the hydrogen molecules [26]. For (H₂)ₙCl⁻, we have concluded, based on ab initio calculations of the (H₂)₂Cl⁻ complex, that it is not necessary to include this 3B term since higher intermolecular distances and weaker interaction components are involved.

Explicit expressions for all the mentioned interactions are gathered in the SI, with the involved parameters being given in Table S1 therein.

Profiles of the H₂-Na⁺ and H₂-Cl⁻ PESs are presented in Fig. 1 for the perpendicular ((θ=π/2) as well as the parallel (θ=0) orientations of the diatomics. It can be seen that the energy minimum corresponds to T- and L-shaped geometries for Na⁺ and Cl⁻, respectively. This is due to a dominant charge-quadrupole interaction that favors perpendicular or parallel orientations for positive or negative ion
<figure><img src="image_6.png" /><figcaption>**Figure 1.** H₂-Na⁺ (upper panel) and H₂-Cl⁻ (lower panel) interaction potentials (in meV) vs. H₂-ion distance (in Å). Dashed and dot-dashed black curves: RigRot PES for parallel and perpendicular orientations, respectively. Orange curves: pseudoatom potentials after averaging the RigRot PES (Av-PsAt). Red and blue curves: adiabatic pseudoatom potentials (Ad-PsAt), correlating with the *para* and *ortho* ground states of H₂, respectively. See text for discussion.</figcaption></figure>

charges, respectively. It is also noticeable the anisotropy of the interactions, i.e., the profiles of the potential curves differ very much between the parallel and the perpendicular orientations.

## Orientationally Averaged Pseudoatom Model (Av-PsAt)

In the Av-PsAt model, the involved interactions are obtained by carrying out averages of the RigRot potentials from Eq. 1 over all the molecular orientations. These terms, only depending on the distances between the interacting partners, are well represented by ILJ functions, as indicated in the SI and Table S1 therein. It is worth noting that the charge-quadrupole interaction cancels upon averaging. The resulting Av-PsAt potentials are depicted in Fig. 1 (orange lines). It can be seen that these interactions are more similar to the perpendicular potential curves, due to the larger weight that configurations close to the perpendicular one have on the final average.

## Orientationally Adiabatic Pseudoatom Model (Ad-PsAt) for p/oH₂

The present Ad-PsAt model is identical to the previous Av-PsAt one except for the representation of the H₂-ion interaction, which is modified here resorting to an adiabatic separation of the angular and radial degrees of freedom of the dimer.

The Hamiltonian for the H₂-ion interaction can be written as, within a rigid-rotor approximation for H₂,

$$ \mathcal{H} = -\frac{\hbar^2}{2\mu} \frac{\partial^2}{\partial R^2} + \mathcal{H}_{(\theta,R)}, \quad (2) $$

where the first term is the kinetic energy for the radial degree of freedom ($\mu$ being the reduced mass of the H₂-ion complex) and

$$ \mathcal{H}_{(\theta,R)} = \frac{\hat{l}^2}{2\mu R^2} + B \hat{j}^2 + V(\theta, R), \quad (3) $$

where $\hat{l}$ is the orbital angular momentum (associated to the vector $\hat{R}$), $\hat{j}$, the rotational angular momentum of H₂, $B=7.355$ meV, the H₂ rotational constant, and finally, $V(\theta, R)$ is the RigRot PES. In the present approach, the Hamiltonian of Eq. 3 is diagonalized with respect to the angular degrees of freedom at each value of the intermolecular distance $R$. The resulting eigenvalues, $V^{ad}(R)$, become effective potentials for the H₂-ion interaction, where the molecule is then treated as a point particle. Basis sets of free rotors are used to compute the Hamiltonian matrix elements, which are defined in a body-fixed reference system as detailed in Ref. [35] (in particular, see Eq. (8) therein). In the calculations, these bases involve rotational functions from $j_{min} = 0$ (1) to $j_{max} = 12$ (11) for the *para* (ortho) states. These bases have well-defined quantum numbers of the total angular momentum, $J$, the inversion parity $\epsilon_i = \pm 1$, and the parity of the H₂ rotational states, $\epsilon_j$ ($\pm 1$ for *para* and *ortho* states, respectively). Therefore, the resulting effective potentials have well defined values of $J$, $\epsilon_i$ and $\epsilon_j$. In particular, potentials representing the interaction with pH₂ and oH₂, $V_p^{ad}(R)$ and $V_o^{ad}(R)$, are chosen from the lowest-energy curves asymptotically correlating with the energy of the $j=0$ and $j=1$ H₂ rotational states, respectively.

The effective potentials for pH₂-Na⁺ and oH₂-Na⁺ are presented in the upper panel of Fig. 1. They correspond to the $J$, $(\epsilon_i, \epsilon_j)$ symmetry blocks $J=0$, (+1,+1) and $J=1$, (+1,-1), respectively. It can be pointed out that the *ortho* curve is almost degenerate with the one of $J=1$ with symmetry (-1, -1) (both curves are within less than 0.01 meV near the equilibrium distance). For the *ortho* states, the projection of the H₂ angular momentum $\hat{j}$ on the H₂-Na⁺ axis suggests that the molecule rotates in a plane perpendicular to that axis (“helicopter” motion). Other adiabatic potentials correlating with oH₂ ($j=1$) are higher in energy (about 50 meV above the $J=1$, ($\pm 1, -1$) curves at the equilibrium distance). A way to test the validity of this adiabatic approximation is to compute the ground state energies of the *para* and *ortho* isomers by solving the equation

$$ \left[ -\frac{\hbar^2}{2\mu} \frac{\partial^2}{\partial R^2} + V_{p/o}^{ad}(R) \right] \Phi_{p/o}^{0}(R) = E_{p/o}^{0} \Phi_{p/o}^{0}(R), \quad (4) $$

and compare them with the energies resulting from a fully coupled calculation. The result is quite satisfactory, as we obtain $E_p^0 = -111.98$ meV and $E_o^0 = -103.01$ meV, in very good agreement with the results of close-coupling computations [24], -111.86 and -102.94 meV, respectively. In Fig. 1 (upper panel), it can be noticed that the $V_p^{ad}(R)$ potential is more attractive than the spherically averaged one (Av-PsAt). It is also worth observing that the separation between the *ortho* and *para* curves ($V_o^{ad}(R) - V_p^{ad}(R)$) decreases as the partners approach. Indeed, the splitting between the *ortho* and *para* states is smaller when the molecule
is bound to the cation ($E_o^0 - E_p^0 = 8.97$ meV) than that in the free molecule ($2B= 14.71$ meV). As a consequence, $oH_2Na^+$ is a more stable complex than $pH_2Na^+$, since energy required to break it (dissociation energy) is larger ($D_o(o)= 103.01+2B= 117.72$ meV vs $D_o(p)= 111.98$ meV).

Adiabatic potentials for the H₂-Cl⁻ complex are shown in the lower panel of Fig. 1. The *para* potential (in red) has $J=0$ with (+1,+1) symmetry, whereas the *ortho* one also has $J=0$ but with (+1,-1) symmetry. In this *ortho* state, the projection of the H₂ angular momentum over the H₂-Cl⁻ axis is zero, implying that the molecule is spinning in a plane parallel to that axis (“ferris wheel” or “cart-wheel” motion). As in the previous system, other potentials correlating with $oH_2(j = 1)$ are higher in energy (42 meV above the $J=0$, (+1,-1) curve at the equilibrium distance). Features mentioned for the H₂-Na⁺ lowest effective potentials also occur for the present ones, though in a more pronounced manner: the $V_p^{ad}(R)$ potential is much more attractive than the Av-PsAt one and the separation between the $V_p^{ad}(R)$ and $V_o^{ad}(R)$ curves becomes very small (about 1 meV) near the minimum of the interaction. Indeed, we obtain very close $pH_2Cl⁻$ and $oH_2Cl⁻$ ground state levels $E_p^0 = -50.79$ meV and $E_o^0 = -49.45$ meV, leading to significantly different dissociation energies $D_o(p)= 50.79$ meV vs $D_o(o)= 49.45 + 2B= 64.16$ meV, respectively. It is worth noting that these dissociation energies are in good agreement with those obtained, using different PESs, by Alexander[33] ($D_o(p)= 51.79$ meV, $D_o(o)= 65.36$ meV) and by Buchachenko et al.[32] ($D_o(p)= 51.06$ meV, $D_o(o)= 49.89 + 2B= 64.60$ meV).

## Results and Discussion: Cluster Energies and Structures

Within the three theoretical models presented in the previous section, ground state energies, $E_n$, and structures of the $(p/oH_2)_nNa^+/Cl^-$ clusters have been obtained by means of DMC calculations, in an implementation provided by Buch and Sandler[16,36,37] where the hydrogen molecules can be described either as rigid rotors (RigRot model) or as pseudоatoms (Av- or Ad-PsAt models). Details for these DMC computations are provided in the SI.

Our goal is to test the validity of the Ad-PsAt model by comparison of the cluster properties with those obtained from the RigRot model, which gives the most realistic results. At the same time, improvements brought by the adiabatic effective potentials, as compared with the Av-PsAt ones, are discussed. As for clusters formed by *ortho* molecules, there is no reference available to which compare the Ad-PsAt results (except for the *n* = 1 case, already discussed), so present results should be taken as reasonable predictions of the behavior of these clusters.

In Fig. 2 we present the energy per molecule, $E_n/n$, of Na⁺- and Cl⁻-doped clusters for the various models considered. For the clusters containing Na⁺ (upper panel), it can be seen that the Ad-PsAt model for $pH_2$ molecules provides energies indistinguishable to those of the RigRot model, whereas orientationally averaged potentials underestimate the absolute values of the energies by about 16-20%. In addition, the Ad-PsAt model for *ortho* isomers predicts that the binding energies of the $(oH_2)_nNa^+$ clusters are about 5-8% larger than those of the *para* clusters. For the aggregates doped with Cl⁻ (lower panel), the Av-PsAt energies involve quite large relative errors (30-40%) with respect to

<figure><img src="image_7.png" /><figcaption>**Figure 2.** Energy per molecule, $E_n/n$, of $(H_2)_nNa^+$ (upper panel) and $(H_2)_nCl^-$ (lower panel) vs. *n*, for the different models considered here: RigRot (black curve with open circles), Av-PsAt (orange with filled diamonds), and Ad-PsAt, for hydrogen molecules in *para* and *ortho* states (red and blue with filled diamonds, respectively).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 3.** Evaporation energies, $E_n^{evap} = E_{n-1} - E_n$ (meV), of $(H_2)_nNa^+$ (upper panel) and $(H_2)_nCl^-$ (lower panel) vs. number of hydrogen molecules, *n*, for the different models considered in this work (labelled as in Fig. 2). See text for discussion.</figcaption></figure>

the RigRot values. This is due to the large anisotropy of the H₂-Cl⁻ PES (Fig. 1). In contrast, the Ad-PsAt-*para* energies do compare quite well with the RigRot ones, despite a tendency of the adiabatic approach to slightly overestimate the binding energies. Moreover, $(oH_2)_nCl^-$ clusters are es-
<figure><img src="image_9.png" /><figcaption>**Figure 4.** Upper panel: Radial density, $D(R)$ (Å⁻³), of the $(H_2)_{11}$Na⁺ cluster as a function of the distance to the ion, $R$ (Å), for the different models studied here: RigRot (black curve), Av-PsAt (orange curve), and Ad-PsAt, for pH₂ and oH₂ (red and blue curves, respectively). Lower panel: For the same cluster, number of H₂ molecules inside a given H₂-Na⁺ distance, $N(R)$, vs. $R$ (Å), for the same models considered. See text for discussion.</figcaption></figure>

timated to be much more stable than those formed by the *para* species (binding energies 25-35% larger).

Cluster evaporation energies, defined as the energy required to remove the most weakly bound molecule from a cluster of size $n$, $E_{n}^{\text{evap}} = E_{n-1} - E_{n}$, are shown in Fig. 3. This quantity allows us to get more insight into the stability of the clusters, even serving in some cases for comparisons with measurements of cluster abundances as functions of their size[³⁸]. It can be seen that the evaporation energies of the Na⁺-doped clusters (upper panel) decrease with $n$ already for the smallest clusters, despite having rather large values. This is due to the role of the (repulsive) 3B interaction as well as the increasing repulsion among the hydrogen molecules when their number increases[²⁶] or, in other words, the packing difficulties of molecules surrounding the ion with a “helicopter motion”. Nevertheless, clusters with six and eight H₂ molecules show a special stability with respect to neighboring sizes and indeed it has been experimentally found that they are more abundant[²⁶]. Beyond these sizes, evaporation energies drop significantly. Regarding the performance of the different approaches, notice again the perfect agreement between the Ad-PsAt-*para* and RigRot evaporation energies and the similar behavior of this quantity for the *ortho* clusters. Taking into account the good agreement of the RigRot approach with the experimental measurements[²⁶], this result gives us confidence on the reliability of the Ad-PsAt approach. Moreover, the Av-PsAt evaporation energies differ from the reference values for small cluster sizes but a good agreement is attained for $n> 6$, indicating that orientational effects are of little importance for the larger clusters.

With regard to the $(H_2)_nCl^-$ clusters (lower panel of Fig.

<figure><img src="image_10.png" /><figcaption>**Figure 5.** Same as in Fig. 4, for the $(H_2)_{15}$Cl⁻ cluster.</figcaption></figure>

3), it is noticeable that the evaporation energies remain almost constant up to $n=8$, suggesting that all H₂ molecules share equivalent bondings with the anion. Indeed, it is reasonable to expect that repulsion among the molecules is not significant for the smaller clusters, as the molecules should tend to align colinearly to the anion. For $n > 8$, evaporation energies decrease but increase again at $n = 12$, where the cluster exhibits a special stability. It can be seen that the Ad-PsAt-*para* evaporation energies agree rather well with the RigRot ones, despite some slight overestimation for $n > 8$. On the other hand, the Av-PsAt model poorly reproduces the behavior of the accurate evaporation energies at almost all sizes. Finally, notice the large values of the $(oH_2)_nCl^-$ evaporation energies in comparison to the ones for the *para* isomer, for almost all sizes. This difference suggests that an enrichment of clusters composed by the *ortho* species is possible in some conditions, as found for the H₂Cl⁻ dimer[¹¹,²⁷].

The drop in the evaporation energies for cluster sizes larger than 8 (Na⁺) and 12 (Cl⁻) suggests that a first shell of molecules directly surrounding the ion has been saturated and that extra molecules are placed more distantly to the ion, where escaping from the cluster requires less energy. To study the cluster structures in more detail, we have obtained the probability density of the cluster state as a function of the distance to the ion, $D(R)$, normalized to the number of molecules,

$$ \int_{0}^{\infty} D(R')R'^{2}dR' = n. \qquad (5) $$

This distribution is shown in the upper panel of Fig. 4 for $(H_2)_{11}$Na⁺ and the different models considered. It can be seen that the Ad-PsAt distributions, both for the *para* and *ortho* isomers, are almost indistinguishable from the RigRot one, showing a first shell between 2 and 3.5 Å and a second one, between 4 and 5.5 Å. Within the Av-PsAt approximation, the first shell is shifted to slightly larger distances. On the other hand, the quantity obtained by integrating $D(R')$ just up to a given distance $R$,
<table><caption>**Table 1.** Coordination number of Na and Cl⁻ (sizes of the first shell), within different models and methods ("Classical" refers to the structure of the minimum of the PES).</caption><thead><tr><th>Model</th><th>Method</th><th>Na<sup>+</sup></th><th>Cl<sup>-</sup></th></tr></thead><tbody><tr><td>RigRot</td><td>Classical</td><td>9</td><td>14</td></tr><tr><td>RigRot</td><td>DMC</td><td>8-9</td><td>12</td></tr><tr><td>Av-PsAt</td><td>DMC</td><td>8-9</td><td>15</td></tr><tr><td>Ad-PsAt-para</td><td>DMC</td><td>8-9</td><td>12</td></tr><tr><td>Ad-PsAt-ortho</td><td>DMC</td><td>8-9</td><td>14</td></tr></tbody></table>

$$ \mathcal{N}(R) = \int_{0}^{R} \mathcal{D}(R') R'^{2} dR', \qquad (6) $$

indicates the number of molecules that are inside a sphere of radius *R* and hence, the sizes of the shells. Indeed, as seen in the lower panel of Fig. 4, *N*(*R*) rises rapidly with *R* up to a distance between the shells (3.5-4.5 Å) where it varies very slowly, suggesting that the first shell has between eight and nine molecules (interpretation of this result follows below). For larger distances, *N*(*R*) increases more rapidly again until it asymptotically reaches 11, the size of the cluster. Notice that the behavior of this function is very similar for all of the models under consideration, leading to the same conclusions about the cluster structure.

Analogous probability distributions for the (H₂)₁₅Cl⁻ cluster are shown in Fig. 5, where a starker difference among the models is clearly noticed. First, distributions of Ad-PsAt model for the *para* species are in good agreement with the RigRot ones, both approaches showing that the first shell of this aggregate has 12 molecules. Notice that, unlike the sodium-doped clusters, the probability density becomes negligible in the region between the shells, suggesting a solid-like structure of the first shell. On the contrary, the Av-PsAt probability density is largely different (upper panel) and the corresponding *N*(*R*) function (lower panel) indicates absence of a second shell. Finally, the Ad-PsAt model for the *ortho* species predicts a first shell with 14 molecules, two more than the coordination number of the *para* isomer.

We have also analyzed the radial distributions of the clusters with different number of molecules and confirmed that the sizes of the first shell coincides with the ones found from the analysis above. The results are gathered in Table 1. In its first row, the corresponding numbers associated to the structure of the absolute minimum of the RigRot PES are also shown. As an example, these classical structures are depicted in Fig. 6 for the clusters of *n* = 11 (Na⁺) and *n* = 15 (Cl⁻) molecules. For the sodium dopant (left panel), it can be seen that the first shell has 9 molecules: 8 of them form a deformed squared antiprism[39], while the additional H₂ locates on top of one of the squared face of the mentioned polyhedron. For chlorine (right panel), 14 molecules coordinate the anion, while the 15th one places at much larger distance to the anion. These results are compared with the RigRot/DMC ones (2nd row in Table 1) to get insight into the role of quantum effects (zero-point energy and anharmonicity) on modifying the predictions uniquely based on the topography of the PESs. On the one hand, quantum effects in sodium clusters cause a delocalization of the 9th molecule of the classical first shell: as seen in Fig. 4, there is a non-negligible probability in the inter-shell region, therefore, a non integer number, between 8

<figure><img src="image_11.png" /><figcaption>**Figure 6.** Classical structures (minima of the RigRot PES) of the clusters (H₂)₁₁Na⁺ and (H₂)₁₅Cl⁻. Na and Cl⁻ ions are shown in red and purple, respectively, and all H₂ are in grey except those assigned to a second coordination shell, which are given in black.</figcaption></figure>

and 9, is assigned to the size of the quantum shell. On the other hand, the chlorine coordination number becomes 12 at the quantum level (RigRot and Ad-PsAt-*para*/DMC). This can be understood from the zero-point energy motion of the wavefunction, which extends to configurations where the H₂-H₂ interaction becomes repulsive, so the first shell becomes saturated with a number of molecules smaller than the classical estimation.

Although our main goal here has been the examination of different approaches describing nuclear quantum effects in the studied clusters, it is worth comparing our (H₂)ₙCl⁻ RigRot PES with the electronic structure study of the same system by Della and Suresh[29]. As detailed in Section 1.c) of the SI, the minimum energies of our PESs and the DFT counterparts[29] can be considered to be in a good agreement, taking into account the various differences between both approaches (pairwise vs. supermolecular, level of calculations and optimization). Cluster structures also agree up to *n* = 14. The main difference is that, in the DFT study, the H₂ molecules remain in the first shell up to *n* = 16, while a second shell starts forming for *n* > 14 in our calculations. We believe that the discrepancy in the coordination number mainly arises from the resulting DFT interactions that, with respect to those adopted in our treatment, provide stronger binding energies, therefore able to better compensate for the increasing H₂-H₂ repulsion among the molecules in the same shell (which is eventually the cause for the formation of a second shell).

## Conclusion

In this work we have proposed effective H₂-ion interaction potentials where the molecule is treated as a point particle (pseudoatom), based on an adiabatic separation of the H₂-ion radial and H₂ rotational degrees of freedom. This has allowed us to build distinct interactions for the *para* and *ortho* spin isomers of hydrogen, by selecting adiabatic curves correlating with the molecule in the ground rotational state (pH₂) or the first excited one (oH₂). The model is applied to the examination of the stability and structure of hydrogen molecules attached to those impurities. Such a study is done by means of quantum Diffusion Monte Carlo calculations and using pairwise Improved Lennard Jones force fields optimized from accurate *ab initio* calculations. For clusters containing pH₂, we have attained an excellent agreement with analogous Monte Carlo calculations where the H₂ rotational degrees of freedom are fully included in the sim-
ulations. Also, a significant improvement with respect to pseudoatom models based on simpler orientational averages has been obtained. This is because simpler pseudoatom models do not capture the anisotropy of the H₂-ion interactions, which indeed plays an important role in the stability and structure of the clusters. Moreover, we have been able to analyze differences in the clusters features when the *para* species is substituted by the *ortho* one. Clusters formed with *oH₂* molecules are found to be more stable and tend to have larger coordination numbers (size of the first shell) than the ones composed by *pH₂*, particularly in the case of the anionic dopant.

We believe that the importance of nuclear quantum effects in the attachment of hydrogen molecules to ions has been demonstrated, specially, the quantization of the H₂ rotational states that are hindered by the anisotropic interaction with the ion. It has been seen that, upon including the nuclear dynamics, the size of the first coordination shell tends to become reduced with respect to predictions solely based on inspection of the equilibrium structure of the potential energy surface. Moreover, shells formed by *oH₂* may accept more molecules than those composed of the *para* species, due to the larger dissociation energy of *oH₂*-ion dimers. These are interesting features to take into account in the modelling of systems for applications in hydrogen storage.

## Supporting Information

The authors have cited additional references within the Supporting Information. [40–49]

## Acknowledgements

E. G.-A. is in the Doctoral School Universidad Autónoma de Madrid (Programme in Condensed Matter Physics, Nanoscience and Biophysics). The work has been funded by Spanish Agencia Estatal de Investigación projects PID2020-114957GB-I00/AEI/10.13039/501100011033 and PID2020-114654GB-I00/AEI/10.13039/501100011033. Collaboration has been supported by the CSIC under i-Link+ program LINKB20041. Allocation of computing time by CESGA (Spain) is also acknowledged.

## Conflict of Interest

There are not conflicts of interest to declare.

**Keywords:** computational chemistry • hydrogen • (ortho/para)hydrogen-ion interactions • nuclear quantum effects • hydrogen storage

## References

[1] Y. Shagam, A. Klein, W. Skomorowski, R. Yun, V. Averbukh, C. P. Koch, E. Narevicius, *Nat. Chem.* **2015**, *7*, 921.

[2] J. R. Goicoechea, O. Roncero, *A&A* **2022**, *664*, A190.

[3] G. Tejeda, J. M. Fernández, S. Montero, D. Blume, J. P. Toennies, *Phys. Rev. Lett.* **2004**, *92*, 223401.

[4] V. M. Akimov, L. I. Kolesnikova, L. Y. Rusin, M. B. Sevryuk, J. P. Toennies, *Russ. J. Phys. Chem. B* **2009**, *3*, 743.

[5] K. Prozument, B. G. Sartakov, A. F. Vilesov, *Phys. Rev. B* **2020**, *101*, 184507.

[6] X. Song, C. Liu, Q. Li, R. J. Hemley, Y. Ma, C. Chen, *Proc. Natl. Acad. Sci. U.S.A.* **2022**, *119*, e2122691119.

[7] W. Fang, J. Chen, Y. Feng, X.-Z. Li, A. Michaelides, *Int. Rev. Phys. Chem.* **2019**, *38*, 35.

[8] P. Chen, X. Wu, J. Lin, K. L. Tan, *Science* **1999**, *285*, 91.

[9] V. Jain, B. Kandasubramanian, *J. Mater. Sci.* **2020**, *55*, 1865–1903.

[10] L. Zhang, M. D. Allendorf, R. Balderas-Xicohténcatl, D. P. Broom, G. S. Fanourgakis, G. E. Froudakis, T. Gennett, K. E. Hurst, S. Ling, C. Milanese, et al., *Progress in Energy* **2022**, *4*, 042013.

[11] T. A. Grinev, A. A. Buchachenko, R. V. Krems, *Chem. Phys. Chem.* **2007**, *8*, 815.

[12] K. Hyeon-Deuk, I. Y. Chang, *Commun. Chem.* **2022**, *5*, 168.

[13] J. W. Leachman, R. T. Jacobsen, S. G. Penoncello, E. W. Lemmon, *J. Phys. Chem. Ref. Data* **2009**, *38*, 721.

[14] V. Dryza, B. Poad, E. Bieske, *Phys. Chem. Chem. Phys.* **2012**, *14*, 14954.

[15] E. G. Noya, L. M. Sesé, R. Ramírez, C. McBride, M. M. Conde, C. Vega, *Molecular Physics* **2011**, *109*, 149.

[16] V. Buch, *J. Chem. Phys.* **1992**, *97*, 726.

[17] D. Marx, M. H. Müser, *J. Phys. Condens. Matter* **1999**, *11*, R117.

[18] S. L. Holmgren, M. Waldman, W. Klemperer, *J. Chem. Phys.* **1977**, *67*, 4414.

[19] H. Li, P.-N. Roy, R. J. Le Roy, *J. Chem. Phys.* **2010**, *133*, 104305.

[20] T. Zeng, H. Li, R. J. Le Roy, P.-N. Roy, *J. Chem. Phys.* **2011**, *135*, 094304.

[21] T. Zeng, H. Li, P.-N. Roy, *J. Phys. Chem. Lett.* **2013**, *4*, 18.

[22] L. Wang, D. Xie, R. J. Le Roy, P.-N. Roy, *J. Chem. Phys.* **2013**, *139*, 034312.

[23] M. Berg, A. Accardi, B. Paulus, B. Schmidt, *J. Chem. Phys.* **2014**, *141*, 074303.

[24] M. Bartolomei, T. González-Lezana, J. Campos-Martínez, M. I. Hernández, F. Pirani, *J. Phys. Chem. A* **2019**, *123*, 8397.

[25] M. Lara-Moreno, P. Halvick, T. Stoecklin, *Phys. Chem. Chem. Phys.* **2020**, *22*, 25552.

[26] S. Kollotzek, J. Campos-Martínez, M. Bartolomei, F. Pirani, L. Tiefenthaler, M. I. Hernández, T. Lázaro, E. Zunzunegui-Bru, T. González-Lezana, J. Bretón, J. Hernández-Rojas, O. Echt, P. Scheier, *Phys. Chem. Chem. Phys.* **2023**, *25*, 462.

[27] F. Dahlmann, C. Lochmann, A. N. Marimuthu, M. Lara-Moreno, T. Stoecklin, P. Halvick, M. Raoult, O. Dulieu, R. Wild, S. Schlemmer, S. Brücken, R. Wester, *J. Chem. Phys.* **2021**, *155*, 241101.

[28] B. Nyulasi, A. Kovács, *Chem. Phys. Lett.* **2006**, *426*, 26.

[29] T. D. Della, C. H. Suresh, *Phys. Chem. Chem. Phys.* **2016**, *18*, 14588.

[30] F. Pirani, M. Alberti, A. Castro, M. Moix Teixidor, D. Cappelletti, *Chem. Phys. Lett.* **2004**, *394*, 37–44.

[31] F. Pirani, S. Brizi, L. Roncaratti, P. Casavecchia,
D. Cappelletti, F. Vecchiocattivi, *Phys. Chem. Chem. Phys.* **2008**, *10*, 5489.

[32] A. A. Buchachenko, T. A. Grinev, J. Kłos, E. J. Bieske, M. M. Szczśniak, G. Chałasiński, *J. Chem. Phys.* **2003**, *119*, 12931.

[33] M. H. Alexander, *J. Chem. Phys.* **2003**, *118*, 9637.

[34] J. Ortiz de Zárate, M. Bartolomei, T. González-Lezana, J. Campos-Martínez, M. I. Hernández, R. P. de Tudela, J. Hernández-Rojas, J. Bretón, F. Pirani, L. Kranabetter, P. Martini, M. Kuhn, F. Laimer, P. Scheier, *Phys. Chem. Chem. Phys.* **2019**, *21*, 15662.

[35] N. Halberstadt, J. A. Beswick, K. C. Janda, *J. Chem. Phys.* **1987**, *87*, 3966.

[36] P. Sandler, V. Buch, J. Sadlej, *J. Chem. Phys.* **1996**, *105*, 10387.

[37] P. Sandler, V. Buch, *General purpose QCLUSTER program for Rigid Body Diffusion Monte Carlo simulation of an arbitrary molecular cluster*, private communication **1999**.

[38] T. González-Lezana, O. Echt, M. Gatchell, M. Bartolomei, J. Campos-Martínez, P. Scheier, *Int. Rev. Phys. Chem.* **2020**, *39*, 465.

[39] Wikipedia, Square antiprism, can be found under https://en.wikipedia.org/wiki/Square_antiprism.

[40] K. Patkowski, W. Cencek, P. Jankowski, K. Szalewicz, J. Mehl, G. Garberoglio, A. H. Harvey, *J. Chem. Phys.* **2008**, *129*, 094304.

[41] A. Halkier, T. Helgaker, P. Jorgensen, W. Klopper, H. Koch, J. Olsen, A. K. Wilson, *Chem. Phys. Lett.* **1998**, *286*, 243-252.

[42] A. Halkier, T. Helgaker, P. Jorgensen, W. Klopper, A. K. Olsen, J. and Wilson, *Chem. Phys. Lett.* **1999**, *302*, 437-446.

[43] R. A. Kendall, T. H. Dunning, R. J. Harrison, *J. Chem. Phys.* **1992**, *96*, 6796.

[44] S. Boys, F. Bernardi, *Mol. Phys.* **1970**, *19*, 553-566.

[45] H.-J. Werner, P. J. Knowles, R. Lindh, F. R. Manby, M. Schutz, P. Celani, T. Korona, G. Rauhut, R. D. Amos, A. et al. Bernhardsson, MOLPRO, Version 2012.1, a Package of Ab Initio Programs **2012**.

[46] J. B. Anderson, *J. Chem. Phys.* **1975**, *63*, 1499.

[47] M. A. Suhm, R. O. Watts, *Phys. Rep.* **1991**, *204*, 293.

[48] J. K. Gregory, D. C. Clary, *Chem. Phys. Lett.* **1994**, *228*, 547.

[49] S. J. Kolmann, J. H. D'Arcy, M. J. T. Jordan, *J. Chem. Phys.* **2013**, *139*.
Entry for the Table of Contents

<figure><img src="image_12.png" /></figure>

In this computational work it is shown that hydrogen clusters doped with ions, formed by molecules in the first excited rotational state (ortho-H₂), are more stable and tend to have larger coordination numbers than clusters composed by molecules in the ground rotational state (para-H₂).