The flexible, polarizable, thole-type interaction potential for water (TTM2-F) revisited

We present a revision of the flexible, polarizable, Thole-type interaction potential for water [J. Chem. Phys.2002, 116, 5115], which allows for condensed-phase simulations. The revised version (TTM2.1-F) of the potential correctly describes the individual water molecular dipole moment and alleviates problems arising at short intermolecular separations that can be sampled in the course of molecular dynamics and Monte Carlo simulations of condensed environments. Furthermore, its parallel implementation under periodic boundary conditions enables the efficient calculation of the macroscopic structural and thermodynamic properties of liquid water, as its performance scales superlinearly with up to a number of 64 processors for a simulation box of 512 molecules. We report the radial distribution functions, average energy, internal geometry, and dipole moment in the liquid as well as the density, dielectric constant, and self-diffusion coefficient at T = 300 K from (NVT) and (NPT) classical molecular dynamics simulations by using the revised version of the potential.     

Re-examining the properties of the aqueous vapor-liquid interface using dispersion corrected density functional theory

First-principles molecular dynamics simulations, in which the forces are computed from electronic structure calculations, have great potential to provide unique insight into structure, dynamics, electronic properties, and chemistry of interfacial systems that is not available from empirical force fields. The majority of current first-principles simulations are driven by forces derived from density functional theory with generalized gradient approximations to the exchange-correlation energy, which do not capture dispersion interactions. We have carried out first-principles molecular dynamics simulations of air-water interfaces employing a particular generalized gradient approximation to the exchange-correlation functional (BLYP), with and without empirical dispersion corrections. We assess the utility of the dispersion corrections by comparison of a variety of structural, dynamic, and thermodynamic properties of bulk and interfacial water with experimental data, as well as other first-principles and force field-based simulations.     

Ab initio molecular dynamics simulations of the liquid/vapor interface of sulfuric acid solutions

Ab initio molecular dynamics simulations of the liquid-vapor interface are presented for thin slabs of 72 water molecules containing a single molecule of sulfuric acid. Trajectories in the 306-330 K range are calculated for two functionals with double- and triple-ζ quality basis sets. Comparisons are made between BLYP and HCTH/120 results for the slab simulations and for bulk simulations of one H(2)SO(4) in a periodic box with 63 waters. Good agreement is found with the available experimental data and the results of other relevant AIMD studies with respect to ionization of the acid, size of the coordination shells, partitioning of the ions with the hydronium exhibiting a surface preference and the anions in the interior, and the orientational distributions for the hydronium ions and for the surface/subsurface water molecules. The major differences in the performance of the two functionals are attributable to the greater basicity of the anion oxygen atoms with the HCTH functional and the more structured aqueous solution with BLYP. The enhanced basicity results in larger aqueous coordination shells for the anion oxygens. The structuring of the BLYP aqueous solution is observed in the corrugation of the water density profile, the higher first peak in g(OO)(r), and a smaller water self-diffusion constant. This structuring with the BLYP functional yields anion hydrogen bonds that endure longer and where the dissociated ions more rapidly and directly segregate in the slab. The simulations indicate that aqueous surfaces containing ionizable diprotic acids can be modeled with rather modest sized systems and be informative.     

An efficient parallelization scheme for molecular dynamics simulations with many-body, flexible, polarizable empirical potentials: application to water

An efficient parallelization scheme for classical molecular dynamics simulations with flexible, polarizable empirical potentials is presented. It is based on the standard Ewald summation technique to handle the long-range electrostatic and induction interactions. The algorithm for this parallelization scheme is designed for systems containing several thousands of polarizable sites in the simulation box. Its performance is evaluated during molecular dynamics simulations under periodic boundary conditions with unit cell sizes ranging from 128 to 512 molecules employing two flexible polarizable water models [DC(F) and TTM2.1-F] containing 1 and 3 polarizable sites, respectively. The time-to-solution for these two polarizable models is compared with the one for a flexible, pairwise-additive water model (TIP4F). The benchmarks were performed on both shared and distributed memory platforms. As a result of the efficient calculation of the induced dipole moments, a superlinear scaling as a function of the number of the processors is observed. To the best of our knowledge, this is the first reported results of parallel scaling and performance for simulations of liquid water with a polarizable potential under periodic boundary conditions.     

Comparisons of classical chemical dynamics simulations of the unimolecular decomposition of classical and quantum microcanonical ensembles

Previous studies have shown that classical trajectory simulations often give accurate results for short-time intramolecular and unimolecular dynamics, particularly for initial non-random energy distributions. To obtain such agreement between experiment and simulation, the appropriate distributions must be sampled to choose initial coordinates and momenta for the ensemble of trajectories. If a molecule's classical phase space is sampled randomly, its initial decomposition will give the classical anharmonic microcanonical (RRKM) unimolecular rate constant for its decomposition. For the work presented here, classical trajectory simulations of the unimolecular decomposition of quantum and classical microcanonical ensembles, at the same fixed total energy, are compared. In contrast to the classical microcanonical ensemble, the quantum microcanonical ensemble does not sample the phase space randomly. The simulations were performed for CH(4), C(2)H(5), and Cl(-)---CH(3)Br using both analytic potential energy surfaces and direct dynamics methods. Previous studies identified intrinsic RRKM dynamics for CH(4) and C(2)H(5), but intrinsic non-RRKM dynamics for Cl(-)---CH(3)Br. Rate constants calculated from trajectories obtained by the time propagation of the classical and quantum microcanonical ensembles are compared with the corresponding harmonic RRKM estimates to obtain anharmonic corrections to the RRKM rate constants. The relevance and accuracy of the classical trajectory simulation of the quantum microcanonical ensemble, for obtaining the quantum anharmonic RRKM rate constant, is discussed.     

Time-dependent density functional theory description of on-site electron repulsion and ligand field effects in the optical spectrum of hexaaquoruthenium(II) in solution

Time-dependent density functional theory (TDDFT) and density functional-based molecular dynamics were used to simulate the finite temperature t(2g)5 e(g) <-- t(2g)6 absorption band of the Ru2+ hexahydrate coordination complex in aqueous solution. The (1)T1 <-- (1)A1 and (1)T2 <-- (1)A1 molecular term splitting of this transition, which is not accounted for by the Kohn-Sham excitation spectrum, is shown to be satisfactorily reproduced by TDDFT at the BLYP/ALDA level of theory. Comparison to the spectrum of the Ru2+ (H2O)6 complex in vacuo computed by similar density functional classical molecular dynamics methods leads to the observation that bulk solvation has a negligible effect on the position and the shape of the absorption profile.     

The influence of temperature and density functional models in ab initio molecular dynamics simulation of liquid water

The performance of density functional theory methods for the modeling of condensed aqueous systems is hard to predict and validation by ab initio molecular simulation of liquid water is absolutely necessary. In order to assess the reliability of these tests, the effect of temperature on the structure and dynamics of liquid water has been characterized with 16 simulations of 20 ps in the temperature range of 280-380 K. We find a pronounced influence of temperature on the pair correlation functions and on the diffusion constant including nonergodic behavior on the time scale of the simulation in the lower temperature range (which includes ambient temperature). These observations were taken into account in a consistent comparison of a series of density functionals (BLYP, PBE, TPSS, OLYP, HCTH120, HCTH407). All simulations were carried out using an ab initio molecular dynamics approach in which wave functions are represented using Gaussians and the density is expanded in an auxiliary basis of plane waves. Whereas the first three functionals show similar behavior, it is found that the latter three functionals yield more diffusive dynamics and less structure.     

Ab initio calculations of the melting temperatures of refractory bcc metals

We present ab initio calculations of the melting temperatures for bcc metals Nb, Ta and W. The calculations combine phase coexistence molecular dynamics (MD) simulations using classical embedded-atom method potentials and ab initio density functional theory free energy corrections. The calculated melting temperatures for Nb, Ta and W are, respectively, within 3%, 4%, and 7% of the experimental values. We compare the melting temperatures to those obtained from direct ab initio molecular dynamics simulations and see if they are in excellent agreement with each other. The small remaining discrepancies with experiment are thus likely due to inherent limitations associated with exchange-correlation energy approximations within density-functional theory.     

Car-Parrinello molecular dynamics simulation of Fe 3+ (aq)

The optimized geometry and energetic properties of Fe(D2O)n 3+ clusters, with n = 4 and 6, have been studied with density-functional theory calculations and the BLYP functional, and the hydration of a single Fe 3+ ion in a periodic box with 32 water molecules at room temperature has been studied with Car-Parrinello molecular dynamics and the same functional. We have compared the results from the CPMD simulation with classical MD simulations, using a flexible SPC-based water model and the same number of water molecules, to evaluate the relative strengths and weaknesses of the two MD methods. The classical MD simulations and the CPMD simulations both give Fe-water distances in good agreement with experiment, but for the intramolecular vibrations, the classical MD yields considerably better absolute frequencies and ion-induced frequency shifts. On the other hand, the CPMD method performs considerably better than the classical MD in describing the intramolecular geometry of the water molecule in the first hydration shell and the average first shell...second shell hydrogen-bond distance. Differences between the two methods are also found with respect to the second-shell water orientations. The effect of the small box size (32 vs 512 water molecules) was evaluated by comparing results from classical simulations using different box sizes; non-negligible effects are found for the ion-water distance and the tilt angles of the water molecules in the second hydration shell and for the O-D stretching vibrational frequencies of the water molecules in the first hydration shell.     

Vibrational relaxation of azide ions in liquid-to-supercritical water

The dynamics of vibrational energy relaxation (VER) of the aqueous azide anion was studied over a wide temperature (300 K ≤ T ≤ 663 K) and density (0.6 g cm(-3) ≤ ρ ≤ 1.0 g cm(-3)) range thereby covering the liquid and the supercritical phase of the water solvent. Femtosecond mid-infrared spectroscopy on the ν(3) band associated with the asymmetric stretching vibration of the azide anion was used to monitor the relaxation dynamics in a time-resolved fashion. The variation of the vibrational relaxation rate constant with temperature and density was found to be rather small. Surprisingly, the simple isolated binary collision model is able to fully reproduce the experimentally observed temperature and density dependence of the relaxation rate provided a local density correction around the vibrationally excited solute based on classical molecular dynamics simulations is used. The simulations further suggest that head-on collisions of the solvent with the terminal nitrogen atoms rather than side-on collisions with the central nitrogen atom of the azide govern the vibrational energy relaxation of this system. Finally, the importance of hydrogen bonding for the VER dynamics in this system is briefly discussed.     

The effects of electronic polarization on water adsorption in metal-organic frameworks: H(2)O in MIL-53(Cr)

The effects of electronic polarization on the adsorption of water in the MIL-53(Cr) metal-organic framework are investigated using molecular dynamics simulations. For this purpose a fully polarizable force field for MIL-53(Cr) was developed which is compatible with the ab initio-based TTM3-F water model. The analysis of the spatial distributions of the water molecules within the MIL-53(Cr) nanopores calculated as a function of loading indicates that polarization effects play an important role in the formation of hydrogen bonds between the water molecules and the hydroxyl groups of the framework. As a result, large qualitative differences are found between the radial distribution functions calculated with non-polarizable and polarizable force fields. The present analysis suggests that polarization effects can significantly impact molecular adsorption in metal-organic frameworks under hydrated conditions.     

An application of flexible constraints in Monte Carlo simulations of the isobaric--isothermal ensemble of liquid water and ice Ih with the polarizable and flexible mobile charge densities in harmonic oscillators model

The method of flexible constraints was implemented in a Monte Carlo code to perform numerical simulations of liquid water and ice Ih in the constant number of molecules, volume, and temperature and constant pressure, instead of volume ensembles, using the polarizable and flexible mobile charge densities in harmonic oscillators (MCDHO) model. The structural and energetic results for the liquid at T=298 K and rho=997 kg m(-3) were in good agreement with those obtained from molecular dynamics. The density obtained at P=1 atm with flexible constraints, rho=1008 kg m(-3), was slightly lower than with the classical sampling of the intramolecular vibrations, rho=1010 kg m(-3). The comparison of the structures and energies found for water hexamers and for ice Ih with six standard empirical models to those obtained with MCDHO, show this latter to perform better in describing water far from ambient conditions: the MCDHO minimum lattice energy, density, and lattice constants were in good agreement with experiment. The average angle HOH of the water molecule in ice was predicted to be slightly larger than in the liquid, yet 1.2% smaller than the experimental value.     

平面五元水分子簇的量子化学研究

选用Gaussian03的B3LYP/6-31G(d,p)、DMol3的BLYP/DNP和deMon的BLYP/TZVP等方法计算了甲烷水合物(结构-1)中平面五元水分子簇的结合能和氢键能,作了基组重叠误差(BSSE)和色散能(dispersion)的修正,估算了次级相互作用的贡献.在DMol3程序中使用了大型数值基组DNP,将基组重叠误差降至最低.在Gaussi-an03的B3LYP/6-31G(d,p)计算中,采用平衡法(Counterpoise)校正基组重叠误差.两种计算方法给出了一致的结果,证实了在使用6-31G(d,p)基组时,一对水分子在平衡距离的基组重叠误差高达8 kJ/mol.为估算色散能的贡献,使用了新近发展的包含色散能的密度泛函的DFT程序deMon计算了五元水分子簇.用多种方法计算出了经基组重叠误差和色散能修正的五元水分子簇的分子间结合能和氢键能的较为精确的势能超曲面,为甲烷和其他气体水合物的分子动力学模拟提供了依据.     

Density functional theory study of uranium(VI) aquo chloro complexes in aqueous solution

Mixed uranyl aquo chloro complexes of the type [UO2(H2O)xCly]2-y (y = 1, 2, 3, 4; x + y = 4, 5) have been optimized at the BLYP, BP86, and B3LYP levels of density functional theory in vacuo and in a polarizable continuum modeling bulk water (PCM) and have been studied at the BLYP level with Car-Parrinello molecular dynamics (MD) simulations in the gas phase and in explicit aqueous solution. Free binding energies were evaluated from static PCM data and from pointwise thermodynamic integration involving constrained MD simulations in water. The computations reveal significant solvent effects on geometric and energetic parameters. Based on the comparison of PCM-optimized or MD-averaged uranyl-ligand bond distances with EXAFS-derived values, the transition between five- and four-coordination about uranyl is indicated to occur at a Cl content of y = 2 or 3.     

Elucidating the vibrational spectra of hydrogen-bonded aggregates in solution: electronic structure calculations with implicit solvent and first-principles molecular dynamics simulations with explicit solvent for 1-hexanol in n-hexane

Fourier transform infrared spectroscopy is a popular method for the experimental investigation of hydrogen-bonded aggregates, but linking spectral information to microscopic information on aggregate size distribution and aggregate architecture is an arduous task. Static electronic structure calculations with an implicit solvent model, Car-Parrinello molecular dynamics (CPMD) using the Becke-Lee-Yang-Parr (BLYP) exchange and correlation energy functionals and classical molecular dynamics simulations for the all-atom version of the optimized parameters for liquid simulations (OPLS-AA) force field were carried out for an ensemble of 1-hexanol aggregates solvated in n-hexane. The initial configurations for these calculations were size-selected from a distribution of aggregates obtained from a large-scale Monte Carlo simulation. The vibrational spectra computed from the static electronic structure calculations for monomers and dimers and from the CPMD simulations for aggregates up to pentamers demonstrate the extent of the contribution of dangling or nondonating hydroxyl groups found in linear and branched aggregates to the "monomeric" peak. Furthermore, the computed spectra show that there is no simple relationship between peak shift and aggregate size nor architecture, but the effect of hydrogen-bond cooperativity is shown to differentiate polymer-like (cooperative) and dimer-like (noncooperative) hydrogen bonds in the vibrational spectrum. In contrast to the static electronic structure calculations and the CPMD simulations, the classical molecular dynamics simulations greatly underestimate the vibrational peak shift due to hydrogen bonding.     

Benchmark ab initio proton affinity and gas-phase basicity of α-alanine based on coupled-cluster theory and statistical mechanics

We determine the proton affinity (PA) and gas-phase basicity (GB) of amino acid α-alanine at a chemically accurate level by performing explicitly-correlated CCSD(T)-F12b/aug-cc-pVDZ geometry optimizations and normal mode vibrational frequency calculations as well as CCSD(T)-F12b/aug-cc-pVTZ energy computations at the possible neutral and protonated geometries. Temperature effects at 298.15 K considering translational, rotational, and vibrational enthalpy and entropy corrections are obtained via standard statistical mechanics utilizing the molecular geometries and the harmonic vibrational energy levels. Both the amino nitrogen (N) and the carbonyl oxygen (O) atoms are proven to be potential protonation sites and a systematic conformational search reveals 3 N- and 9 O-protonated conformers in the 0.00–7.88 and 25.43–30.43 kcal/mol energy ranges at 0 K, respectively. The final computed PA and GB values at (0)298.15 K in case of N-protonation are (214.47)216.80 and 207.07 kcal/mol, respectively, whereas the corresponding values for O-protonation are (189.04)190.63 and 182.31 kcal/mol. The results of the benchmark high-level coupled-cluster computations are utilized to assess the accuracy of several lower-level cost-effective methods such as MP2 and density functional theory with various functionals (SOGGA11-X, M06-2X, PBE0, B3LYP, M06, TPSS).     

Theoretical characterization of the (H2O)21 cluster: application of an n-body decomposition procedure

Two low-energy minima of (H2O)21 with very different H-bonding arrangements have been investigated with the B3LYP density functional and RIMP2 methods, as well as with the TIP4P, Dang-Chang, AMOEBA, and TTM2-F force fields. The AMOEBA and TTM2-F model potentials give an energy ordering that agrees with the results of the electronic structure calculations, while the TIP4P and Dang-Chang models give the opposite ordering. Insight into the role of many-body polarization for establishing the relative stability of the two isomers is provided by an n-body decomposition of the energies calculated using the various theoretical methods.     

Hydronium ion complex of 18-crown-6: where are the protons? A density functional study of static and dynamic properties

A quantum-chemical study employing the BLYP density functional is reported for the complex of H3O+ with 18-crown-6. According to a Car-Parrinello molecular dynamics (CPMD) study at 340 K, the complex is quite flexible, and is characterized by three quasi-linear (two-center) hydrogen-bond interactions for most of the time. On a time scale of 10 ps, frequent inversions of H3O+ are observed, as well as two 120 degrees rotations switching the hydrogen bonds from one set of crown-ether O atoms to the other. These results are consistent with density-functional studies of stationary points on the potential energy surface, which show how the crown "catalyzes" the guest's inversion. Two close-lying minima are characterized, as well as two distinct transition states connecting them, either via H3O+ inversion or rotation, with barriers of 1.0 and 4.6 kcal/mol, respectively, at the BLYP/II'//BLYP/6-31G level. Orbital interactions between lone pairs on ether O atoms and hydronium sigma(OH) antibonding orbitals are important factors for the directionality of the hydrogen bonds.