Pair potential for helium from symmetry-adapted perturbation theory calculations and from supermolecular data

An efficient method by combining different basis sets and SAPT levels

In symmetry-adapted perturbation theory (SAPT), accurate calculations on non-covalent interaction (NCI) for large complexes with more than 50 atoms are time-consuming using large basis sets. More efficient ones with smaller basis sets usually result in poor prediction in terms of dispersion and overall energies. In this study, we propose two composite methods with baseline calculated at SAPT2/aug-cc-pVDZ and SAPT2/aug-cc-pVTZ with dispersion term corrected at SAPT2+ level using bond functions and smaller basis set with δ MP2 corrections respectively. Benchmark results on representative NCI data sets, such as S22, S66, and so forth, show significant improvements on the accuracy compared to the original SAPT Silver standard and comparable to SAPT Gold standard in some cases with much less computational cost.

Efficient Density-Fitted Explicitly Correlated Dispersion and Exchange Dispersion Energies

The leading-order dispersion and exchange-dispersion terms in symmetry-adapted perturbation theory (SAPT), E_disp⁽²⁰⁾ and E_exch-disp⁽²⁰⁾, suffer from slow convergence to the complete basis set limit. To alleviate this problem, explicitly correlated variants of these corrections, E_disp⁽²⁰⁾-F12 and E_exch-disp⁽²⁰⁾-F12, have been proposed recently. However, the original formalism (M., Kodrycka , J. Chem. Theory Comput. 2019, 15, 5965-5986), while highly successful in terms of improving convergence, was not competitive to conventional orbital-based SAPT in terms of computational efficiency due to the need to manipulate several kinds of two-electron integrals. In this work, we eliminate this need by decomposing all types of two-electron integrals using robust density fitting. We demonstrate that the error of the density fitting approximation is negligible when standard auxiliary bases such as aug-cc-pVXZ/MP2FIT are employed. The new implementation allowed us to study all complexes in the A24 database in basis sets up to aug-cc-pV5Z, and the E_disp⁽²⁰⁾-F12 and E_exch-disp⁽²⁰⁾-F12 values exhibit vastly improved basis set convergence over their conventional counterparts. The well-converged E_disp⁽²⁰⁾-F12 and E_exch-disp⁽²⁰⁾-F12 numbers can be substituted for conventional E_disp⁽²⁰⁾ and E_exch-disp⁽²⁰⁾ ones in a calculation of the total SAPT interaction energy at any level (SAPT0, SAPT2+3, ...). We show that the addition of F12 terms does not improve the accuracy of low-level SAPT treatments. However, when the theory errors are minimized in high-level SAPT approaches such as SAPT2+3(CCD)δMP2, the reduction of basis set incompleteness errors thanks to the F12 treatment substantially improves the accuracy of small-basis calculations.

SAPT codes for calculations of intermolecular interaction energies

Symmetry-adapted perturbation theory (SAPT) is a method for calculations of intermolecular (noncovalent) interaction energies. The set of SAPT codes that is described here, the current version named SAPT2020, includes virtually all variants of SAPT developed so far, among them two-body SAPT based on perturbative, coupled cluster, and density functional theory descriptions of monomers, three-body SAPT, and two-body SAPT for some classes of open-shell monomers. The properties of systems governed by noncovalent interactions can be predicted only if potential energy surfaces (force fields) are available. SAPT is the preferred approach for generating such surfaces since it is seamlessly connected to the asymptotic expansion of interaction energy. SAPT2020 includes codes for automatic development of such surfaces, enabling generation of complete dimer surfaces with a rigid monomer approximation for dimers containing about one hundred atoms. These codes can also be used to obtain surfaces including internal degrees of freedom of monomers.

Explicitly Correlated Dispersion and Exchange Dispersion Energies in Symmetry-Adapted Perturbation Theory

The individual interaction energy terms in symmetry-adapted perturbation theory (SAPT) not only have different physical interpretations but also converge to their complete basis set (CBS) limit values at quite different rates. Dispersion energy is notoriously the slowest converging interaction energy contribution, and exchange dispersion energy, while smaller in absolute value, converges just as poorly in relative terms. To speed up the basis set convergence of the lowest-order SAPT dispersion and exchange dispersion energies, we borrow the techniques from explicitly correlated (F12) electronic structure theory and develop practical expressions for the closed-shell E_disp⁽²⁰⁾-F12 and E_exch-disp⁽²⁰⁾-F12 contributions. While the latter term has been derived and implemented for the first time, the former correction was recently proposed by Przybytek [ J. Chem. Theory Comput. 2018 , 14 , 5105 - 5117 ] using an Ansatz with a full optimization of the explicitly correlated amplitudes. In addition to reimplementing the fully optimized variant of E_disp⁽²⁰⁾-F12, we propose three approximate Ansätze that substantially improve the scaling of the method and at the same time avoid the numerical instabilities of the unrestricted optimization. The performance of all four resulting flavors of E_disp⁽²⁰⁾-F12 and E_exch-disp⁽²⁰⁾-F12 is first tested on helium, neon, argon, water, and methane dimers, with orbital and auxiliary basis sets up to aug-cc-pV5Z and aug-cc-pV5Z-RI, respectively. The double- and triple-ζ basis set calculations are then extended to the entire A24 database of noncovalent interaction energies and compared with CBS estimates for E_disp⁽²⁰⁾ and E_exch-disp⁽²⁰⁾ computed using conventional SAPT with basis sets up to aug-cc-pV6Z with midbond functions. It is shown that the F12 treatment is highly successful in improving the basis set convergence of the SAPT terms, with the F12 calculations in an X-tuple ζ basis about as accurate as conventional calculations in bases with cardinal numbers (X + 2) for E_disp⁽²⁰⁾ and either (X + 1) or (X + 2) for E_exch-disp⁽²⁰⁾. While the full amplitude optimization affords the highest accuracy for both corrections, the much simpler and numerically stable optimized diagonal Ansatz is a very close second. We have also tested the performance of the simple F12 correction based on the second-order Møller-Plesset perturbation theory, SAPT-F12(MP2) [ Frey , J. A. ; Chem. Rev. 2016 , 116 , 5614 - 5641 ] and observed that it is also quite successful in speeding up the basis set convergence of conventional E_disp⁽²⁰⁾ + E_exch-disp⁽²⁰⁾, albeit with some outliers.

Nine questions on energy decomposition analysis

The paper collects the answers of the authors to the following questions: Is the lack of precision in the definition of many chemical concepts one of the reasons for the coexistence of many partition schemes? Does the adoption of a given partition scheme imply a set of more precise definitions of the underlying chemical concepts? How can one use the results of a partition scheme to improve the clarity of definitions of concepts? Are partition schemes subject to scientific Darwinism? If so, what is the influence of a community's sociological pressure in the "natural selection" process? To what extent does/can/should investigated systems influence the choice of a particular partition scheme? Do we need more focused chemical validation of Energy Decomposition Analysis (EDA) methodology and descriptors/terms in general? Is there any interest in developing common benchmarks and test sets for cross-validation of methods? Is it possible to contemplate a unified partition scheme (let us call it the "standard model" of partitioning), that is proper for all applications in chemistry, in the foreseeable future or even in principle? In the end, science is about experiments and the real world. Can one, therefore, use any experiment or experimental data be used to favor one partition scheme over another? © 2019 Wiley Periodicals, Inc.

Rotational Echoes as a Tool for Investigating Ultrafast Collisional Dynamics of Molecules

We show that recently discovered rotational echoes of molecules provide an efficient tool for studying collisional molecular dynamics in high-pressure gases. Our study demonstrates that rotational echoes enable the observation of extremely fast collisional dissipation, at timescales of the order of a few picoseconds, and possibly shorter. The decay of the rotational alignment echoes in CO_{2} gas and CO_{2}-He mixture up to 50 bar was studied experimentally, delivering collision rates that are in good agreement with the theoretical expectations. The suggested measurement protocol may be used in other high-density media, and potentially in liquids.

Cold atom-atom-anion three-body recombination of 4He4HexLi- (x = 6 or 7) systems

Atom-atom-anion three-body recombination (TBR) in mixed ⁴He and ^xLi^- (x = 6 or 7) is investigated in the adiabatic hyperspherical representation by quantum mechanically solving the Schrödinger equation. The distributions of product states following these TBR processes are found to be relatively different for the two systems when the collision energy is less than roughly 0.6 mK × k_B or 0.3 mK × k_B for ⁴He⁴He⁶Li^- and ⁴He⁴He⁷Li^- systems, respectively, with k_B being the Boltzmann constant. For ⁴He⁴He⁶Li^- systems, the rate of recombination into (v=0) l = 0⁴He⁶Li^- molecular anions is the largest with v and l denoting the rovibrational quantum numbers, while the TBR rate that leads to the formation of l = 1⁴He⁶Li^- molecular anions is a little smaller than that of neutral ⁴He₂ molecules. For ⁴He⁴He⁷Li^- systems, neutral ⁴He₂ molecules tend to be the most products, following the yields of l = 0 and 1 ⁴He⁷Li^- molecular anions. However, in spite of these distinctly different distributions, the products of molecular anions, the sum of l = 0 and 1 ⁴He^xLi^- products, are relatively larger than that of neutral ⁴He₂ molecules for both the two systems.

Fully quantum calculation of the second and third virial coefficients of water and its isotopologues from ab initio potentials

Path-Integral Monte Carlo methods were applied to calculate the second, B(T), and the third, C(T), virial coefficients for water. A fully quantum approach and state-of-the-art flexible-monomer pair and three-body potentials were used. Flexible-monomer potentials allow calculations for any isotopologue; we performed calculations for both H2O and D2O. For B(T) of H2O, the quantum effect contributes 25% of the value at 300 K and is not entirely negligible even at 1000 K, in accordance with recent literature findings. The effect of monomer flexibility, while not as large as some claims in the literature, is significant compared to the experimental uncertainty. It is of opposite sign to the quantum effect, smaller in magnitude than the latter below 500 K, and varies from 2% at 300 K to 10% at 700 K. When monomer flexibility is accounted for, results from the CCpol-8sf pair potential are in excellent agreement with the available experimental data and provide reliable B(T) values at temperatures outside the range of experimental data. The flexible-monomer MB-pol pair potential yields B(T) values that are slightly too high compared to experiment. For C(T), our calculations confirm earlier findings that the use of three-body potential is necessary for meaningful predictions. However, due to various uncertainties of the potentials used, especially the three-body ones, we were not able to establish benchmark values of C(T), although our results are in qualitative agreement with available experimental data. The quantum effect, never before included for water, reduces the magnitude of the classical value for H2O by a factor of 2.5 at 300 K and is not entirely negligible even at 1000 K.

Repulsion-dispersion parameters for the modelling of organic molecular crystals containing N, O, S and Cl

In lattice energy models that combine ab initio and empirical components, it is important to ensure consistency between these components so that meaningful quantitative results are obtained. A method for deriving parameters of atom-atom repulsion dispersion potentials for crystals, tailored to different ab initio models, is presented. It is based on minimization of the sum of squared deviations between experimental and calculated structures and energies. The solution algorithm is designed to avoid convergence to local minima in the parameter space by combining a deterministic low-discrepancy sequence for the generation of multiple initial parameter guesses with an efficient local minimization algorithm. The proposed approach is applied to derive transferable exp-6 potential parameters suitable for use in conjunction with a distributed multipole electrostatics model derived from isolated molecule charge densities calculated at the M06/6-31G(d,p) level of theory. Data for hydrocarbons, azahydrocarbons, oxohydrocarbons, organosulphur compounds and chlorohydrocarbons are used for the estimation. A good fit is achieved for the new set of parameters with a mean absolute error in sublimation enthalpies of 4.1 kJ mol-1 and an average rmsd15 of 0.31 Å. The parameters are found to perform well on a separate cross-validation set of 39 compounds.

Pair Potential with Submillikelvin Uncertainties and Nonadiabatic Treatment of the Halo State of the Helium Dimer

The pair potential for helium is computed with accuracy improved by an order of magnitude relative to the best previous determination. For the well region, its uncertainties are now below 1 millikelvin. The main improvement is due to the use of explicitly correlated wave functions at the nonrelativistic Born-Oppenheimer (BO) level of theory. The diagonal BO and the relativistic corrections are obtained from large full configuration interaction calculations. Nonadiabatic perturbation theory is used to predict the properties of the halo state of the helium dimer. Its binding energy and the average value of the interatomic distance are found to be 138.9(5) neV and 47.13(8) Å. The binding energy agrees with its first experimental determination of 151.9(13.3) neV [Zeller et al., Proc. Natl. Acad. Sci. U.S.A. 113, 14651 (2016)PNASA60027-842410.1073/pnas.1610688113].

Cold atom-atom-ion three-body recombination of 4He-4He-X- (X = H or D)

The atom-atom-ion three-body recombination (TBR) of mixed ⁴He and X^- (X = H or D) systems is investigated by solving the Schrödinger equation using the adiabatic hyperspherical representation method. It is shown that the dominant products after a TBR in the ultracold limit (E ≤ 0.1 mK × k_B) are different for the two systems. For the ⁴He⁴HeH^- system, the rate of TBR into the ⁴HeH^- ion is nearly two orders of magnitude larger than that of TBR into the neutral ⁴He₂ molecule. In contrast, the yield of ⁴He₂ is a little higher than that of ⁴HeD^- for the ⁴He⁴HeD^- system. Furthermore, since the adiabatic potentials become more attractive and the nonadiabatic couplings become much stronger by substituting D for H in the ⁴He⁴HeH^- system, the total TBR rate for the ⁴He⁴HeD^- system is nearly two orders of magnitude larger than that for the ⁴He⁴HeH^- system.

Quantum Halo States in Helium Tetramers

The universality of quantum halo states enables a comparison of systems from different fields of physics, as demonstrated in two- and three-body clusters. In the present work, we studied weakly bound helium tetramers in order to test whether some of these four-body realistic systems qualify as halos. Their ground-state binding energies and structural properties were thoroughly estimated using the diffusion Monte Carlo method with pure estimators. Helium tetramer properties proved to be less sensitive on the potential model than previously evaluated trimer properties. We predict the existence of realistic four-body halo ⁴He₂³He₂, whereas ⁴He₄ and ⁴He₃³He are close to the border and thus can be used as prototypes of quasi-halo systems. Our results could be tested by the experimental determination of the tetramers' structural properties using a setup similar to the one developed for the study of helium trimers.

Ground-state properties of weakly bound helium-alkali trimers

Weakly bound triatomic molecules consisting of two helium atoms and one alkali metal atom are studied by means of the diffusion Monte Carlo method. We determined the stability of ⁴He₂A, ⁴He³HeA, and ³He₂A, where A is one of the alkali atoms Li, Na, K, Rb, or Cs. Some of the trimers with ³He are predicted to be self-bound for the first time, but this is observed to be dependent on the He-A interaction potential model. In addition to the ground-state energy of the trimers, we determined their density, radial, and angular distributions. Many of them are spatially very extended, which qualifies them as quantum halo states.

Imaging the He2 quantum halo state using a free electron laser

Quantum tunneling is a ubiquitous phenomenon in nature and crucial for many technological applications. It allows quantum particles to reach regions in space which are energetically not accessible according to classical mechanics. In this "tunneling region," the particle density is known to decay exponentially. This behavior is universal across all energy scales from nuclear physics to chemistry and solid state systems. Although typically only a small fraction of a particle wavefunction extends into the tunneling region, we present here an extreme quantum system: a gigantic molecule consisting of two helium atoms, with an 80% probability that its two nuclei will be found in this classical forbidden region. This circumstance allows us to directly image the exponentially decaying density of a tunneling particle, which we achieved for over two orders of magnitude. Imaging a tunneling particle shows one of the few features of our world that is truly universal: the probability to find one of the constituents of bound matter far away is never zero but decreases exponentially. The results were obtained by Coulomb explosion imaging using a free electron laser and furthermore yielded He2's binding energy of [Formula: see text] neV, which is in agreement with most recent calculations.

Beyond Born-Mayer: Improved Models for Short-Range Repulsion in ab Initio Force Fields

Short-range repulsion within intermolecular force fields is conventionally described by either Lennard-Jones (A/r(12)) or Born-Mayer (A exp(-Br)) forms. Despite their widespread use, these simple functional forms are often unable to describe the interaction energy accurately over a broad range of intermolecular distances, thus creating challenges in the development of ab initio force fields and potentially leading to decreased accuracy and transferability. Herein, we derive a novel short-range functional form based on a simple Slater-like model of overlapping atomic densities and an iterated stockholder atom (ISA) partitioning of the molecular electron density. We demonstrate that this Slater-ISA methodology yields a more accurate, transferable, and robust description of the short-range interactions at minimal additional computational cost compared to standard Lennard-Jones or Born-Mayer approaches. Finally, we show how this methodology can be adapted to yield the standard Born-Mayer functional form while still retaining many of the advantages of the Slater-ISA approach.

Next-Generation Force Fields from Symmetry-Adapted Perturbation Theory

Symmetry-adapted perturbation theory (SAPT) provides a unique set of advantages for parameterizing next-generation force fields from first principles. SAPT provides a direct, basis-set superposition error free estimate of molecular interaction energies, a physically intuitive energy decomposition, and a seamless transition to an asymptotic picture of intermolecular interactions. These properties have been exploited throughout the literature to develop next-generation force fields for a variety of applications, including classical molecular dynamics simulations, crystal structure prediction, and quantum dynamics/spectroscopy. This review provides a brief overview of the formalism and theory of SAPT, along with a practical discussion of the various methodologies utilized to parameterize force fields from SAPT calculations. It also highlights a number of applications of SAPT-based force fields for chemical systems of particular interest. Finally, the review ends with a brief outlook on the future opportunities and challenges that remain for next-generation force fields based on SAPT.

Scalable improvement of SPME multipolar electrostatics in anisotropic polarizable molecular mechanics using a general short-range penetration correction up to quadrupoles

We propose a general coupling of the Smooth Particle Mesh Ewald SPME approach for distributed multipoles to a short-range charge penetration correction modifying the charge-charge, charge-dipole and charge-quadrupole energies. Such an approach significantly improves electrostatics when compared to ab initio values and has been calibrated on Symmetry-Adapted Perturbation Theory reference data. Various neutral molecular dimers have been tested and results on the complexes of mono- and divalent cations with a water ligand are also provided. Transferability of the correction is adressed in the context of the implementation of the AMOEBA and SIBFA polarizable force fields in the TINKER-HP software. As the choices of the multipolar distribution are discussed, conclusions are drawn for the future penetration-corrected polarizable force fields highlighting the mandatory need of non-spurious procedures for the obtention of well balanced and physically meaningful distributed moments. Finally, scalability and parallelism of the short-range corrected SPME approach are addressed, demonstrating that the damping function is computationally affordable and accurate for molecular dynamics simulations of complex bio- or bioinorganic systems in periodic boundary conditions.

The van der Waals interactions in rare-gas dimers: the role of interparticle interactions

We investigate the potential energy curves of rare-gas dimers with various ranges and strengths of interparticle interactions (nuclear–electron, electron–electron, and nuclear–nuclear interactions).

Universality in molecular halo clusters

The ground state of weakly bound dimers and trimers with a radius extending well into the classically forbidden region is explored, with the goal to test the predicted universality of quantum halo states. The focus of the study is molecules consisting of T↓, D↓, ^{3}He, ^{4}He, and alkali atoms, where the interaction between particles is much better known than in the case of nuclei, which are traditional examples of quantum halos. The study of realistic systems is supplemented by model calculations in order to analyze how low-energy properties depend on the interaction potential. The use of variational and diffusion Monte Carlo methods enabled a very precise calculation of both the size and binding energy of the trimers. In the quantum halo regime, and for large values of scaled binding energies, all clusters follow almost the same universal line. As the scaled binding energy decreases, Borromean states separate from tango trimers.

Determination of structure and properties of molecular crystals from first principles

CONSPECTUS: Until recently, it had been impossible to predict structures of molecular crystals just from the knowledge of the chemical formula for the constituent molecule(s). A solution of this problem has been achieved using intermolecular force fields computed from first principles. These fields were developed by calculating interaction energies of molecular dimers and trimers using an ab initio method called symmetry-adapted perturbation theory (SAPT) based on density-functional theory (DFT) description of monomers [SAPT(DFT)]. For clusters containing up to a dozen or so atoms, interaction energies computed using SAPT(DFT) are comparable in accuracy to the results of the best wave function-based methods, whereas the former approach can be applied to systems an order of magnitude larger than the latter. In fact, for monomers with a couple dozen atoms, SAPT(DFT) is about equally time-consuming as the supermolecular DFT approach. To develop a force field, SAPT(DFT) calculations are performed for a large number of dimer and possibly also trimer configurations (grid points in intermolecular coordinates), and the interaction energies are then fitted by analytic functions. The resulting force fields can be used to determine crystal structures and properties by applying them in molecular packing, lattice energy minimization, and molecular dynamics calculations. In this way, some of the first successful determinations of crystal structures were achieved from first principles, with crystal densities and lattice parameters agreeing with experimental values to within about 1%. Crystal properties obtained using similar procedures but empirical force fields fitted to crystal data have typical errors of several percent due to low sensitivity of empirical fits to interactions beyond those of the nearest neighbors. The first-principles approach has additional advantages over the empirical approach for notional crystals and cocrystals since empirical force fields can only be extrapolated to such cases. As an alternative to applying SAPT(DFT) in crystal structure calculations, one can use supermolecular DFT interaction energies combined with scaled dispersion energies computed from simple atom-atom functions, that is, use the so-called DFT+D approach. Whereas the standard DFT methods fail for intermolecular interactions, DFT+D performs reasonably well since the dispersion correction is used not only to provide the missing dispersion contribution but also to fix other deficiencies of DFT. The latter cancellation of errors is unphysical and can be avoided by applying the so-called dispersionless density functional, dlDF. In this case, the dispersion energies are added without any scaling. The dlDF+D method is also one of the best performing DFT+D methods. The SAPT(DFT)-based approach has been applied so far only to crystals with rigid monomers. It can be extended to partly flexible monomers, that is, to monomers with only a few internal coordinates allowed to vary. However, the costs will increase relative to rigid monomer cases since the number of grid points increases exponentially with the number of dimensions. One way around this problem is to construct force fields with approximate couplings between inter- and intramonomer degrees of freedom. Another way is to calculate interaction energies (and possibly forces) "on the fly", i.e., in each step of lattice energy minimization procedure. Such an approach would be prohibitively expensive if it replaced analytic force fields at all stages of the crystal predictions procedure, but it can be used to optimize a few dozen candidate structures determined by other methods.

The structure of the asymmetric helium trimer (3)He(4)He2

Despite their apparent simplicity, the properties of the two helium trimers, (4)He3 and (3)He(4)He2, are still not completely understood. In particular, the existence of a bound state of the asymmetric trimer (3)He(4)He2 was established many years ago, using different theoretical approaches, and later it was experimentally detected. However its structural properties have not been thoroughly investigated so far, probably because an accurate theoretical description of this very weakly bound system is computationally quite demanding. In this work we give for the first time an accurate and complete theoretical description of the geometrical structure of this fragile system using quantum Monte Carlo techniques employing the TTY potential and compare its properties with those of (4)He2 and (4)He3. We compute average values of distances and angles, along with the angle-angle distribution function: a two-dimensional probability distribution well suited to discuss the shape of a trimer. Our analysis shows that the lighter isotope is very diffuse and can be found at large distances from the other two atoms, but also close to the center of mass of the system in nearly linear configurations. For this system the concept of "equilibrium structure" is meaningless and all kinds of three-atom configurations must be taken into account in its description.

Symmetry-adapted perturbation theory study for some magnesium complexes

A systematic theoretical study on Mg–ligand interactions has been carried out employing both ab initio correlated wave function and density functional methods. The interactions of the Mg(CH3N2)2moiety with BF, CO, N2, NH3, and H2O ligands have been investigated by performing calculations at the B3LYP, MP2, MP4, and CCSD(T)/6–311++G(3df,3pd) levels of theory. Results indicate that the interaction energies of the Mg(CH3N2)2–L complexes increase in the order NH3 > H2O > BF > CO > N2. Symmetry-adapted perturbation theory (SAPT) analysis has been carried out to understand the nature of the forces involved in the bonding. The SAPT results indicate that the stabilities of the Mg–L interactions are attributed mainly to electrostatic effects, while induction and dispersion forces also play a significant role. The evaluated SAPT interaction energies for the Mg(CH3N2)2–L complexes are generally in good agreement with those obtained using the supermolecule CCSD(T) methods, suggesting that SAPT is a proper method to study the intermolecular interactions in these complexes. The results also suggest an explanation for the unique role of Mg2+as a carrier of water molecules that mediate enzymatic hydrolysis reactions.

Force field development for actinyl ions via quantum mechanical calculations: an approach to account for many body solvation effects

Advances in computational algorithms and methodologies make it possible to use highly accurate quantum mechanical calculations to develop force fields (pair-wise additive intermolecular potentials) for condensed phase simulations. Despite these advances, this approach faces numerous hurdles for the case of actinyl ions, AcO2(n+) (high-oxidation-state actinide dioxo cations), mainly due to the complex electronic structure resulting from an interplay of s, p, d, and f valence orbitals. Traditional methods use a pair of molecules (“dimer”) to generate a potential energy surface (PES) for force field parametrization based on the assumption that many body polarization effects are negligible. We show that this is a poor approximation for aqueous phase uranyl ions and present an alternative approach for the development of actinyl ion force fields that includes important many body solvation effects. Force fields are developed for the UO2(2+) ion with the SPC/Fw, TIP3P, TIP4P, and TIP5P water models and are validated by carrying out detailed molecular simulations on the uranyl aqua ion, one of the most characterized actinide systems. It is shown that the force fields faithfully reproduce available experimental structural data and hydration free energies. Failure to account for solvation effects when generating PES leads to overbinding between UO2(2+) and water, resulting in incorrect hydration free energies and coordination numbers. A detailed analysis of arrangement of water molecules in the first and second solvation shell of UO2(2+) is presented. The use of a simple functional form involving the sum of Lennard-Jones + Coulomb potentials makes the new force field compatible with a large number of available molecular simulation engines and common force fields.