On the Contribution of Vibrational Anharmonicity to the Binding Energies of Water Clusters

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.

Amino Acids Compete with Ammonia in Sulfuric Acid-Based Atmospheric Aerosol Prenucleation: The Case of Glycine and Serine

We present a computational investigation of the sulfuric acid, glycine, serine, ammonia, and water system to understand if this system can form prenucleation clusters, which are precursors to larger aerosols in the atmosphere. We have performed a comprehensive configurational search of all possible clusters in this system, starting with the four different monomers and zero to five waters. Accurate Gibbs free energies of formation have been calculated with the DLPNO-CCSD(T)/complete basis set (CBS) method on ωb97xd/6-31++G** geometries. For the dry dimers of sulfuric acid, the weakest base, serine, is found to form the most stable complex, which is a consequence of the strong di-ionic complex formed between the bisulfate ion and the protonated serine cation. For the dry dimers without sulfuric acid, the glycine-serine complex is more stable than the glycine-ammonia or serine-ammonia complexes, stemming from the detailed structure and not related to base strength. For the larger complexes, sulfuric acid deprotonates and the proton is shifted to glycine, serine, or ammonia. The two amino acids and ammonia are almost interchangeable and there is no easy way to predict which molecule will be protonated without the calculated results. Assuming reasonable starting concentrations and a closed system of sulfuric acid, glycine, serine, ammonia, and five waters, we predict the concentrations of all possible complexes at two temperatures spanning the troposphere. The most negative ΔG° values are a function of the detailed molecular interactions of these clusters. These details are more important than the base strength of ammonia, glycine, and serine.

The water hexamer: three-body interactions, structures, energetics, and OH-stretch spectroscopy at finite temperature

Using a newly developed and recently parameterized classical empirical simulation model for water that involves explicit three-body interactions, we determine the eleven most stable isomers of the water hexamer. We find that the lowest energy isomer is one of the cage structures, in agreement with far-IR and microwave experiments. The energy ordering for the binding energies is cage > glove > book > bag > chair > boat > chaise, and energies relative to the cage are in good agreement with CCSD(T) calculations. The three-body contributions to the cage, book, and chair are also in reasonable agreement with CCSD(T) results. The energy of each isomer results from a delicate balance involving the number of hydrogen bonds, the strain of these hydrogen bonds, and cooperative and anti-cooperative three-body interactions, whose contribution we can understand simply from the form of the three-body interactions in the simulation model. Oxygen-oxygen distances in the cage and book isomers are in good agreement with microwave experiments. Hydrogen-bond distances depend on both donor and acceptor, which can again be understood from the three-body model. Fully anharmonic OH-stretch spectra are calculated for these low-energy structures, and compared with shifted harmonic results from ab initio and density functional theory calculations. Replica-exchange molecular dynamics simulations were performed from 40 to 194 K, which show that the cage isomer has the lowest free energy from 0 to 70 K, and the book isomer has the lowest free energy from 70 to 194 K. OH-stretch spectra were calculated between 40 and 194 K, and results at 40, 63, and 79 K were compared to recent experiments, leading to re-assignment of the peaks in the experimental spectra. We calculate local OH-stretch cumulative spectral densities for different donor-acceptor types and compare to analogous results for liquid water.

Ionization of doped helium nanodroplets: complexes of C60 with water clusters

Water clusters are known to undergo an autoprotonation reaction upon ionization by photons or electron impact, resulting in the formation of (H(2)O)(n)H(3)O(+). Ejection of OH cannot be quenched by near-threshold ionization; it is only partly quenched when clusters are complexed with inert gas atoms. Mass spectra recorded by electron ionization of water-doped helium droplets show that the helium matrix also fails to quench OH loss. The situation changes drastically when helium droplets are codoped with C(60). Charged C(60)-water complexes are predominantly unprotonated; C(60)(H(2)O)(4)(+) and (C(60))(2)(H(2)O)(4)(+) appear with enhanced abundance. Another intense ion series is due to C(60)(H(2)O)(n)OH(+); dehydrogenation is proposed to be initiated by charge transfer between the primary He(+) ion and C(60). The resulting electronically excited C(60)(+*) leads to the formation of a doubly charged C(60)-water complex either via emission of an Auger electron from C(60)(+*), or internal Penning ionization of the attached water complex, followed by charge separation within {C(60)(H(2)O)(n)}(2+). This mechanism would also explain previous observations of dehydrogenation reactions in doped helium droplets. Mass-analyzed ion kinetic energy scans reveal spontaneous (unimolecular) dissociation of C(60)(H(2)O)(n)(+). In addition to the loss of single water molecules, a prominent reaction channel yields bare C(60)(+) for sizes n=3, 4, or 6. Ab initio Hartree-Fock calculations for C(60)-water complexes reveal negligible charge transfer within neutral complexes. Cationic complexes are well described as water clusters weakly bound to C(60)(+). For n=3, 4, or 6, fissionlike desorption of the entire water complex from C(60)(H(2)O)(n)(+) energetically competes with the evaporation of a single water molecule.

Anomalously strong effect of the ion sign on the thermochemistry of hydrogen bonded aqueous clusters of identical chemical composition

The sign preference of hydrogen bonded aqueous ionic clusters X±(H₂O)_i (n =1–5, X = F; Cl; Br) has been investigated using the Density Functional Theory and ab initio MP2 method. The present study indicates the anomalously large difference in formation free energies between cations and anions of identical chemical composition. The effect of vibrational anharmonicity on stepwise Gibbs free energy changes has been investigated, and possible uncertainties associated with the harmonic treatment of vibrational spectra have been discussed.

Association patterns in (HF)(m)(H2O)(n) (m + n = 2-8) clusters

In an attempt to understand the phase behavior of aqueous hydrogen fluoride, the clustering in the mixture is investigated at the molecular level. The study is performed at the mPW1B95/6-31+G(d,p) level of theory. Several previous studies attempted to describe the dissociation of HF in water, but in this investigation, the focus is only on the association patterns that are present in this binary mixture. A total of 214 optimized geometries of (HF)n(H2O)m clusters, with m + n as high as 8, were investigated. For each cluster combination, several different conformations are investigated, and the preferred conformations are presented. Using multiple linear regressions, the average strengths of the four possible H-bonding interactions are obtained. The strongest H-bond interaction is reported to be the H2O...H-F interaction. The most probable distributions of mixed clusters as a function of composition are also deduced. It is found that the larger (HF)n(H2O)m clusters are favored both energetically and entropically compared to the ones that are of size m + n < or = 3. Also, the clusters with equimolar contributions of HF and H2O are found to have the strongest interactions.

Are recent water models obtained by fitting diffraction data consistent with infrared/Raman and x-ray absorption spectra?

X-ray absorption (XA) spectra have been computed based on water structures obtained from a recent fit to x-ray and neutron diffraction data using models ranging from symmetrical to asymmetrical local coordination of the water molecules [A. K. Soper, J. Phys.: Condens. Matter 17, S3273 (2005)]. It is found that both the obtained symmetric and asymmetric structural models of water give similar looking XA spectra, which do not match the experiment. The fitted models both contain unphysical structures that are allowed by the diffraction data, where, e.g., hydrogen-hydrogen interactions may occur. A modification to the asymmetric model, in which the non-hydrogen-bonded OH intramolecular distance is allowed to become shorter while the bonded OH distance becomes longer, improves the situation somewhat, but the overall agreement is still unsatisfactory. The electric field (E-field) distributions and infrared (IR) spectra are also calculated using two established theoretical approaches, which, however, show significant discrepancies in their predictions for the asymmetric structural models. Both approaches predict the Raman spectrum of the symmetric model fitted to the diffraction data to be significantly blueshifted compared to experiment. At the moment no water model exists that can equally well describe IR/Raman, x-ray absorption spectroscopy, and diffraction data.

Theoretical calculations: Can Gibbs free energy for intermolecular complexes be predicted efficiently and accurately?

The theoretical study has been performed to refine the procedure for calculations of Gibbs free energy with a relative accuracy of less than 1 kcal/mol. Three benchmark intermolecular complexes are examined via several quantum-chemical methods, including the second-order Moller-Plesset perturbation (MP2), coupled cluster (CCSD(T)), and density functional (BLYP, B3LYP) theories augmented by Dunnings correlation-consistent basis sets. The effects of electron correlation, basis set size, and anharmonicity are systematically analyzed, and the results are compared with available experimental data. The results of the calculations suggest that experimental accuracy can be reached only by extrapolation of MP2 and CCSD(T) total energies to the complete basis set. The contribution of anharmonicity to the zero point energy and TDeltaSint values is fairly small. The new, economic way to reach chemical accuracy in the calculations of the thermodynamic parameters of intermolecular interactions is proposed. In addition, interaction energy (De) and free energy change (DeltaA) for considered species have been evaluated by Carr-Parrinello molecular dynamics (CPMD) simulations and static BLYP-plane wave calculations. The free energy change along the reaction paths were determined by the thermodynamic integration/"Blue Moon Ensemble" technique. Comparison between obtained values, and available experimental and conventional ab initio results has been made. We found that the accuracy of CPMD simulations is affected by several factors, including statistical uncertainty and convergence of constrained forces (TD integration), and the nature of DFT (density functional theory) functional. The results show that CPMD technique is capable of reproducing interaction and free energy with an accuracy of 1 kcal/mol and 2-3 kcal/mol respectively.

Analysis of Nuclear Quantum Effects on Hydrogen Bonding

The impact of nuclear quantum effects on hydrogen bonding is investigated for a series of hydrogen fluoride (HF)n clusters and a partially solvated fluoride anion, F-(H2O). The nuclear quantum effects are included using the path integral formalism in conjunction with the Car-Parrinello molecular dynamics (PICPMD) method and using the second-order vibrational perturbation theory (VPT2) approach. For the HF clusters, a directional change in the impact of nuclear quantum effects on the hydrogen-bonding strength is observed as the clusters evolve toward the condensed phase. Specifically, the inclusion of nuclear quantum effects increases the F-F distances for the (HF)n=2-4 clusters and decreases the F-F distances for the (HF)n>4 clusters. This directional change occurs because the enhanced electrostatic interactions between the HF monomers become more dominant than the zero point energy effects of librational modes as the size of the HF clusters increases. For the F-(H2O) system, the inclusion of nuclear quantum effects decreases the F-O distance and strengthens the hydrogen bonding interaction between the fluoride anion and the water molecule because of enhanced electrostatic interactions. The vibrationally averaged 19F shielding constant for F-(H2O) is significantly lower than the value for the equilibrium geometry, indicating that the electronic density on the fluorine decreases as a result of the quantum delocalization of the shared hydrogen. Deuteration of this system leads to an increase in the vibrationally averaged F-O distance and nuclear magnetic shielding constant because of the smaller degree of quantum delocalization for deuterium.

Theoretical Study on the Structure of the BrO Hydrates

The hydrates of bromine monoxide, BrO(H2O)n, n = 1-4, have been studied by means of ab initio calculations at the B3LYP/aug-cc-pVTZ level of theory. These systems could be formed in the troposphere and participate in chemical reactions involved in the depletion of ozone. Several conformations are obtained and discussed for each of the hydrates mentioned. Two rather different intermolecular interactions are found, namely, conventional hydrogen bonding and Br...O associations. In contrast with a more traditional point of view in which hydrogen bonds could be assumed as the preferential interaction for the formation of these complexes, it is the Br...O association which yields the most stable conformations. Equilibrium geometries, harmonic frequencies, and relative energies have been calculated for the bromine monoxide hydrates for the first time. The theoretical binding energies indicate that the stabilization of the hydrates increases with the number of water molecules added. Cooperative effects are suggested to play a significant role in this stabilization. An analysis of relevant properties depending on the electron density in the bond critical points of the Br...O associations has been done for the first time, showing characteristic features of this interaction in comparison with the hydrogen bonds formed.

Prediction of accurate anharmonic experimental vibrational frequencies for water clusters, (H2O)n, n=2-5

Accurate anharmonic experimental vibrational frequencies for water clusters consisting of 2-5 water molecules have been predicted on the basis of comparing different methods with MP2/aug-cc-pVTZ calculated and experimental anharmonic frequencies. The combination of using HF/6-31G* scaled frequencies for intramolecular modes and anharmonic frequencies for intermolecular modes gives excellent agreement with experiment for the water dimer and trimer and are as good as the expensive anharmonic MP2 calculations. The water trimer, the cyclic Ci and S4 tetramers, and the cyclic pentamer all have unique peaks in the infrared spectrum between 500 and 800 cm-1 and between 3400 and 3700 cm-1. Under the right experimental conditions these different clusters can be uniquely identified using high-resolution IR spectroscopy.

Probing the mechanism of hypoxia selectivity of copper bis(thiosemicarbazonato) complexes: DFT calculation of redox potentials and absolute acidities in solution

Density functional theory (DFT) calculations have been performed using the uB3LYP/6-31++G(d,p) model to calculate the solution phase one-electron reduction potentials (E(calc)) and absolute pKa values of a series of copper bis(thiosemicarbazonato) complexes. The effects of solvation in water and dimethylsulfoxide (DMSO) are incorporated as a self-consistent reaction field (SCRF) using the integral equation formalism polarisable continuum model (IEFPCM) and are found to be essential for quantitative agreement with an average error in E(calc) of -0.02 V compared to experiment. The bonding and spin densities are examined through the use of Natural Bond Order analysis and the results used to rationalise the calculated and observed reduction potentials. Calculated estimates of pKa values of several copper(II) species are presented and their implications for the mechanisms of transport and trapping within hypoxic cells are considered. Reduction is found to be a prerequisite for protonation of the complexes which suggests their transport in the blood stream as neutral species, and the mechanistic sequence is identified as a sequential electrochemical-chemical (EC) process. The complex equilibria of protonation, reoxidation and dissociation are discussed and the copper(I) diprotonated, cationic complex of diacetyl bis(4-methyl-3-thiosemicarbazonato)copper(II), Cu(I)ATSMH2(+), is identified as a possible candidate for the initial species trapped in hypoxic cells.

Density Functional Studies of Actinyl Aquo Complexes Studied Using Small-Core Effective Core Potentials and a Scalar Four-Component Relativistic Method

The title compounds, [AnO(2)(H(2)O)(5)](n)(+), n = 1 or 2 and An = U, Np, and Pu, are studied using relativistic density functional theory (DFT). Three rather different relativistic methods are used, small-core effective core potentials (SC-ECP), a scalar four-component all-electron relativistic method, and the zeroeth-order regular approximation. The methods provide similar results for a variety of properties, giving confidence in their accuracy. Spin-orbit and multiplet corrections to the An(VI)/An(V) reduction potential are added in an approximate fashion but are found to be essential. Bulk solvation effects are modeled with continuum solvation models (CPCM, COSMO). These models are tested by comparing explicit (cluster), continuum, and mixed cluster/continuum solvation models as applied to various properties. The continuum solvation models are shown to accurately account for the effects of the solvent, provided that at least the first coordination sphere is included. Reoptimizing the structures in the presence of the bulk solvent is seen to be important for the equatorial bond lengths but less relevant for energetics. Explicit inclusion of waters in the second coordination sphere has a modest influence on the energetics. For the first time, free energies of solvation are calculated for all six [AnO(2)(H(2)O)(5)](n)(+) species. The calculated numbers are within the experimental error margins, and the experimental trend is reproduced correctly. By comparison of different relativistic methods, it is shown that an accurate relativistic description leads to marked improvements over the older large-core ECP (LC-ECP) method for bond lengths, vibrational frequencies, and, in particular, the An(VI)/An(V) reduction potential. Two approximate DFT methods are compared, B3LYP, a hybrid DFT method, and PBE, a generalized gradient approximation. Either method yields An(VI)/An(V) reduction potentials of comparable quality. Overall, the experimental reduction potentials are accurately reproduced by the calculations.

Impact of Nuclear Quantum Effects on the Molecular Structure of Bihalides and the Hydrogen Fluoride Dimer

The structural impact of nuclear quantum effects is investigated for a set of bihalides, [XHX](-), X = F, Cl, and Br, and the hydrogen fluoride dimer. Structures are calculated with the vibrational self-consistent-field (VSCF) method, the second-order vibrational perturbation theory method (VPT2), and the nuclear-electronic orbital (NEO) approach. In the VSCF and VPT2 methods, the vibrationally averaged geometries are calculated for the Born-Oppenheimer electronic potential energy surface. In the NEO approach, the hydrogen nuclei are treated quantum mechanically on the same level as the electrons, and mixed nuclear-electronic wave functions are calculated variationally with molecular orbital methods. Electron-electron and electron-proton dynamical correlation effects are included in the NEO approach using second-order perturbation theory (NEO-MP2). The nuclear quantum effects are found to alter the distances between the heavy atoms by 0.02-0.05 A for the systems studied. These effects are of similar magnitude as the electron correlation effects. For the bihalides, inclusion of the nuclear quantum effects with the NEO-MP2 or the VSCF method increases the X-X distance. The bihalide X-X distances are similar for both methods and are consistent with two-dimensional grid calculations and experimental values, thereby validating the use of the computationally efficient NEO-MP2 method for these types of systems. For the hydrogen fluoride dimer, inclusion of nuclear quantum effects decreases the F-F distance with the NEO-MP2 method and increases the F-F distance with the VSCF and VPT2 methods. The VPT2 F-F distances for the hydrogen fluoride dimer and the deuterated form are consistent with the experimentally determined values. The NEO-MP2 F-F distance is in excellent agreement with the distance obtained experimentally for a model that removes the large amplitude bending motions. The analysis of these calculations provides insight into the significance of electron-electron and electron-proton correlation, anharmonicity of the vibrational modes, and nonadiabatic effects for hydrogen-bonded systems.

Water Clusters on Graphite: Methodology for Quantum Chemical A Priori Prediction of Reaction Rate Constants

The properties, interactions, and reactions of cyclic water clusters (H(2)O)(n=1-5) on model systems for a graphite surface have been studied using pure B3LYP, dispersion-augmented density functional tight binding (DFTB-D), and integrated ONIOM(B3LYP:DFTB-D) methods. Coronene C(24)H(12) as well as polycircumcoronenes C(96)H(24) and C(216)H(36) in monolayer, bilayer, and trilayer arrangements were used as model systems to simulate ABA bulk graphite. Structures, binding energies, and vibrational frequencies of water clusters on mono- and bilayer graphite models have been calculated, and structural changes and frequency shifts due to the water cluster-graphite interactions are discussed. ONIOM(B3LYP:DFTB-D) with coronene and water in the high level and C(96)H(24) in the low level mimics the effect of extended graphite pi-conjugation on the water-graphite interaction very reasonably and suggests that water clusters only weakly interact with graphite surfaces, as suggested by the fact that water is an excellent graphite lubricant. We use the ONIOM(B3LYP:DFTB-D) method to predict rate constants for model pathways of water dissociative adsorption on graphite. Quantum chemical molecular dynamics (QM/MD) simulations of water clusters and water addition products on the C(96)H(24) graphite model are presented using the DFTB-D method. A three-stage strategy is devised for a priori investigations of high temperature corrosion processes of graphite surfaces due to interaction with water molecules and fragments.

Global Search for Minimum Energy (H2O)n Clusters, n = 3−5

The Gaussian-3 (G3) model chemistry method has been used to calculate the relative deltaG(o) values for all possible conformers of neutral clusters of water, (H2O)n, where n = 3-5. A complete 12-fold conformational search around each hydrogen bond produced 144, 1728, and 20,736 initial starting structures of the water trimer, tetramer, and pentamer. These structures were optimized with PM3, followed by HF/6-31G* optimization, and then with the G3 model chemistry. Only two trimers are present on the G3 potential energy hypersurface. We identified 5 tetramers and 10 pentamers on the potential energy and free-energy hypersurfaces at 298 K. None of these 17 structures were linear; all linear starting models folded into cyclic or three-dimensional structures. The cyclic pentamer is the most stable isomer at 298 K. On the basis of this and previous studies, we expect the cyclic tetramers and pentamers to be the most significant cyclic water clusters in the atmosphere.