A theoretical study of aqueous solvation of K comparing ab initio, polarizable, and fixed-charge models

The hydration of K(+) is studied using a hierarchy of theoretical approaches, including ab initio Born-Oppenheimer molecular dynamics and Car-Parrinello molecular dynamics, a polarizable force field model based on classical Drude oscillators, and a nonpolarizable fixed-charge potential based on the TIP3P water model. While models based more directly on quantum mechanics offer the possibility to account for complex electronic effects, polarizable and fixed-charges force fields allow for simulations of large systems and the calculation of thermodynamic observables with relatively modest computational costs. A particular emphasis is placed on investigating the sensitivity of the polarizable model to reproduce key aspects of aqueous K(+), such as the coordination structure, the bulk hydration free energy, and the self diffusion of K(+). It is generally found that, while the simple functional form of the polarizable Drude model imposes some restrictions on the range of properties that can simultaneously be fitted, the resulting hydration structure for aqueous K(+) agrees well with experiment and with more sophisticated computational models. A counterintuitive result, seen in Car-Parrinello molecular dynamics and in simulations with the Drude polarizable force field, is that the average induced molecular dipole of the water molecules within the first hydration shell around K(+) is slightly smaller than the corresponding value in the bulk. In final analysis, the perspective of K(+) hydration emerging from the various computational models is broadly consistent with experimental data, though at a finer level there remain a number of issues that should be resolved to further our ability in modeling ion hydration accurately.     

Polarizable molecular dynamics simulation of Zn(II) in water using the AMOEBA force field

The hydration free energy, structure, and dynamics of the zinc divalent cation are studied using a polarizable force field in molecular dynamics simulations. Parameters for the Zn(2+) are derived from gas-phase ab initio calculation of Zn(2+)-water dimer. The Thole-based dipole polarization is adjusted based on the Constrained Space Orbital Variations (CSOV) calculation while the Symmetry Adapted Perturbation Theory (SAPT) approach is also discussed. The vdW parameters of Zn(2+) have been obtained by comparing the AMOEBA Zn(2+)-water dimerization energy with results from several theory levels and basis sets over a range of distances. Molecular dynamics simulations of Zn(2+) solvation in bulk water are subsequently performed with the polarizable force field. The calculated first-shell water coordination number, water residence time and free energy of hydration are consistent with experimental and previous theoretical values. The study is supplemented with extensive Reduced Variational Space (RVS) and Electron Localization Function (ELF) computations in order to unravel the nature of the bonding in Zn(2+)(H(2)O)(n) (n=1,6) complexes and to analyze the charge transfer contribution to the complexes. Results show that the importance of charge transfer decreases as the size of Zn-water cluster grows due to anticooperativity and to changes in the nature of the metal-ligand bonds. Induction could be dominated by polarization when the system approaches condensed-phase and the covelant effects are eliminated from the Zn(II)-water interaction. To construct an "effective" classical polarizable potential for Zn(2+) in bulk water, one should therefore avoid over-fitting to the ab initio charge transfer energy of Zn(2+)-water dimer. Indeed, in order to avoid overestimation of condensed-phase many-body effects, which is crucial to the transferability of polarizable molecular dynamics, charge transfer should not be included within the classical polarization contribution and should preferably be either incorporated in to the pairwise van der Waals contribution or treated explicitly.     

Ab initio molecular dynamics using hybrid density functionals

Ab initio molecular dynamics simulations with hybrid density functionals have so far found little application due to their computational cost. In this work, an implementation of the Hartree-Fock exchange is presented that is specifically targeted at ab initio molecular dynamics simulations of medium sized systems. We demonstrate that our implementation, which is available as part of the CP2K/Quickstep program, is robust and efficient. Several prescreening techniques lead to a linear scaling cost for integral evaluation and storage. Integral compression techniques allow for in-core calculations on systems containing several thousand basis functions. The massively parallel implementation respects integral symmetry and scales up to hundreds of CPUs using a dynamic load balancing scheme. A time-reversible multiple time step scheme, exploiting the difference in computational efficiency between hybrid and local functionals, brings further time savings. With extensive simulations of liquid water, we demonstrate the ability to perform, for several tens of picoseconds, ab initio molecular dynamics based on hybrid functionals of systems in the condensed phase containing a few thousand Gaussian basis functions.     

CHARMM fluctuating charge force field for proteins: II protein/solvent properties from molecular dynamics simulations using a nonadditive electrostatic model

A fluctuating charge (FQ) force field is applied to molecular dynamics simulations for six small proteins in explicit polarizable solvent represented by the TIP4P-FQ potential. The proteins include 1FSV, 1ENH, 1PGB, 1VII, 1H8K, and 1CRN, representing both helical and beta-sheet secondary structural elements. Constant pressure and temperature (NPT) molecular dynamics simulations are performed on time scales of several nanoseconds, the longest simulations yet reported using explicitly polarizable all-atom empirical potentials (for both solvent and protein) in the condensed phase. In terms of structure, the FQ force field allows deviations from native structure up to 2.5 A (with a range of 1.0 to 2.5 A). This is commensurate to the performance of the CHARMM22 nonpolarizable model and other currently existing polarizable models. Importantly, secondary structural elements maintain native structure in general to within 1 A (both helix and beta-strands), again in good agreement with the nonpolarizable case. In qualitative agreement with QM/MM ab initio dynamics on crambin (Liu et al. Proteins 2001, 44, 484), there is a sequence dependence of average condensed phase atomic charge for all proteins, a dependence one would anticipate considering the differing chemical environments around individual atoms; this is a subtle quantum mechanical feature captured in the FQ model but absent in current state-of-the-art nonpolarizable models. Furthermore, there is a mutual polarization of solvent and protein in the condensed phase. Solvent dipole moment distributions within the first and second solvation shells around the protein display a shift towards higher dipole moments (increases on the order of 0.2-0.3 Debye) relative to the bulk; protein polarization is manifested via the enhanced condensed phase charges of typical polar atoms such as backbone carbonyl oxygens, amide nitrogens, and amide hydrogens. Finally, to enlarge the sample set of proteins, gas-phase minimizations and 1 ps constant temperature simulations are performed on various-sized proteins to compare to earlier work by Kaminsky et al. (J Comp Chem 2002, 23, 1515). The present work establishes the feasibility of applying a fully polarizable force field for protein simulations and demonstrates the approach employed in extending the CHARMM force field to include these effects.     

Calibration of force-field dependency in free energy landscapes of peptide conformations by quantum chemical calculations

The free energy landscapes of peptide conformations were calibrated by ab initio quantum chemical calculations, after the enhanced conformational diversity search using the multicanonical molecular dynamics simulations. Three different potentials of mean force for an isolated dipeptide were individually obtained by the multicanonical molecular dynamics simulations using the conventional force fields, AMBER parm94, AMBER parm96, and CHARMm22. Each potential of mean force was then calibrated based upon the umbrella sampling algorithm from the adiabatic energy map that was calculated separately by the ab initio molecular orbital method, and all of the calibrated potentials of mean force coincided well. The calibration method was also applied to the simulations of a peptide dimer in explicit water models, and it was shown that the calibrated free energy landscapes did not depend on the force field used in the classical simulations, as far as the conformational space was sampled well. The current calibration method fuses the classical free energy calculation with the quantum chemical calculation, and it should generally make simulations for biomolecular systems much more reliable when combining with enhanced conformational sampling.     

Developing ab initio quality force fields from condensed phase quantum-mechanics/molecular-mechanics calculations through the adaptive force matching method

A new method called adaptive force matching (AFM) has been developed that is capable of producing high quality force fields for condensed phase simulations. This procedure involves the parametrization of force fields to reproduce ab initio forces obtained from condensed phase quantum-mechanics/molecular-mechanics (QM/MM) calculations. During the procedure, the MM part of the QM/MM is iteratively improved so as to approach ab initio quality. In this work, the AFM method has been tested to parametrize force fields for liquid water so that the resulting force fields reproduce forces calculated using the ab initio MP2 and the Kohn-Sham density functional theory with the Becke-Lee-Yang-Parr (BLYP) and Becke three-parameter LYP (B3LYP) exchange correlation functionals. The AFM force fields generated in this work are very simple to evaluate and are supported by most molecular dynamics (MD) codes. At the same time, the quality of the forces predicted by the AFM force fields rivals that of very expensive ab initio calculations and are found to successfully reproduce many experimental properties. The site-site radial distribution functions (RDFs) obtained from MD simulations using the force field generated from the BLYP functional through AFM compare favorably with the previously published RDFs from Car-Parrinello MD simulations with the same functional. Technical aspects of AFM such as the optimal QM cluster size, optimal basis set, and optimal QM method to be used with the AFM procedure are discussed in this paper.     

Behavior of the Eigen form of hydronium at the air/water interface

Surface affinity of hydro-nium was explored using umbrella sampling molecular dynamics simulations with a refined polarizable potential. The polarizable interaction potential of H(3)O(+) was reparametrized against accurate ab initio calculations for geometries including a water molecule approaching the Eigen cation from its oxygen side. Although there is no true hydrogen bonding with H(3)O(+) acting as an acceptor, respecting in the force field the very shallow ab initio minimum corresponding to this interaction leads to a decrease in surface propensity of hydronium compared to previous results. Qualitatively, the mild surface affinity and strong surface orientation of hydronium is, nevertheless, robustly predicted by various computational approaches, as well as by spectroscopic experiments.     

Gradients of the polarization energy in the effective fragment potential method

The effective fragment potential (EFP) method is an ab initio based polarizable classical method in which the intermolecular interaction parameters are obtained from preparative ab initio calculations on isolated molecules. The polarization energy in the EFP method is modeled with asymmetric anisotropic dipole polarizability tensors located at the centroids of localized bond and lone pair orbitals of the molecules. Analytic expressions for the translational and rotational gradients (forces and torques) of the EFP polarization energy have been derived and implemented. Periodic boundary conditions (the minimum image convention) and switching functions have also been implemented for the polarization energy, as well as for other EFP interaction terms. With these improvements, molecular dynamics simulations can be performed with the EFP method for various chemical systems.     

Towards accurate solvation dynamics of divalent cations in water using the polarizable amoeba force field: From energetics to structure

Molecular dynamics simulations were performed using a modified amoeba force field to determine hydration and dynamical properties of the divalent cations Ca2+ and Mg2+. The extension of amoeba to divalent cations required the introduction of a cation specific parametrization. To accomplish this, the Thole polarization damping model parametrization was modified based on the ab initio polarization energy computed by a constrained space orbital variation energy decomposition scheme. Excellent agreement has been found with condensed phase experimental results using parameters derived from gas phase ab initio calculations. Additionally, we have observed that the coordination of the calcium cation is influenced by the size of the periodic water box, a recurrent issue in first principles molecular dynamics studies.     

Embedded,graph-theoretically defined many-body approximations for wavefunction-in-DFT and DFT-in-DFT: Applications to gas- and condensed-phase ab initio molecular dynamics,and potential surfaces for quantum nuclear effects

We present a graph-theoretic approach to adaptively compute many-body approximations in an efficient manner to perform (a) accurate post-Hartree–Fock (HF) ab initio molecular dynamics (AIMD) at density functional theory (DFT) cost for medium- to large-sized molecular clusters, (b) hybrid DFT electronic structure calculations for condensed-phase simulations at the cost of pure density functionals, (c) reduced-cost on-the-fly basis extrapolation for gas-phase AIMD and condensed phase studies, and (d) accurate post-HF-level potential energy surfaces at DFT cost for quantum nuclear effects. The salient features of our approach are ONIOM-like in that (a) the full system (cluster or condensed phase) calculation is performed at a lower level of theory (pure DFT for condensed phase or hybrid DFT for molecular systems), and (b) this approximation is improved through a correction term that captures all many-body interactions up to any given order within a higher level of theory (hybrid DFT for condensed phase; CCSD or MP2 for cluster), combined through graph-theoretic methods. Specifically, a region of chemical interest is coarse-grained into a set of nodes and these nodes are then connected to form edges based on a given definition of local envelope (or threshold) of interactions. The nodes and edges together define a graph, which forms the basis for developing the many-body expansion. The methods are demonstrated through (a) ab initio dynamics studies on protonated water clusters and polypeptide fragments, (b) potential energy surface calculations on one-dimensional water chains such as those found in ion channels, and (c) conformational stabilization and lattice energy studies on homogeneous and heterogeneous surfaces of water with organic adsorbates using two-dimensional periodic boundary conditions.     

Ab initio molecular dynamics study of formate ion hydration

We perform ab initio molecular dynamics simulations of the aqueous formate ion. The mean number of water molecules in the first solvation shell, or the hydration number, of each formate oxygen is found to be consistent with recent experiments. Our ab initio pair correlation functions, however, differ significantly from many classical force field results and hybrid quantum mechanics/molecular mechanics predictions. They yield roughly one less hydrogen bond between each formate oxygen and water than force field or hybrid methods predict. Both the BLYP and PW91 exchange correlation functionals give qualitatively similar results. The time dependence of the hydration numbers are examined, and Wannier function techniques are used to analyze electronic configurations along the molecular dynamics trajectory.     

First principles and classical molecular dynamics simulations of solvated benzene

We have performed extensive ab initio and classical molecular dynamics (MD) simulations of benzene in water in order to examine the unique solvation structures that are formed. Qualitative differences between classical and ab initio MD simulations are found and the importance of various technical simulation parameters is examined. Our comparison indicates that nonpolarizable classical models are not capable of describing the solute-water interface correctly if local interactions become energetically comparable to water hydrogen bonds. In addition, a comparison is made between a rigid water model and fully flexible water within ab initio MD simulations which shows that both models agree qualitatively for this challenging system.     

Surface relaxation in water clusters: evidence from theoretical analysis of the oxygen 1s photoelectron spectrum

We present a theoretical interpretation of the oxygen 1s photoelectron spectrum published by Ohrwall et al. [J. Chem. Phys. 123, 054310 (2005)]. A water cluster that contains 200 molecules was simulated at 215 K using the polarizable AMOEBA force field. The force field predicts longer O...O distances at the cluster surface than in the bulk. Comparisons to ab initio molecular dynamics (MD) simulations indicate that the force field underestimates the degree of surface relaxation. By comparing cluster lineshape models, computed from MD simulations, to the experimental spectrum we find further evidence of surface relaxation.     

Comparing polarizable force fields to ab initio calculations reveals nonclassical effects in condensed phases

In a recent work [Giese and York J. Chem. Phys. 120, 9903 (2004)] showed that many-body force field models based solely on pairwise Coulomb screening cannot simultaneously reproduce both gas-phase and condensed-phase polarizability limits. In particular, polarizable force fields applied to bifurcated water chains have been demonstrated to be overpolarized with respect to ab initio methods. This behavior was ascribed to the neglect of coupling between many-body exchange and polarization. In the present article we reproduce those results using different ab initio levels of theory and a polarizable model based on the chemical-potential equalization principle. Moreover we show that, when hydrogen-bond (H-bond) forming systems are considered, an additional nonclassical effect, i.e., intermolecular charge transfer, must be taken into account. Such effect, contrarily to that of coupling between many-body exchange and polarization, makes classical polarizable force fields underpolarized. In the case of water at standard conditions, being H-bonded geometries much more probable than the bifurcated ones, intermolecular charge transfer is the dominant effect.     

Ab initio parametrized polarizable force field for rutile-type SnO<Subscript>2</Subscript>

We report a new, polarizable classical force field for the rutile-type phase of SnO₂, casserite. This force field has been parametrized using results from ab initio (density functional theory) calculations as a basis for fitting. The force field was found to provide structural, dynamical and thermodynamic properties of tin oxide that compare well with both ab initio and experimental results at ambient and high pressures.     

Urea: an ab initio and force field study of the gas and solid phases

We have studied the gaseous and solid phases of urea using both quantum mechanics calculation and force field simulation methods. Our ab initio calculations confirmed experimental observations that urea structure is planar in the crystal, but nonplanar in the gas phase. Based on electron structure analysis, we suggest that the significant difference between these two structures in different environments can be qualitatively explained by two resonance structures. The planar structure is more polarized than the nonplanar one, and the former is stabilized in the solid phases due to strong electrostatic interactions. We found classical force field method is incapable to represent such strong polarization effect. Using molecular dynamics simulations with a force field optimized for condensed phases, we calculated the crystalline structures of urea in the temperature range of 12 to 293 K. The densities as well as cell parameters are within 2% deviation from the experimental data in the temperature range.