An accurate and efficient computational protocol for obtaining the complete basis set limits of the binding energies of water clusters at the MP2 and CCSD(T) levels of theory: Application to (H2O)m, m = 2-6, 8, 11, 16, and 17

Systematic analysis of electronic barrier heights and widths for concerted proton transfer in cyclic hydrogen bonded clusters: (HF)n, (HCl)n and (H2O)n where n = 3, 4, 5

The MP2 and CCSD(T) methods are paired with correlation consistent basis sets as large as aug-cc-pVQZ to optimize the structures of the cyclic minima for (HF)n, (HCl)n and (H2O)n where n = 3-5, as well as the corresponding transition states (TSs) for concerted proton transfer (CPT). MP2 and CCSD(T) harmonic vibrational frequencies confirm the nature of each minimum and TS. Both conventional and explicitly correlated CCSD(T) computations are employed to assess the electronic dissociation energies and barrier heights for CPT near the complete basis (CBS) limit for all 9 clusters. Results for (HF)n are consistent with prior studies identifying Cnh and Dnh point group symmetry for the minima and TSs, respectively. Our computations also confirm that CPT proceeds through Cs TS structures for the C1 minima of (H2O)3 and (H2O)5, whereas the process goes through a TS with D2d symmetry for the S4 global minimum of (H2O)4. This work corroborates earlier findings that the minima for (HCl)3, (HCl)4 and (HCl)5 have C3h, S4 and C1 point group symmetry, respectively, and that the Cnh structures are not minima for n = 4 and 5. Moreover, our computations show the TSs for CPT in (HCl)3, (HCl)4 and (HCl)5 have D3h, D2d, and C2 point group symmetry, respectively. At the CCSD(T) CBS limit, (HF)4 and (HF)5 have the smallest electronic barrier heights for CPT (≈15 kcal mol-1 for both), followed by the HF trimer (≈21 kcal mol-1). The barriers are appreciably higher for the other clusters (around 27 kcal mol-1 for (H2O)4 and (HCl)3; roughly 30 kcal mol-1 for (H2O)3, (H2O)5 and (HCl)4; up to 38 kcal mol-1 for (HCl)5). At the CBS limit, MP2 significantly underestimates the CCSD(T) barrier heights (e.g., by ca. 2, 4 and 7 kcal mol-1 for the pentamers of HF, H2O and HCl, respectively), whereas CCSD overestimates these barriers by roughly the same magnitude. Scaling the barrier heights and dissociation energies by the number of fragments in the cluster reveals strong linear relationships between the two quantities and with the magnitudes of the imaginary vibrational frequency for the TSs.

Insight into the Binding of Argon to Cyclic Water Clusters from Symmetry-Adapted Perturbation Theory

This work systematically examines the interactions between a single argon atom and the edges and faces of cyclic H2O clusters containing three-five water molecules (Ar(H2O)n=3-5). Full geometry optimizations and subsequent harmonic vibrational frequency computations were performed using MP2 with a triple-ζ correlation consistent basis set augmented with diffuse functions on the heavy atoms (cc-pVTZ for H and aug-cc-pVTZ for O and Ar; denoted as haTZ). Optimized structures and harmonic vibrational frequencies were also obtained with the two-body-many-body (2b:Mb) and three-body-many-body (3b:Mb) techniques; here, high-level CCSD(T) computations capture up through the two-body or three-body contributions from the many-body expansion, respectively, while less demanding MP2 computations recover all higher-order contributions. Five unique stationary points have been identified in which Ar binds to the cyclic water trimer, along with four for (H2O)4 and three for (H2O)5. To the best of our knowledge, eleven of these twelve structures have been characterized here for the first time. Ar consistently binds more strongly to the faces than the edges of the cyclic (H2O)n clusters, by as much as a factor of two. The 3b:Mb electronic energies computed with the haTZ basis set indicate that Ar binds to the faces of the water clusters by at least 3 kJ mol-1 and by nearly 6 kJ mol-1 for one Ar(H2O)5 complex. An analysis of the interaction energies for the different binding motifs based on symmetry-adapted perturbation theory (SAPT) indicates that dispersion interactions are primarily responsible for the observed trends. The binding of a single Ar atom to a face of these cyclic water clusters can induce perturbations to the harmonic vibrational frequencies on the order of 5 cm-1 for some hydrogen-bonded OH stretching frequencies.

Water Charge Transfer Accelerates Criegee Intermediate Reaction with H2O- Radical Anion at the Aqueous Interface

Criegee intermediates (CIs) are important carbonyl oxides that may react with atmospheric trace chemicals and impact the global climate. The CI reaction with water has been widely studied and is a main channel for trapping CIs in the troposphere. Previous experimental and computational reports have largely focused on reaction kinetic processes in various CI-water reactions. The molecular-level origin of CI's interfacial reactivity at the water microdroplet surface (e.g., as found in aerosols and clouds) is unclear. In this study, by employing the quantum mechanical/molecular mechanical (QM/MM) Born-Oppenheimer molecular dynamics with the local second-order Møller-Plesset perturbation theory, our computational results reveal a substantial water charge transfer up to ∼20% per water, which creates the surface H₂O⁺/H₂O^- radical pairs to enhance the CH₂OO and anti-CH₃CHOO reactivity with water: the resulting strong CI-H₂O^- electrostatic attraction at the microdroplet surface facilitates the nucleophilic attack to the CI carbonyl by water, which may counteract the apolar hindrance of the substituent to accelerate the CI-water reaction. Our statistical analysis of the molecular dynamics trajectories further resolves a relatively long-lived bound CI(H₂O^-) intermediate state at the air/water interface, which has not been observed in gaseous CI reactions. This work provides insights into what may alter the oxidizing power of the troposphere by the next larger CIs than simple CH₂OO and implicates a new perspective on the role of interfacial water charge transfer in accelerating molecular reactions at aqueous interfaces.

The primary gas phase hydration shell of hydroxide

The number of water molecules in hydroxide's primary hydration shell has been long debated to be three from the interpretation of experimental data and four from theoretical studies. Here, we provide direct evidence for the presence of a fourth water molecule in hydroxide's primary hydration shell from a combined study based on high-resolution cryogenic experimental photoelectron spectroscopy and high-level quantum chemical computations. Well-defined spectra of OH^-(H₂O)_n clusters (n = 2 to 5) yield accurate electron binding energies, which are, in turn, used as key signatures of the underlying molecular conformations. Although the smaller OH^-(H₂O)₃ and OH^-(H₂O)₄ clusters adopt close-lying conformations with similar electron binding energies that are hard to distinguish, the OH^-(H₂O)₅ cluster clearly has a predominant conformation with a four-coordinated hydroxide binding motif, a finding that unambiguously determines the gas phase coordination number of hydroxide to be four.

An extensive assessment of the performance of pairwise and many-body interaction potentials in reproducing ab initio benchmark binding energies for water clusters n = 2-25

We assess the performance of 7 pairwise additive (TIP3P, TIP4P, TIP4P-ice, TIP5P, OPC, SPC, SPC/E) and 8 families of many-body potentials (q-AQUA, HIPPO, AMOEBA, EFP, TTM, WHBB, MB-pol, MB-UCB) in reproducing high-level ab initio benchmark values, CCSD(T) or MP2 at the complete basis set (CBS) limit for the binding energy and the many-body expansion (MBE) of water clusters n = 2-11, 16-17, 20, 25. By including a large range of cluster sizes having dissimilar hydrogen bonding networks, we obtain an understanding of how these potentials perform for different hydrogen bonding arrangements that are mostly outside of their parameterization range. While it is appropriate to compare the results of ab initio based many-body potentials directly to the electronic binding energies (D_e's), the pairwise additive ones are compared to the enthalpies at T = 298 K, ΔH(298 K), as the latter class of force fields are parametrized to reproduce enthalpies (implicitly accounting for zero-point energy corrections) rather than binding energies. We find that all pairwise additive potentials considered overestimate the reference ΔH values for the n = 2-25 clusters by >13%. For the water dimer (n = 2) in particular, the errors are in the range 83-119% for the pairwise additive potentials studied since these are based on an effective rather than the true 2-body interaction specifically designed as a means of partially accounting for the missing many-body terms. This stronger 2-body interaction is achieved by an enhanced monomer dipole moment that mimics its increase from the gas phase monomer to the condensed phase value. Indeed, for cluster sizes n ≥ 4 the percent deviations become slightly smaller (albeit all exceeding 13%). In contrast, we find that the many-body potentials perform more accurately in reproducing the electronic binding energies (D_e's) throughout the entire cluster range (n = 2-25), all reproducing the ab initio benchmark binding energies within ±7% of the respective CBS values. We further assess the ability of a subset of the many-body potentials (MB-UCB, q-AQUA, MB-pol, and TTM2.1-F) to also reproduce the magnitude of the ab initio many-body energy terms for water cluster sizes n = 7, 10, 16 and 17. The potentials show an overall good agreement with the available benchmark values. However, we identify characteristic differences upon comparing the many-body terms at both the ab initio-optimized geometries and the respective potential-optimized geometries to the reference ab initio values. Additionally, by applying this analysis to a wide range of cluster sizes, trends in the MBE of the potentials with increasing cluster size can be identified. Finally, in an attempt to draw a parallel between the pairwise additive and many-body potentials, we report the analysis of the individual molecular dipole moments for water clusters with 1 to ∼4 solvation shells with the TTM2.1-F potential. We find that the internally solvated water molecules have in general a larger molecular dipole moment ranging from 2.6-3.0 D. This justifies the use of an enhanced, with respect to the gas-phase value, molecular dipole moment for the pairwise additive potentials, which is intended to fold in the many body terms into an effective (enhanced) pairwise interaction through the choice of the charges. These results have important implications for the development of future generations of efficient, transferable, and highly accurate classical interaction potentials in both the pairwise additive and many-body categories.

A critical comparison of CH⋯π versus π⋯π interactions in the benzene dimer: obtaining benchmarks at the CCSD(T) level and assessing the accuracy of lower scaling methods

We have established CCSD(T)/CBS (Complete Basis Set) limits for 3 stationary points on the benzene dimer potential energy surface, corresponding to the π⋯π (parallel displaced or PD(C_2h), minimum) and CH⋯π (T-shaped or T(C_2v), transition state) and tilted T-shaped (or TT(C_s), minimum) bonding scenarios considering both the structure and binding energy. The CCSD(T)/CBS binding energies are -2.65 ± 0.02 (PD), -2.74 ± 0.03 (T), and -2.83 ± 0.01 kcal mol^-1 (TT). To this end, the CH⋯π is ∼0.2 kcal mol^-1 stronger than the π⋯π interaction, whereas the tilting of the CH donating benzene molecule with respect to the other benzene is worth 0.1 kcal mol^-1. As previously discussed in the literature, the MP2 level of theory does not provide a close match for either the energy or structure, yet the SCS-MP2 yields structures in excellent agreement with respect to the CCSD(T) result. It is found that the SCS-MI-MP2 also gives optimized structures very close to SCS-MP2 (within ∼0.01 Å of the benchmark). Despite the closer match in structure, the spin-biased MP2 methods (SCS-, SCS-MI-, and SOS-MP2) incorrectly predict the relative stabilities of the isomers. That said, none of the spin biased MP2 methods offers a good compromise between energy and structure for the systems examined. Finally, the CCSD(T)/CBS benchmarks were used to assess the performance of 13 DFT functionals selected from different rungs of Jacob's ladder. Several functionals such as TPSS-D3, B3LYP-D3, B97-D, B97-D3, and B2PLYP-D3 provided a good description of the binding energies for both CH⋯π and π⋯π interactions, yielding values within 6% of the CCSD(T)/CBS benchmark values. Unlike the MP2 methods, these functionals correctly predict the relative stability of the PD(C_2h) and T(C_2v) dimers. Further, we find that there is no systematic improvement as Jacob's ladder is ascended (increased complexity of functional). The best functionals that result in a good compromise between structure and energy accuracy are B97-D3 and B2PLYP-D3 for both the CH⋯π and π⋯π interaction. Despite the impressive performance of these functionals, a challenge that remains is ensuring the transferability of these density functionals in accurately describing the interaction between dimers of larger aromatic molecules, the latter requiring high-level benchmarks for these systems.

Adsorption free energy of phenol onto coronene: Solvent and temperature effects

Molecular modeling can considerably speed up the discovery of materials with high adsorption capacity for wastewater treatment. Despite considerable efforts in computational studies, the molecular modeling of adsorption processes has several limitations in reproducing experimental conditions. Handling the environmental effects (solvent effects) and the temperature effects are part of the important limitations in the literature. In this work, we address these two limitations using the adsorption of phenol onto coronene as case study. In the proposed model, for the solvent effects, we used a hybrid solvation model, with n explicit water molecules and implicit solvation. We increasingly used n=1 to n=12 explicit water molecules. To account for the temperature effects, we evaluated the adsorption efficiency using the adsorption free energy for temperatures varying from 200 to 400K. We generated initial configurations using classical molecular dynamics, before further optimisation at the ωB97XD/aug-cc-pVDZ level of theory. Polarisable continuum solvation model (PCM) is used for the implicit solvation. The adsorption free energy is evaluated to be -1.3kcal/mol at room temperature. It has been found that the adsorption free energy is more negative at low temperatures. Above 360K, the adsorption free energy is found to be positive.

Spontaneous Formation of Hydrogen Peroxide in Water Microdroplets

There is accumulating evidence that many chemical reactions are accelerated by several orders of magnitude in micrometer-sized aqueous or organic liquid droplets compared to their corresponding bulk liquid phase. However, the molecular origin of the enhanced rates remains unclear as in the case of spontaneous appearance of 1 μM hydrogen peroxide in water microdroplets. In this Letter, we consider the range of ionization energies and whether interfacial electric fields of a microdroplet can feasibly overcome the high energy step from hydroxide ions (OH^-) to hydroxyl radicals (OH^•) in a primary H₂O₂ mechanism. We find that the vertical ionization energies (VIEs) of partially solvated OH^- ions are greatly lowered relative to the average VIE in the bulk liquid, unlike the case of the Cl^- anion which shows no reduction in the VIEs regardless of solvation environment. Overall reduced hydrogen-bonding and undercoordination of OH^- are structural features that are more readily present at the air-water interface, where the energy scale for ionization can be matched by statistically probable electric field values.

Experimental and theoretical characterization of chelidonic acid structure

Chelidonic acid (4-oxo-4H-pyran-2,6-dicarboxylic acid) is present in plants of Papaveraceae family, especially in Chelidonium majus. Due to its anticancer, antibacterial, hepatoprotective, and antioxidant properties, it has been used in medical treatments. In this work, the X-ray structure of methanol solvate of chelidonic acid was determined. Layers of chelidonic acid are held by hydrogen bonds via COOH and C = O fragments and additionally bridged by methanol. The formed H-bond network between two acid units is different from typical –COOH dimers observed, e.g., in crystals of isophtalic acid. The molecular structure of 2,6-dimethyl-γ-pyrone (2Me4PN) and chelidonic acid, a 2,6-dicarboxylic derivate of γ-pyrone (4PN), was verified in silico using density functional theory (DFT-B3LYP) combined with large correlation-consistent basis sets. The impact of –CH3 and –COOH substituents on 4PN ring structure, dipole moments, geometric/magnetic indexes of aromaticity, and NBO charges was assessed following unconstrained geometry optimization in the gas phase, chloroform, methanol, DMSO, and water with solvent effect introduced using the polarized continuous model (PCM). H-bond network formed in chelidonic acid–methanol complex was analyzed and their interaction energy estimated. Theoretical modeling enabled prediction of accurate structural parameters, dipole moments, and geometric/magnetic indexes of aromaticity of the studied 4PN, 2Me4PN, and chelidonic acid molecules.

Modeling of laser-pulse induced small water cluster-(H2O)N (N = 1-10) decomposition on suitable metal cluster catalysts

Understanding the microscopic mechanisms of electronic excitation in water clusters is a very important and challenging problem in a series of solar energy applications, such as solar water evaporation, photolysis, etc. Here we employ real time-time-dependent density functional theory (RT-TDDFT) and Ehrenfest dynamics to investigate the photodissociation dynamic process of (H₂O)_N=1-10 clusters and photoinduced charge transfer in them. The research presented here confirms that the plane tetramer, (H₂O)₄, is the most difficult one to be dissociated under laser irradiation in the ten clusters for its high (S₄) symmetry; the overall order of the ease of decomposition is as follows: (H₂O)_6-p > (H₂O)₈ > (H₂O)_6-c > (H₂O)₇ > (H₂O)₁₀ > (H₂O)₁ > (H₂O)₃ > (H₂O)₂ > (H₂O)₉ > (H₂O)₅ > (H₂O)₄. Plasmon catalyst-induced water splitting is a promising and feasible way to efficiently convert solar to chemical energy via reducing the laser amplitude threshold significantly; and among the Ag₆, Au₆, Cu₆, Al₆ chains and several Cu₆ clusters with O_h symmetry, the Cu₆ chain seems to be the most cost-effective one. This article aims at unraveling the fundamental mechanisms and providing valuable physical insights into the behavior of water splitting to pave the way for the theoretical and experimental design of the photolysis process.

Hydration Shell Model for Expeditious and Reliable Individual Hydrogen Bond Energies in Large Water Clusters

Recently, we have developed and tested a method, based on the molecular tailoring approach (MTA-based) to directly estimate the individual hydrogen bond (HB) energies in molecular clusters. Application of this...

Polarizable Water Potential Derived from a Model Electron Density

A new empirical potential for efficient, large scale molecular dynamics simulation of water is presented. The HIPPO (Hydrogen-like Intermolecular Polarizable POtential) force field is based upon the model electron density of a hydrogen-like atom. This framework is used to derive and parametrize individual terms describing charge penetration damped permanent electrostatics, damped polarization, charge transfer, anisotropic Pauli repulsion, and damped dispersion interactions. Initial parameter values were fit to Symmetry Adapted Perturbation Theory (SAPT) energy components for ten water dimer configurations, as well as the radial and angular dependence of the canonical dimer. The SAPT-based parameters were then systematically refined to extend the treatment to water bulk phases. The final HIPPO water model provides a balanced representation of a wide variety of properties of gas phase clusters, liquid water, and ice polymorphs, across a range of temperatures and pressures. This water potential yields a rationalization of water structure, dynamics, and thermodynamics explicitly correlated with an ab initio energy decomposition, while providing a level of accuracy comparable or superior to previous polarizable atomic multipole force fields. The HIPPO water model serves as a cornerstone around which similarly detailed physics-based models can be developed for additional molecular species.

Energies, structures, and harmonic frequencies of small water clusters from the direct random phase approximation

The binding energies, structures, and vibrational frequencies of water clusters up to 20 molecules are computed at the direct random phase approximation (RPA) level of theory and compared to theoretical benchmarks. Binding energies of the WATER27 set, which includes neutral and positively and negatively charged clusters, are predicted to be too low in the complete basis set limit by an average of 7 kcal/mol (9%) and are worse than the results from the best density functional theory methods or from the Møller-Plesset theory. The RPA shows significant basis set size dependence for binding energies. The order of the relative energies of the water hexamer and dodecamer isomers is predicted correctly by the RPA. The mean absolute deviation for angles and distances for neutral clusters up to the water hexamer are 0.2° and 0.6 pm, respectively, using quintuple-ζ basis sets. The relative energetic order of the hexamer isomers is preserved upon optimization. Vibrational frequencies for these systems are underestimated by several tens of wavenumbers for large basis sets, and deviations increase with the basis set size. Overall, the direct RPA method yields accurate structural parameters but systematically underestimates binding energies and shows strong basis set size dependence.

A unified cost-effective method for the construction of reliable potential energy surfaces for H2S and H2O clusters

A DFT-based methodology has been used to construct the potential energy surface of H2S clusters up to pentamers. Geometrical parameters and energetics show very good agreement with the existing experimental and high-level theoretical results. Distinct stable conformers of three dimers, six trimers, eleven tetramers and twenty-three pentamers have been identified. Both S-HS H-bond and SS interactions are identified in dimers, trimers and pentamers, while no SS interactions could be found in any of the 11 tetramer conformers. The binding energies of the most stable dimer, trimer, tetramer and pentamer are -1.66, -5.21, -8.57 and -12.54 kcal mol-1, respectively. The PES has been found to be exceedingly flat and the energy gap between the most and the least stable conformers was found to be only 0.09, 2.13, 1.65 and 1.13 kcal mol-1, from the dimer to the pentamer, respectively. The proposed method has also been used for water clusters up to the pentamer. The results obtained were found to agree closely with the existing results. Only one conformer was found for the water dimer, whereas four, five and fifteen conformers were obtained for the trimer, tetramer and pentamer, respectively. Atoms in molecular calculations were found to corroborate with the geometric and energetic results for both clusters.

Guest-Host Interactions in Clathrate Hydrates: Benchmark MP2 and CCSD(T)/CBS Binding Energies of CH4, CO2, and H2S in (H2O)20 Cages

We present benchmark binding energies of naturally occurring gas molecules CH₄, CO₂, and H₂S in the small cage, namely, the pentagonal dodecahedron (5¹²) (H₂O)₂₀, which is one of the constituent cages of the 3 major lattices (structures I, II, and H) of clathrate hydrates. These weak interactions require higher levels of electron correlation and converge slowly with an increasing basis set to the complete basis set (CBS) limit, necessitating the use of large basis sets up to the aug-cc-pV5Z and subsequent correction for basis set superposition error (BSSE). For the host hollow (H₂O)₂₀ cages, we have identified a most stable isomer with binding energy of -200.8 ± 2.1 kcal/mol at the CCSD(T)/CBS limit (-199.2 ± 0.5 kcal/mol at the MP2/CBS limit). Additionally, we report converged second order Møller-Plesset (MP2) CBS binding energies for the encapsulation of guests in the (H₂O)₂₀ cage of -4.3 ± 0.1 for CH₄@(H₂O)₂₀, -6.6 ± 0.1 for CO₂@(H₂O)₂₀, and -8.5 ± 0.1 kcal/mol for H₂S@(H₂O)₂₀, respectively. For CH₄@(H₂O)₂₀, exhibiting the weakest encapsulation affinity among the three, we report CCSD(T)/aug-cc-pVTZ binding energies and, based on them, a CCSD(T)/CBS estimate of -4.75 ± 0.1 kcal/mol. To the best of our knowledge, the CCSD(T)/aug-cc-pVTZ calculation for CH₄@(H₂O)₂₀ is the largest one reported to date (168 valence electrons, 1978 basis functions, and the correlation of 84 doubly occupied and 1873 virtual orbitals) and required a scalable implementation of the (T) module on 6144 nodes (350 208 cores) of the "Cori" supercomputer at the National Energy Research Supercomputing Center (NERSC) for a total execution time of 195 min (for the (T) part). These efficient scalable implementations of highly correlated methods offer the capability to obtain long-lasting benchmarks of intermolecular interactions in complex systems. They also provide a path toward parametrizing classical potentials needed to study the dynamical and transport properties in these complex systems as well as assess the accuracy of lower scaling electronic structure methods such as density functional theory (DFT) and MP2 including its spin-biased variants.

Quantum Chemistry Calculations for Metabolomics

A primary goal of metabolomics studies is to fully characterize the small-molecule composition of complex biological and environmental samples. However, despite advances in analytical technologies over the past two decades, the majority of small molecules in complex samples are not readily identifiable due to the immense structural and chemical diversity present within the metabolome. Current gold-standard identification methods rely on reference libraries built using authentic chemical materials ("standards"), which are not available for most molecules. Computational quantum chemistry methods, which can be used to calculate chemical properties that are then measured by analytical platforms, offer an alternative route for building reference libraries, i.e., in silico libraries for "standards-free" identification. In this review, we cover the major roadblocks currently facing metabolomics and discuss applications where quantum chemistry calculations offer a solution. Several successful examples for nuclear magnetic resonance spectroscopy, ion mobility spectrometry, infrared spectroscopy, and mass spectrometry methods are reviewed. Finally, we consider current best practices, sources of error, and provide an outlook for quantum chemistry calculations in metabolomics studies. We expect this review will inspire researchers in the field of small-molecule identification to accelerate adoption of in silico methods for generation of reference libraries and to add quantum chemistry calculations as another tool at their disposal to characterize complex samples.

Determining Partial Atomic Charges for Liquid Water: Assessing Electronic Structure and Charge Models

Partial atomic charges provide an intuitive and efficient way to describe the charge distribution and the resulting intermolecular electrostatic interactions in liquid water. Many charge models exist and it is unclear which model provides the best assignment of partial atomic charges in response to the local molecular environment. In this work, we systematically scrutinize various electronic structure methods and charge models (Mulliken, natural population analysis, CHelpG, RESP, Hirshfeld, Iterative Hirshfeld, and Bader) by evaluating their performance in predicting the dipole moments of isolated water, water clusters, and liquid water as well as charge transfer in the water dimer and liquid water. Although none of the seven charge models is capable of fully capturing the dipole moment increase from isolated water (1.85 D) to liquid water (about 2.9 D), the Iterative Hirshfeld method performs best for liquid water, reproducing its experimental average molecular dipole moment, yielding a reasonable amount of intermolecular charge transfer, and showing modest sensitivity to the local water environment. The performance of the charge model is dependent on the choice of the density functional and the quantum treatment of the environment. The computed molecular dipole moment of water generally increases with the percentage of the exact Hartree-Fock exchange in the functional, whereas the amount of charge transfer between molecules decreases. For liquid water, including two full solvation shells of surrounding water molecules (within about 5.5 Å of the central water) in the quantum chemical calculation converges the charges of the central water molecule. Our final pragmatic quantum chemical charge-assigning protocol for liquid water is the Iterative Hirshfeld method with M06-HF/aug-cc-pVDZ and a quantum region cutoff radius of 5.5 Å.

Structural Stability of the CO2 @sI Hydrate: a Bottom-Up Quantum Chemistry Approach on the Guest-Cage and Inter-Cage Interactions

Through reliable first-principles computations, we have demonstrated the impact of CO₂ molecules enclathration on the stability of sI clathrate hydrates. Given the delicate balance between the interaction energy components (van der Waals, hydrogen bonds) present on such systems, we follow a systematic bottom-up approach starting from the individual 5¹² and 5¹² 6² sI cages, up to all existing combinations of two-adjacent sI crystal cages to evaluate how such clathrate-like models perform on the evaluation of the guest-host and first-neighbors inter-cage effects, respectively. Interaction and binding energies of the CO₂ occupation of the sI cages were computed using DF-MP2 and different DFT/DFT-D electronic structure methodologies. The performance of selected DFT functionals, together with various semi-classical dispersion corrections schemes, were validated by comparison with reference ab initio DF-MP2 data, as well as experimental data from x-ray and neutron diffraction studies available. Our investigation confirms that the inclusion of the CO₂ in the cage/s is an energetically favorable process, with the CO₂ molecule preferring to occupy the large 5¹² 6² sI cages compared to the 5¹² ones. Further, the present results conclude on the rigidity of the water cages arrangements, showing the importance of the inter-cage couplings in the cluster models under study. In particular, the guest-cage interaction is the key factor for the preferential orientation of the captured CO₂ molecules in the sI cages, while the inter-cage interactions seems to cause minor distortions with the CO₂ guest neighbors interactions do not extending beyond the large 5¹² 6² sI cages. Such findings on these clathrate-like model systems are in accord with experimental observations, drawing a direct relevance to the structural stability of CO₂ @sI clathrates.

Temperature and Nuclear Quantum Effects on the Stretching Modes of the Water Hexamer

The water hexamer has many low-lying isomers, e.g., ring, book, cage, and prism, shifting from two- to three-dimensional structures. We show that this dimensionality change is accompanied by a drop in the quantum nature of the cluster, as manifested in the red shift of the quantal OH stretching modes as compared with their classical counterparts. We obtain this "nuclear quantum effect" (NQE) as the mean deviation between the OH stretch frequencies from velocity autocorrelation Fourier transforms from classical trajectories on a high-level water potential (MB-pol) as compared with scaled harmonic frequencies from high-level quantum chemistry calculations. With a universal scaling factor, the predicted OH frequencies agree with experiment to a mean absolute deviation ≤10 cm-1, which allows unequivocal isomer assignments. By assuming temperature-independent NQEs, we produce the temperature dependence of the cage isomer OH stretch spectrum below 70 K, where it is the dominant structure. All bands widen and blue-shift with increasing temperature, most conspicuously the reddest mode, which thus constitutes a "vibrational thermometer".

He Inclusion in Ice-like and Clathrate-like Frameworks: A Benchmark Quantum Chemistry Study of Guest-Host Interactions

Energetics and structural properties of selected type and size He@hydrate frameworks, e.g., from regular structured ice channels to clathrate-like cages, are presented from first-principles quantum chemistry methods. The scarcity of information on He@hydrates makes such complexes challenging targets, while their computational study entails an interesting and arduous task. Some of them have been synthesized in the laboratory, which motivates further investigations on their stability. Hence, the main focus is to examine the performance and accuracy of different wave function-based electronic structure methods, such as MP2, CCSD(T), their explicitly correlated (F12) and domain-based local pair-natural orbital (DLPNO) analogs, as well as modern and conventional density functional theory (DFT) approaches, and analytical model potentials available. Different structures are considered, starting from the "simplest system" formed by a noble gas atom (such as He) and one water molecule, followed by the study of the "fundamental units" present in all ice-like and clathrate-like frameworks (such as pentamers and hexamers) and finally the description of interactions in the "building blocks" of three-dimensional (3D) ice channels (e.g., horizontal and perpendicular ice II and I_h) and clathrate-like cages, such as the 5¹² present in the most common sI, sII, and sH clathrate-hydrate structures. The idea is to provide well-converged DLPNO-CCSD(T) and DFMP2/CBS reference datasets that in turn are used to validate how DFT functionals (in total, 29 approaches from generalized-gradient approximation (GGA), meta-GGA, to hybrid and range-separated functionals, including dispersion correction treatments, were checked) and analytical semiempirical/ab initio-based potentials perform compared with high-level alternatives. Within all tested approaches, those best-performing were identified and classified. Most of the DFT/DFT-D functionals, as well as available analytical pairwise model potentials, face difficulties in describing both hydrogen-bonded water frameworks and dispersion bound He-water interactions. Including dispersion corrections yields an overall well-balanced performance for LCωPBE-D3BJ and PBE0-D4 functionals. Such benchmark datasets can benefit research into the development of new cheminformatics models, as can serve to guide and cross-check methodologies, lending increased predicted power to future molecular simulations for investigating the role of structures and phase transitions from nanoscale clusters to macroscopic crystalline structures.

Infrared spectroscopy of neutral water clusters at finite temperature: Evidence for a noncyclic pentamer

Infrared spectroscopic study of neutral water clusters is crucial to understanding of the hydrogen-bonding networks in liquid water and ice. Here we report infrared spectra of size-selected neutral water clusters, (H₂O) _n (n = 3-6), in the OH stretching vibration region, based on threshold photoionization using a tunable vacuum ultraviolet free-electron laser. Distinct OH stretch vibrational fundamentals observed in the 3,500-3,600-cm^-1 region of (H₂O)₅ provide unique spectral signatures for the formation of a noncyclic pentamer, which coexists with the global-minimum cyclic structure previously identified in the gas phase. The main features of infrared spectra of the pentamer and hexamer, (H₂O) _n (n = 5 and 6), span the entire OH stretching band of liquid water, suggesting that they start to exhibit the richness and diversity of hydrogen-bonding networks in bulk water.

Thermally Induced Hydrogen-Bond Rearrangements in Small Water Clusters and the Persistent Water Tetramer

Small water clusters absorb heat and catalyze pivotal atmospheric reactions. Yet, experiments produced conflicting results on water cluster distribution under atmospheric conditions. Additionally, it is unclear which "phase transitions" such clusters exhibit, at what temperatures, and what are their underlying molecular mechanisms. We find that logarithmically small tails in the radial probability densities of (H₂O) _n clusters (n = 2 - 6) provide direct testimony for such transitions. Using the best available water potential (MB-pol), an advanced thermostating algorithm (g-BAOAB), and sufficiently long trajectories, we map the "bifurcation", "melting", and (hitherto unexplored) "vaporization" transitions, finding that both melting and vaporization proceed via a "monomer on a ring" conformer, exhibiting huge distance fluctuations at the vaporization temperatures (T _v). T _v may play a role in determining the atmospheric cluster size distribution such that the dimer and tetramer, with their exceptionally low/high T _v values, are under/over-represented in these distributions, as indeed observed in nondestructive mass spectrometric measurements.

Atlas of putative minima and low-lying energy networks of water clusters n = 3-25

We report a database consisting of the putative minima and ∼3.2 × 10⁶ local minima lying within 5 kcal/mol from the putative minima for water clusters of sizes n = 3-25 using an improved version of the Monte Carlo temperature basin paving (MCTBP) global optimization procedure in conjunction with the ab initio based, flexible, polarizable Thole-Type Model (TTM2.1-F, version 2.1) interaction potential for water. Several of the low-lying structures, as well as low-lying penta-coordinated water networks obtained with the TTM2.1-F potential, were further refined at the Møller-Plesset second order perturbation (MP2)/aug-cc-pVTZ level of theory. In total, we have identified 3 138 303 networks corresponding to local minima of the clusters n = 3-25, whose Cartesian coordinates and relative energies can be obtained from the webpage https://sites.uw.edu/wdbase/. Networks containing penta-coordinated water molecules start to appear at n = 11 and, quite surprisingly, are energetically close (within 1-3 kcal/mol) to the putative minima, a fact that has been confirmed from the MP2 calculations. This large database of water cluster minima spanning quite dissimilar hydrogen bonding networks is expected to influence the development and assessment of the accuracy of interaction potentials for water as well as lower scaling electronic structure methods (such as different density functionals). Furthermore, it can also be used in conjunction with data science approaches (including but not limited to neural networks and machine and deep learning) to understand the properties of water, nature's most important substance.

Quantum Monte Carlo Study of the Water Dimer Binding Energy and Halogen-π Interactions

Halogen-π systems are involved with competition between halogen bonding and π-interaction. Using the diffusion quantum Monte Carlo (DMC) method, we compare the interaction of benzene with diatomic halogens (X₂: Cl₂/Br₂) with the typical hydrogen bonding in the water dimer, taking into account explicit correlations of up to three bodies. The benzene-Cl₂/Br₂ binding energies (13.07 ± 0.42/16.62 ± 0.02 kJ/mol) attributed to both halogen bonding and dispersion are smaller than but comparable to the typical hydrogen bonding in the water dimer binding energy (20.88 ± 0.27 kJ/mol). All of the above values are in good agreement with those from the coupled-cluster with single, double, and noniterative triple excitations (CCSD(T)) results at the complete basis set limit (benzene-Cl₂/Br₂: 12.78/16.17 kJ/mol; water dimer: 21.0 kJ/mol).

Performance of polarization-consistent vs. correlation-consistent basis sets for CCSD(T) prediction of water dimer interaction energy

Detailed study of Jensen’s polarization-consistent vs. Dunning’s correlation-consistent basis set families performance on the extrapolation of raw and counterpoise-corrected interaction energies of water dimer using coupled cluster with single, double, and perturbative correction for connected triple excitations (CCSD(T)) in the complete basis set (CBS) limit are reported. Both 3-parameter exponential and 2-parameter inverse-power fits vs. the cardinal number of basis set, as well as the number of basis functions were analyzed and compared with one of the most extensive CCSD(T) results reported recently. The obtained results for both Jensen- and Dunning-type basis sets underestimate raw interaction energy by less than 0.136 kcal/mol with respect to the reference value of − 4.98065 kcal/mol. The use of counterpoise correction further improves (closer to the reference value) interaction energy. Asymptotic convergence of 3-parameter fitted interaction energy with respect to both cardinal number of basis set and the number of basis functions are closer to the reference value at the CBS limit than other fitting approaches considered here. Separate fits of Hartree-Fock and correlation interaction energy with 3-parameter formula additionally improved the results, and the smallest CBS deviation from the reference value is about 0.001 kcal/mol (underestimated) for CCSD(T)/aug-cc-pVXZ calculations. However, Jensen’s basis set underestimates such value to 0.012 kcal/mol. No improvement was observed for using the number of basis functions instead of cardinal number for fitting.

Surface Penetration without Enrichment: Simulations Show Ion Surface Propensities Consistent with Both Elevated Surface Tension and Surface Sensitive Spectroscopy

Molecular dynamics was used to investigate ion surface propensities in NaCl, KBr, and CsI solutions with an MP2-based force field. Although NaCl is found to be strongly repelled from the liquid-vapor interface, softer ions, such as I^-, penetrate closely to the interface. Despite the surface penetration, the concentration of CsI near the interface is still lower than that in the bulk, thus leading to no surface enrichment. The salt concentration is found to affect relative surface propensities of the ions. More significant surface penetration is observed at higher salt concentrations. Softer ions at higher concentrations form a complex multilayer arrangement that can not be characterized as a simple surface bilayer. The simulated ion distributions explain the spectroscopic evidence of surface perturbation by soft ions with a negative surface excess consistent with an increased surface tension of salt solutions.

Data-Driven Acceleration of the Coupled-Cluster Singles and Doubles Iterative Solver

Solving the coupled-cluster (CC) equations is a cost-prohibitive process that exhibits poor scaling with system size. These equations are solved by determining the set of amplitudes (t) that minimize the system energy with respect to the coupled-cluster equations at the selected level of truncation. Here, a novel approach to predict the converged coupled-cluster singles and doubles (CCSD) amplitudes, thus the coupled-cluster wave function, is explored by using machine learning and electronic structure properties inherent to the MP2 level. Features are collected from quantum chemical data, such as orbital energies, one-electron Hamiltonian, Coulomb, and exchange terms. The data-driven CCSD (DDCCSD) is not an alchemical method because the actual iterative coupled-cluster equations are solved. However, accurate energetics can also be obtained by bypassing solving the CC equations entirely. Our preliminary data show that it is possible to achieve remarkable speedups in solving the CCSD equations, especially when the correct physics are encoded and used for training of machine learning models.

Ring-Stacking Water Clusters: Morphology and Stabilities

The structures and interaction energies of water clusters with ring stacking motifs are studied by using ab initio calculations. The structures of the water clusters are constructed by stacking either single rings or multi-rings of tetramer, pentamer, and hexamer. We found that, in the single-ring-stacking motif, the most stable isomers exhibit an alternative clockwise-anticlockwise stacking pattern. We also show that four-layer single-ring-stacking isomers are not energetically favorable in comparison with those of two-layer multi-ring-stacking isomers. The relative stability of the isomers is also analyzed in terms of H-bond strength and elastic distortions of the water molecules.

Structures, relative stability and binding energies of neutral water clusters, (H2O)2–30

We have revised the structures of neutral water clusters, (H₂O)_n=2–30, with the affordable M06-2X functional, presenting up to 25 isomers for each cluster size.

Harmonizing accuracy and efficiency: A pragmatic approach to fragmentation of large molecules

Fragmentation methods offer an attractive alternative for ab initio treatment of large molecules and molecular clusters. However, balancing the accuracy and efficiency of these methods is a tight-rope-act. With this in view, we present an algorithm for automatic molecular fragmentation within Molecular Tailoring Approach (MTA) achieving this delicate balance. The automated code is tested out on a variety of molecules and clusters at the Hartree-Fock (HF)- and Møller-Plesset second order perturbation theory as well as density functional theory employing augmented Dunning basis sets. The results show remarkable accuracy and efficiency vis-à-vis the respective full calculations. Thus the present work forms an important step toward the development of an MTA-based black box code for implementation of HF as well as correlated quantum chemical calculations on large molecular systems.

Assessing Intermolecular Interactions in Guest-Free Clathrate Hydrate Systems

Recently, empty hydrate structures sI, sII, sH, and others have been proposed as low-density ice structures by both experimental observations and computer simulations. Some of them have been synthesized in the laboratory, which motivates further investigations on the stability of such guest-free clathrate structures. Using semiempirical and ab initio-based water models, as well as dispersion-corrected density functional theory approaches, we predict their stability, including cooperative many-body effects, in comparison with reference data from converged wave function-based DF-MP2 electronic structure calculations. We show that large basis sets and counterpoise corrections are required to improve convergence in the interaction/binding energies for such systems. Therefore, extrapolation schemes based on triple/quadruple and quadruple/quintuple ζ quality basis sets are used to reach high accuracy. Eleven different water structures corresponding to dodecahedron, edge sharing, face sharing, and fused cubes, as a part of the WATER27 database, as well as cavities from the sI, sII, and sH clathrate hydrates formed by 20, 24, 28, and 36 water molecules, are employed, and new benchmark energies are reported. Using these benchmark sets of interaction energies, we assess the performance of both analytical models and direct DFT calculations for such clathrate-like systems. In particular, seven popular water models (TIP4P/ice, TIP4P/2005, q-TIP4P/F, TTM2-F, TTM3-F, TTM4-F, and MB-pol) available in the literature, and nine density functional approximations (3 meta-GGAs, 3 hybrids, and 3 range separated functionals) are used to investigate their accuracy. By including dispersion corrections, our results show that errors in the interaction energies are reduced for most of the DFT functionals. Despite the difficulties faced by current water models and DFT functionals to accurately describe the interactions in such water systems, we found some general trends that could serve to extend their applicability to larger systems.

Understanding the structure and hydrogen bonding network of (H2O)32 and (H2O)33: an improved Monte Carlo temperature basin paving (MCTBP) method and quantum theory of atoms in molecules (QTAIM) analysis

Large water clusters are of particular interest because of their connection to liquid water and the intricate hydrogen bonding networks they possess.

Breaking the bottleneck: Use of molecular tailoring approach for the estimation of binding energies at MP2/CBS limit for large water clusters

A pragmatic method based on the molecular tailoring approach (MTA) for estimating the complete basis set (CBS) limit at Møller-Plesset second order perturbation (MP2) theory accurately for large molecular clusters with limited computational resources is developed. It is applied to water clusters, (H2O)n (n = 7, 8, 10, 16, 17, and 25) optimized employing aug-cc-pVDZ (aVDZ) basis-set. Binding energies (BEs) of these clusters are estimated at the MP2/aug-cc-pVNZ (aVNZ) [N = T, Q, and 5 (whenever possible)] levels of theory employing grafted MTA (GMTA) methodology and are found to lie within 0.2 kcal/mol of the corresponding full calculation MP2 BE, wherever available. The results are extrapolated to CBS limit using a three point formula. The GMTA-MP2 calculations are feasible on off-the-shelf hardware and show around 50%-65% saving of computational time. The methodology has a potential for application to molecular clusters containing ∼100 atoms.

Toward an Accurate and Inexpensive Estimation of CCSD(T)/CBS Binding Energies of Large Water Clusters

Owing to the steep scaling behavior, highly accurate CCSD(T) calculations, the contemporary gold standard of quantum chemistry, are prohibitively difficult for moderate- and large-sized water clusters even with the high-end hardware. The molecular tailoring approach (MTA), a fragmentation-based technique is found to be useful for enabling such high-level ab initio calculations. The present work reports the CCSD(T) level binding energies of many low-lying isomers of large (H2O)n (n = 16, 17, and 25) clusters employing aug-cc-pVDZ and aug-cc-pVTZ basis sets within the MTA framework. Accurate estimation of the CCSD(T) level binding energies [within 0.3 kcal/mol of the respective full calculation (FC) results] is achieved after effecting the grafting procedure, a protocol for minimizing the errors in the MTA-derived energies arising due to the approximate nature of MTA. The CCSD(T) level grafting procedure presented here hinges upon the well-known fact that the MP2 method, which scales as O(N(5)), can be a suitable starting point for approximating to the highly accurate CCSD(T) [that scale as O(N(7))] energies. On account of the requirement of only an MP2-level FC on the entire cluster, the current methodology ultimately leads to a cost-effective solution for the CCSD(T) level accurate binding energies of large-sized water clusters even at the complete basis set limit utilizing off-the-shelf hardware.

Challenging Dogmas: Hydrogen Bond Revisited

Hydrogen bond directionality in the water dimer is explained on the basis of symmetry-adapted intermolecular perturbation theory which directly separates the intermolecular interaction energy into four physically interpretable components: electrostatics, exchange-repulsion, dispersion, and induction. Analysis of these four main contributions to the binding energy allows a deeper understanding of the dominant factors ruling the mutual arrangement of the two monomers. A preference for the linear configuration is shown to be due to a subtle interplay of all four energy components. While the first-order terms, electrostatic and exchange-repulsion, almost perfectly cancel each other near the equilibrium geometry of the dimer, the importance of the second- and higher-order terms, induction and dispersion, becomes evident.

How Anion Chaotrope Changes the Local Structure of Water: Insights from Photoelectron Spectroscopy and Theoretical Modeling of SCN(-) Water Clusters

The behavior of charged solute molecules in aqueous solutions is often classified using the concept of kosmotropes ("structure makers") and chaotropes ("structure breakers"). There is a growing consensus that the key to kosmotropic/chaotropic behaviors lies in the local solvent region, but the exact microscopic basis for such differentiation is not well-understood. This issue is examined in this work by analyzing size selective solvation of a well-known chaotrope, a negatively charged SCN(-) molecule. Combining experimental photoelectron spectroscopy measurements with theoretical modeling, we examine evolution of solvation structure up to eight waters. We observe that SCN(-) indeed fits the description of weakly hydrated ion, and its solvation is heavily driven by stabilization of water-water interaction network. However, the impact on water structure is more subtle than that associated with "structure breaker". In particular, we observe that the solvation structure of SCN(-) preserves the "packing" structure of the water network but changes local directionality of hydrogen bonds in the local solvent region. The resulting effect is closer to that of "structure weakener", where solute can be readily accommodated into the native water network, at the cost of compromising its stability due to constraints on hydrogen bonding directionality.