Ab Initio Potential Energy Surface for NaCl–H2 with Correct Long-Range Behavior

We report a full dimensional ab initio potential energy surface for NaCl–H2 based on precise fitting of a large data set of CCSD(T)/aug-cc-pVTZ energies. A major goal of this fit is to describe the very long-range interaction accurately. This is done in this instance via the dipole–quadrupole interaction. The NaCl dipole and the H2 quadrupole are available through previous works over a large range of internuclear distances. We use these to obtain exact effect charges on each atom. Diffusion Monte Carlo calculations are done for the ground vibrational state using the new potential.


INTRODUCTION
Molecular hydrogen is the fundamental component in the interstellar medium and is prevalent in most interstellar environments.Collisions involving molecular hydrogen interacting with itself and other rotationally and vibrationally excited molecules are the essential components of astrophysics to understand the evolution of planetary gases. 1,2Since replicating the precise conditions of astrophysical environments in the laboratory is challenging, most collision studies involving molecular hydrogen and other molecules rely on quantum scattering calculations, necessitating an accurate understanding of the relevant interaction potential energy surface (PES).In recent years, researchers have addressed the potential energy surface of various systems involving molecular hydrogen, including H 2 + H 2 , 3 SiO 2 −H 2 , 4 CN−H 2 , 5 CO−H 2 , 6 and others. 7,8Apart from these predominant molecules, the formation of molecules composed of elements with lower abundances also occurs.For instance, recently, sodium chloride has been discovered in Europa, 9 in the disk around Orion Source I, 10 as well as in other locations. 11,12Hence, the computation of a full-dimensional PES for the NaCl−H 2 system is required to understand the collisions between them.
Previously, experimental and theoretical studies have been conducted to investigate the behavior of H 2 on crystalline solid surfaces and in aqueous solutions of NaCl.For instance, Ewing, Heidberg, and others conducted multiple experiments to explore the adsorption of H 2 on annealed NaCl films and NaCl(110) using infrared (IR) spectroscopy 13−16 and He atom scattering. 17Additionally, Monte Carlo simulations and perturbation theory calculations were utilized to investigate the structure of monolayer and bilayer films of H 2 molecules adsorbed on NaCl(001) at different temperatures. 18Recently, Zhu et al. developed an accurate model of H 2 solubility in aqueous NaCl solution. 19To the best of our knowledge, no NaCl−H 2 PES exists addressing the interaction between NaCl and H 2 .
The accuracy of the PES for many atom systems, including biomolecules, condensed matter, and macromolecules, depends on the precise depiction of both short-range and longrange interactions.Therefore, accurately describing the longrange interaction potential is crucial to PES calculation 20 and understanding phenomena such as molecular collisions, 21,22 dispersion, 23 and extended charge transfer. 24,25To achieve this goal, modeling long-range interactions using electrostatic, induction, and dispersion components yields computationally efficient and precise results. 26escribing long-range interactions is a challenge for the large class of atom-centered machine learning methods that use a strict cutoff for the interaction range to represent potentials. 27,28Quoting from Behler and co-workers in 2021, "Machine learning potentials•••[that] rely on local properties••• are unable to take global changes in the electronic structure into account, which result from long-range charge transfer or different charge states." 29The authors then present ad hoc proposals to account for long-range interactions in their fourthgeneration neural network approach.Namely, they include classical electrostatic interactions that are damped to zero in the short-range.This approach goes back to the origins of small molecule potentials 30 and more modern versions of it have been incorporated into global methods that do not use range cutoffs, e.g., permutationally invariant polynomial (PIP)based potentials.Examples include the PIP potentials for CH 5 +31 and (H 2 O) 2 . 32While these do fit data in the longrange, it was also clear that when available, switching to accurate long-range interactions made sense.For example, for the 2-body interaction in a many-body expansion of the water potential q-AQUA, the long-range interaction is switched to the dipole−dipole interaction 33 using an accurate, flexible dipole moment surface. 34In order to avoid "ripples" in the final PES, it is important to verify that the machine-learned (ML) PES overlaps the long-range analytical potential accurately in the region of the switch, as was done in q-AQUA. 33n the absence of a high-quality long-range interaction for flexible monomers, we have utilized two component fits. 4One component is a precise fit to data in the short-range, and the second one is a fit to data in the long-range.In the example in ref 4 for SiO + H 2 , the long-range data, i.e., CCSD(T) energies, extend to 11.1 Å.The RMS fitting error for the longrange data is 0.05 cm −1 .This level of precision is just not feasible for a single fit to the entire data set.
In this paper, we present a full-dimensional potential energy surface for the NaCl−H 2 system, including the correct longrange behavior.The structure of the paper is outlined as follows: Section 2 provides a detailed description of the theoretical and computational aspects, including ab initio calculations, the fitting procedure, long-range interaction, and DMC calculations.The results and properties of the fitted potential energy surface are discussed in Section 3, followed by the summary and conclusions in Section 4.

COMPUTATIONAL DETAILS
The interaction potential for the NaCl−H 2 system in the electronic ground state has been computed using electronic energies generated predominantly on a six dimensional grid defined in Jacobi coordinates (r 1 , r 2 , R, θ 1 , θ 2 , ϕ) as shown in Figure 1.In this context, R is the distance between the centers of mass of NaCl and H 2 , while r 1 and r 2 indicate the respective bond lengths of the NaCl and H 2 molecule.Furthermore, θ 1 and θ 2 depict in-plane orientation angles between r 1 and R, and r 2 and R, while ϕ is the out-of-plane dihedral angle.In the ab initio calculations, R is varied in the range of 2.6−12.6Å, while the bond distances are within the intervals of 2.16 Å ≤ r 1 ≤ 2.66 and 0.50 Å ≤ r 2 ≤ 1.15 Å.The angular coordinates are within the ranges of 0°≤ θ 1 ≤ 360°and 0°≤ θ 2 (ϕ) ≤ 180°w here θ 1 = θ 2 = 0°corresponds to the collinear arrangement Na−Cl−H−H.
All electronic energy calculations were performed using the coupled-cluster method with single and double excitations and perturbative treatment of triple [CCSD(T)] using the MOLPRO suite of quantum chemistry software programs. 35,36he ab initio calculations employed an augmented correlationconsistent polarized valence triple-ζ basis set (aug-cc-pVTZ). 37 The NaCl−H 2 potential is given by the "plug and play" form, i.e., it is written as the sum of an interaction potential plus potentials for isolated NaCl and H 2 , as we have done for several potentials.4,38 NaCl and H 2 monomers were fit using the simple sum of powers of the Morse variables.Note that here we employed a straightforward 1-d fit for the NaCl and H 2 monomers at the CCSD(T)/aug-cc-pVTZ level (abbreviated below as aVTZ). Hover, users have the flexibility to provide any monomer potential.At each NaCl−H 2 configuration, the interaction energy is given by the electronic energy of the NaCl−H 2 minus the electronic energy of NaCl and the electronic energy of H 2 .Clearly, this interaction potential goes to zero asymptotically as NaCl and H 2 separate.For the previous analogous potential for SiO−H 2 , CCSD(T) energies were obtained for R as large as 11 Å such that the interaction potential is less than 1 cm −1 .4 Note that the long-range interaction for NaCl−H 2 is stronger than that for SiO + H 2 .To obtain a precise fit to these small asymptotic interaction energies, the PES fit in the short range (V SR ) was joined with the analytical long-range interaction (V LR ) with a smooth switching function as follows where Note that the choice of switching range R i and R f was motivated by the intention to capture a region where the direct CCSD(T) interaction energies exhibit a strong agreement with the analytical dipole−quadrupole interaction.
A large data set of ab initio points was calculated at the CCSD(T) level and used to fit the short-range PES.The PES has been fitted in 6D using an invariant polynomial method.(3) where c n n ... The total polynomial order for the fits is 7 resulting in a total of 918 linear coefficients.These were obtained using standard linear least-squares, utilizing the MSA software 2.0, 39,40 which extends the original code 41,42 to produce and also incorporate gradients data.The leading interaction potential between NaCl and H 2 in the long-range is given by the dipole−quadrupole interaction 26 V R where μ is the dipole moment of NaCl and Θ is the quadrupole moment of the H 2 molecule.Instead of using this directly, we use the equivalent expression employing Coulomb's law where q i are the point charges on NaCl and H 2 and d ij are the corresponding distances between the charges on the different molecules.The charges on (q Na and q Cl ) are determined using the dipole moment q r ( . ) obtained from ref 43.The charges on the H 2 molecule are determined by noting that the quadrupole moment can be written as −2q H at the center of the molecule as shown in Figure 2. Therefore, the charges of the H 2 molecule at a given bond length are determined using the quadrupole moment obtained from ref 44.Using the given equation q x q a q q a q a q r Figure 3 shows the variation of the NaCl dipole moment and the corresponding charge and the variation of the H 2 quadrupole moment and the effective charge with respect to the corresponding internuclear distances.We use eq 5 instead of eq 4 as the former is general and can be used for any term in the multipole expansion.We do note that other shorter range terms, such as the NaCl quadrupole and H 2 quadrupole interaction, are not considered since the NaCl quadrupole moment as a function of the NaCl distance is not known.−47 were performed to obtain the quantum zero point energy (ZPE) and wave function of the NaCl−H 2 system using the PES fit and to examine the quality of a PES in extended regions of the configuration space.The fundamental objective of DMC calculation is to solve the time-dependent Schrodinger wave equation in imaginary time τ.DMC calculations begin with an initial guess of the ground-state wave function represented by a population of N(0) equally weighted Gaussian random walkers.Subsequently, these walkers diffuse randomly subject to the potential in imaginary time following a Gaussian distribution.Each walker has the possibility to persist (potentially giving birth to a new walker) or be eliminated; the population is regulated by birth-death processes, as described below

DMC. Unbiased diffusion Monte Carlo (DMC) calculations
) 8)   where E i is the energy of the i-th walker, E r is the reference energy used to stabilize the diffusion system in its ground state, and Δτ is the step size in imaginary time.To keep the number of random walkers around the initial value N( 0), E r is modified at the end of each time step as per the following equation.
where N(τ) is the number of walkers at the imaginary time step τ, ⟨V(τ)⟩ is the average potential energy of all the alive walkers and α is the parameter that regulates the fluctuations in the number of walkers and the reference energy.The average of E r over a sufficiently long time period provides an estimate of the ZPE.Here, we performed five DMC calculations, each with an imaginary time step of Δτ = 5 au and α = 0.1.For every DMC calculation, we propagated a set of 20,000 random walkers from the global minimum for 30,000 time steps (≈3.62 ps), out of which 20,000 steps are used to equilibrate the walkers, and the reference energies in the remaining 10,000 steps are used to compute the ZPE.The Cartesian coordinates

The Journal of Physical Chemistry A
of the walkers at the last ten steps were recorded, and their interatomic distances were determined.

Precision of the Fits.
The potential energy surfaces are fitted using the ab initio data, with a cut-off based on the maximum interaction energy.The precision of the PES fits was determined by calculating the root-mean-square (RMS) fitting errors as a function of the maximum interaction energy.In this work, we choose two interaction energy cut-offs of 20,000 and 2000 cm −1 as shown in Table 1.The RMS fitting errors in the short-range fit are 1.42 cm −1 for PES I and 1.23 cm −1 for PES II.Note that our emphasis is on sampling the ab initio interaction energies essential for achieving an accurate fit.Although opting for an interaction energy range up to 2000 cm −1 is a suitable choice, computing additional points in the higher energy range, up to 20,000 cm −1 , is done for the completeness and to prevent the holes in the PES.Correlation plots are given for the two fits in Figure 4.As seen, there is nearly perfect correlation for both fits over the respective ranges of energies.
Table 2 presents the optimized geometrical structure and dissociation energy of the minimum for the PES fits with zero at the relaxed isolated NaCl and H 2 asymptote, comparing it with the MOLPRO calculated CCSD(T)/aVTZ optimized structures.Hereafter, PES refers to PES I.The global minimum is not for the T-shaped structure, and evidently this indicates that "chemical" effects beyond electrostatics are responsible for the global minimum.It might be interesting for those with more expertise than we have to analyze the source(s) responsible for this structure.
In Figure 5, several 1D cuts for the interaction energy between NaCl−H 2 are shown at different (θ 1 , θ 2 , ϕ) values from PES and are compared with the ab initio calculations using the aVTZ basis sets.The bond distances of NaCl and H 2 are fixed at their equilibrium values, r 1 = 2.413 Å and r 2 = 0.743 Å for Figure 5a−c, whereas for Figure 5d, bond lengths are fixed to the values obtained from the optimized minimum geometry listed in Table 2.As seen, over the large range shown, there is excellent agreement between direct CCSD(T) energies and the PES, including in the long-range, R greater than 11 Å where the PES is given by the analytical dipole− quadrupole interaction.Note in panel (b) there is a small sudden variation in the PES energy at 10 Å where the switching to the analytical interaction begins.This is actually due to a correction to the PES that slightly underestimates the ab initio energies between around 8−10 Å.The switch restores excellent agreement with the ab initio energies.
Figure 6 shows the potential energy cut as a function of ϕ.All other coordinates are fixed in the NaCl−H 2 equilibrium configuration.As seen, there is a large energy change of roughly 1400 cm −1 , suggesting that the complex will maintain a mostly planar configuration.

Long-Range Interaction.
A close examination of the dipole−quadrupole interaction is shown in Figure 7, along with direct CCSD(T) energies.As seen, the maximum deviation of the long-range interaction from the direct CCSD(T) energies is less than 1 cm −1 for R = 11 Å, which becomes almost exact with a further increase in R.
3.3.Zero-Point Wave Function.Next, we present results from the DMC calculations.When we conducted the unconstrained DMC, the PES was found to be "hole-free", meaning that the PES does not contain any configurations with unphysical negative energies.The average ZPE obtained from   The Journal of Physical Chemistry A five independent DMC simulations using the PES is 2778.16± 1.72 cm −1 .Figure 8 shows histograms of the wave function vs the NaCl and H 2 internuclear distances, r 1 and r 2 , respectively, and R. As seen, the results for r 1 and r 2 closely resemble harmonic oscillator wave functions for the ground state.For R, where the motion is strongly anharmonic, the wave function displays an asymmetry which has more amplitude at large Rvalues than at small values.The Journal of Physical Chemistry A

SUMMARY AND CONCLUSIONS
We reported a new potential energy surface for NaCl + H 2 from precise fitting of thousands of CCSD(T)/aug-cc-pVTZ energies.This PES is notable, as it is extended to infinite separation of NaCl and H 2 by means of an accurate dipole− quadrupole interaction for the flexible monomers.The PES is expressed as the sum of an interaction potential plus flexible monomer potentials, using the "plug and play" strategy adopted previously. 38he electronic energies can be accessed at ref 48 and the PES is available as Supporting Information.

The Journal of Physical Chemistry A
and R i = 10.0 Å, R f = 11.0Å, and b

1 6
is the linear coefficient and y i are the Morse variables of form exp(−d i /λ).The parameter λ is subject to user specification, which is set to 3.0a 0 for V SR .The internuclear distances d i are denoted as d 1 = d NaCl , d 2 = d NaH , d 3 = d NaH , d 4 = d ClH , d 5 = d ClH , and d 6 = d HH .

Figure 2 .
Figure 2. Electric quadrupole moment of the H 2 molecule.

Figure 3 .
Figure 3. Dipole moment and calculated effective charge as a function of bond length (a) dipole moment and calculated effective charge for NaCl, and (b) quadrupole moment and calculated effective charge for H 2 molecule, with circles denoting the values at the equilibrium bond lengths.The dipole moment and quadrupole moment data are taken from refs 43 and 44.

Figure 8 .
Figure 8. Histogram of cuts of the ground state DMC wave function vs the NaCl and H 2 internuclear distances, r 1 and r 2 , respectively, and R.

Table 1 .
RMS Fitting Error of Two PES Fits as a Function of the Maximum Energy Cut-Off a aEnergies are in cm −1 .

Table 2 .
Optimized Structures of NaCl−H 2 and Dissociation Energies of D e (cm −1 ) a a Distances are in Å and angles are in degrees.