Weakly bound water structure, bond valence saturation and water dynamics at the goethite (100) surface/aqueous interface: ab initio dynamical simulations.

BACKGROUND: Many important geochemical and biogeochemical reactions occur in the mineral/formation water interface of the highly abundant mineral, goethite [α-Fe(OOH)]. Ab initio molecular dynamics (AIMD) simulations of the goethite α-FeOOH (100) surface and the structure, water bond formation and dynamics of water molecules in the mineral/aqueous interface are presented. Several exchange correlation functionals were employed (PBE96, PBE96 + Grimme, and PBE0) in the simulations of a (3 × 2) goethite surface with 65 absorbed water molecules in a 3D-periodic supercell (a = 30 Å, FeOOH slab ~12 Å thick, solvation layer ~18 Å thick). RESULTS: The lowest energy goethite (100) surface termination model was determined to have an exposed surface Fe3+ that was loosely capped by a water molecule and a shared hydroxide with a neighboring surface Fe3+. The water molecules capping surface Fe3+ ions were found to be loosely bound at all DFT levels with and without Grimme corrections, indicative that each surface Fe3+ was coordinated with only five neighbors. These long bonds were supported by bond valence theory calculations, which showed that the bond valence of the surface Fe3+ was saturated and surface has a neutral charge. The polarization of the water layer adjacent to the surface was found to be small and affected only the nearest water. Analysis by density difference plots and localized Boys orbitals identified three types of water molecules: those loosely bound to the surface Fe3+, those hydrogen bonded to the surface hydroxyl, and bulk water with tetrahedral coordination. Boys orbital analysis showed that the spin down lone pair orbital of the weakly absorbed water interact more strongly with the spin up Fe3+ ion. These weakly bound surface water molecules were found to rapidly exchange with the second water layer (~0.025 exchanges/ps) using a dissociative mechanism. CONCLUSIONS: Water molecules adjacent to the surface were found to only weakly interact with the surface and as a result were readily able to exchange with the bulk water. To account for the large surface Fe-OH2 distances in the DFT calculations it was proposed that the surface Fe3+ atoms, which already have their bond valence fully satisfied with only five neighbors, are under-coordinated with respect to the bulk coordination. Graphical abstract All first principle calculations, at all practically achievable levels, for the goethite 100 aqueous interface support a long bond and weak interaction between the exposed surface Fe3+ and water molecules capping the surface. This result is supported by bond valence theory calculations and is indicative that each surface Fe3+ is coordinated with only 5 neighbors.

QM/MM simulation (B3LYP) of the RNase A cleavage-transesterification reaction supports a triester A(N) + D(N) associative mechanism with an O2' H internal proton transfer.

The mechanism of the backbone cleavage-transesterification step of the RNase A enzyme remains controversial even after 60 years of study. We report quantum mechanics/molecule mechanics (QM/MM) free energy calculations for two optimized reaction paths based on an analysis of all structural data and identified by a search for reaction coordinates using a reliable quantum chemistry method (B3LYP), equilibrated structural optimizations, and free energy estimations. Both paths are initiated by nucleophilic attack of the ribose O2' oxygen on the neighboring diester phosphate bond, and both reach the same product state (PS) (a O3'-O2' cyclic phosphate and a O5' hydroxyl terminated fragment). Path 1, resembles the widely accepted dianionic transition-state (TS) general acid (His119)/base (His12) classical mechanism. However, this path has a barrier (25 kcal/mol) higher than that of the rate-limiting hydrolysis step and a very loose TS. In Path 2, the proton initially coordinating the O2' migrates to the nonbridging O1P in the initial reaction path rather than directly to the general base resulting in a triester (substrate as base) AN + DN mechanism with a monoanionic weakly stable intermediate. The structures in the transition region are associative with low barriers (TS1 10, TS2 7.5 kcal/mol). The Path 2 mechanism is consistent with the many results from enzyme and buffer catalyzed and uncatalyzed analog reactions and leads to a PS consistent with the reactive state for the following hydrolysis step. The differences between the consistently estimated barriers in Path 1 and 2 lead to a 10(11) difference in rate strongly supporting the less accepted triester mechanism.

Extending molecular simulation time scales: Parallel in time integrations for high-level quantum chemistry and complex force representations.

Parallel in time simulation algorithms are presented and applied to conventional molecular dynamics (MD) and ab initio molecular dynamics (AIMD) models of realistic complexity. Assuming that a forward time integrator, f (e.g., Verlet algorithm), is available to propagate the system from time ti (trajectory positions and velocities xi = (ri, vi)) to time ti + 1 (xi + 1) by xi + 1 = fi(xi), the dynamics problem spanning an interval from t0[ellipsis (horizontal)]tM can be transformed into a root finding problem, F(X) = [xi - f(x(i - 1)]i = 1, M = 0, for the trajectory variables. The root finding problem is solved using a variety of root finding techniques, including quasi-Newton and preconditioned quasi-Newton schemes that are all unconditionally convergent. The algorithms are parallelized by assigning a processor to each time-step entry in the columns of F(X). The relation of this approach to other recently proposed parallel in time methods is discussed, and the effectiveness of various approaches to solving the root finding problem is tested. We demonstrate that more efficient dynamical models based on simplified interactions or coarsening time-steps provide preconditioners for the root finding problem. However, for MD and AIMD simulations, such preconditioners are not required to obtain reasonable convergence and their cost must be considered in the performance of the algorithm. The parallel in time algorithms developed are tested by applying them to MD and AIMD simulations of size and complexity similar to those encountered in present day applications. These include a 1000 Si atom MD simulation using Stillinger-Weber potentials, and a HCl + 4H2O AIMD simulation at the MP2 level. The maximum speedup (serial execution/timeparallel execution time) obtained by parallelizing the Stillinger-Weber MD simulation was nearly 3.0. For the AIMD MP2 simulations, the algorithms achieved speedups of up to 14.3. The parallel in time algorithms can be implemented in a distributed computing environment using very slow transmission control protocol/Internet protocol networks. Scripts written in Python that make calls to a precompiled quantum chemistry package (NWChem) are demonstrated to provide an actual speedup of 8.2 for a 2.5 ps AIMD simulation of HCl + 4H2O at the MP2/6-31G* level. Implemented in this way these algorithms can be used for long time high-level AIMD simulations at a modest cost using machines connected by very slow networks such as WiFi, or in different time zones connected by the Internet. The algorithms can also be used with programs that are already parallel. Using these algorithms, we are able to reduce the cost of a MP2/6-311++G(2d,2p) simulation that had reached its maximum possible speedup in the parallelization of the electronic structure calculation from 32 s/time step to 6.9 s/time step.

Ion association in AlCl3 aqueous solutions from constrained first-principles molecular dynamics.

The Car-Parrinello-based molecular dynamics (CPMD) method was used to investigate the ion-pairing behavior between Cl(-) and Al(3+) ions in an aqueous AlCl(3) solution containing 63 water molecules. A series of constrained simulations was carried out at 300 K for up to 16 ps each, with the internuclear separation (r(Al-Cl)) between the Al(3+) ion and one of the Cl(-) ions held constant. The calculated potential of mean force (PMF) of the Al(3+)-Cl(-) ion pair shows a global minimum at r(Al-Cl) = 2.3 Å corresponding to a contact ion pair (CIP). Two local minima assigned to solvent-separated ion pairs (SSIPs) are identified at r(Al-Cl) = 4.4 and 6.0 Å. The positions of the free energy minima coincide with the hydration-shell intervals of the Al(3+) cation, suggesting that the Cl(-) ion is inclined to reside in regions with low concentrations of water molecules, that is, between the first and second hydration shells of Al(3+) and between the second shell and the bulk. A detailed analysis of the solvent structure around the Al(3+) and Cl(-) ions as a function of r(Al-Cl) is presented. The results are compared to structural data from X-ray measurements and unconstrained CPMD simulations of single Al(3+) and Cl(-) ions and AlCl(3) solutions. The dipole moments of the water molecules in the first and second hydration shells of Al(3+) and in the bulk region and those of Cl(-) ions were calculated as a function of r(Al-Cl). Major changes in the electronic structure of the system were found to result from the removal of Cl(-) from the first hydration shell of the Al(3+) cation. Finally, two unconstrained CPMD simulations of aqueous AlCl(3) solutions corresponding to CIP and SSIP configurations were performed (17 ps, 300 K). Only minor structural changes were observed in these systems, confirming their stability.

Near-Quantitative Agreement of Model-Free DFT-MD Predictions with XAFS Observations of the Hydration Structure of Highly Charged Transition-Metal Ions.

First-principles dynamics simulations (DFT, PBE96, and PBE0) and electron scattering calculations (FEFF9) provide near-quantitative agreement with new and existing XAFS measurements for a series of transition-metal ions interacting with their hydration shells via complex mechanisms (high spin, covalency, charge transfer, etc.). This analysis does not require either the development of empirical interparticle interaction potentials or structural models of hydration. However, it provides consistent parameter-free analysis and improved agreement with the higher-R scattering region (first- and second-shell structure, symmetry, dynamic disorder, and multiple scattering) for this comprehensive series of ions. DFT+GGA MD methods provide a high level of agreement. However, improvements are observed when exact exchange is included. Higher accuracy in the pseudopotential description of the atomic potential, including core polarization and reducing core radii, was necessary for very detailed agreement. The first-principles nature of this approach supports its application to more complex systems.

Parallel implementation of Γ-point pseudopotential plane-wave DFT with exact exchange.

Semi-local functionals commonly used in density functional theory (DFT) studies of solids usually fail to reproduce localized states such as trapped holes, polarons, excitons, and solitons. This failure is ascribed to self-interaction which creates a Coulomb barrier to localization. Pragmatic approaches in which the exchange correlation functionals are augmented with small amount of exact exchange (hybrid-DFT, e.g., B3LYP and PBE0) have shown to promise in rectifying this type of failure, as well as producing more accurate band gaps and reaction barriers. The evaluation of exact exchange is challenging for large, solid state systems with periodic boundary conditions, especially when plane-wave basis sets are used. We have developed parallel algorithms for implementing exact exchange into pseudopotential plane-wave DFT program and we have implemented them in the NWChem program package. The technique developed can readily be employed in Γ-point plane-wave DFT programs. Furthermore, atomic forces and stresses are straightforward to implement, making it applicable to both confined and extended systems, as well as to Car-Parrinello ab initio molecular dynamic simulations. This method has been applied to several systems for which conventional DFT methods do not work well, including calculations for band gaps in oxides and the electronic structure of a charge trapped state in the Fe(II) containing mica, annite.

Structure and dynamics of the hydration shells of the Zn(2+) ion from ab initio molecular dynamics and combined ab initio and classical molecular dynamics simulations.

Results of ab initio molecular dynamics (AIMD) simulations (density functional theory+PBE96) of the dynamics of waters in the hydration shells surrounding the Zn(2+) ion (T approximately 300 K, rho approximately 1 gm/cm(3)) are compared to simulations using a combined quantum and classical molecular dynamics [AIMD/molecular mechanical (MM)] approach. Both classes of simulations were performed with 64 solvating water molecules ( approximately 15 ps) and used the same methods in the electronic structure calculation (plane-wave basis set, time steps, effective mass, etc.). In the AIMD/MM calculation, only six waters of hydration were included in the quantum mechanical (QM) region. The remaining 58 waters were treated with a published flexible water-water interaction potential. No reparametrization of the water-water potential was attempted. Additional AIMD/MM simulations were performed with 256 water molecules. The hydration structures predicted from the AIMD and AIMD/MM simulations are found to agree in detail with each other and with the structural results from x-ray data despite the very limited QM region in the AIMD/MM simulation. To further evaluate the agreement of these parameter-free simulations, predicted extended x-ray absorption fine structure (EXAFS) spectra were compared directly to the recently obtained EXAFS data and they agree in remarkable detail with the experimental observations. The first hydration shell contains six water molecules in a highly symmetric octahedral structure is (maximally located at 2.13-2.15 A versus 2.072 A EXAFS experiment). The widths of the peak of the simulated EXAFS spectra agree well with the data (8.4 A(2) versus 8.9 A(2) in experiment). Analysis of the H-bond structure of the hydration region shows that the second hydration shell is trigonally bound to the first shell water with a high degree of agreement between the AIMD and AIMD/MM calculations. Beyond the second shell, the bonding pattern returns to the tetrahedral structure of bulk water. The AIMD/MM results emphasize the importance of a quantum description of the first hydration shell to correctly describe the hydration region. In these calculations the full d(10) electronic structure of the valence shell of the Zn(2+) ion is retained. The simulations show substantial and complex charge relocation on both the Zn(2+) ion and the first hydration shell. The dipole moment of the waters in the first hydration shell is 3.4 D (3.3 D AIMD/MM) versus 2.73 D bulk. Little polarization is found for the waters in the second hydration shell (2.8 D). No exchanges were seen between the first and the second hydrations shells; however, many water transfers between the second hydration shell and the bulk were observed. For 64 waters, the AIMD and AIMD/MM simulations give nearly identical results for exchange dynamics. However, in the larger particle simulations (256 waters) there is a significant reduction in the second shell to bulk exchanges.

Vertical ionization potentials of nucleobases in a fully solvated DNA environment.

Vertical ionization potentials (IPs) of nucleobases embedded in a fully solvated DNA fragment (12-mer B-DNA fragment + 22 sodium counterions + 5760 water molecules equilibrated to 298 K) have been calculated using a combined quantum mechanical molecular mechanics (QM/MM) approach. Calculations of the vertical IP of the anion Cl(-) are reported that support the accuracy of the application of a QM/MM method to this problem. It is shown that the pi nucleotide HOMO origin for the emitted electron is localized on the base by the hydration structure surrounding the DNA in a way similar to that recently observed for pyrimidine nucleotides in aqueous solutions (Slavicek, P.; et al. J. Am. Chem. Soc. 2009, 131, 6460). In a first step, a high level of theory, CCSD(T)/aug-cc-pVDZ, was used to calculate the vertical IP of each of the four single bases isolated in the QM region while the remaining DNA fragment, counterions, and water solvent molecules were included in the MM region. The calculated vertical IPs show a large positive shift of 3.2-3.3 eV compared to the corresponding gas-phase values. This shift is similar for all four DNA bases. The origin of the large increase in vertical IPs of nucleobases is found to be the long-range electrostatic interactions with the solvation structure outside the DNA helix. Thermal fluctuations in the fluid can result in IP changes of roughly 1 eV on a picosecond time scale. IPs of pi-stacked and H-bonded clusters of DNA bases were also calculated using the same QM/MM model but with a lower level of theory, B3LYP/6-31G(d=0.2). An IP shift of 4.02 eV relative to the gas phase is found for a four-base-pair B-DNA duplex configuration. The primary goal of this work was to estimate the influence of long-range solvation interactions on the ionization properties of DNA bases rather than provide highly precise IP evaluations. The QM/MM model presented in this work provides an attractive method to treat the difficult problem of incorporating a detailed long-range structural model of physiological conditions into investigations of the electronic processes in DNA.

First principles simulation of the bonding, vibrational, and electronic properties of the hydration shells of the high-spin Fe(3+) ion in aqueous solutions.

Results of parameter-free first principles simulations of a spin up 3d(5) Fe(3+) ion hydrated in an aqueous solution (64 waters, 30 ps, 300 K) are reported. The first hydration shell associated with the first maximum of the radial distribution function, g(FeO)(r), at d(Fe-O(I)) = 2.11-2.15 A, contains 6 waters with average d(OH) = 0.99 A, in good agreement with observations. A second shell with average coordination number 13.3 can be identified with average shell radius of d(Fe-O(II)) = 4.21-4.32 A. The waters in this hydration shell are coordinated to the first shell via a trigonal H-bond network with d(O(I)-O(II)) = 2.7-2.9 A, also in agreement with experimental measurements. The first shell tilt angle average is 33.4 degrees as compared to the reported value of 41 degrees . Wannier-Boys orbitals (WBO) show an interaction between the unoccupied 3d orbitals of the Fe(3+) valence (spin up, 3d(5)) and the occupied spin down lone pair orbitals of first shell waters. The effect of the spin ordering of the Fe(3+) ion on the WBO is not observed beyond the first shell. From this local bond analysis and consistent with other observations, the electronic structure of waters in the second shell is similar to that of a bulk water even in this strongly interacting system. H-bond decomposition shows significant bulk-like structure within the second shell for Fe(3+). The vibrational density of states shows a first shell red shift of 230 cm(-1) for the v(1),2v(2),v(3) overtone, in reasonable agreement with experimental estimates for trivalent cations (300 cm(-1)). No exchanges between first and second shell were observed. Waters in the second shell exchanged with bulk waters via dissociative and associative mechanisms. Results are compared with an AIMD study of Al(3+) and 64 waters. For Fe(3+) the average first shell tilt angle is larger and the tilt angle distribution wider. H-bond decomposition shows that second shell to second shell H-bonding is enhanced in Fe(3+) suggesting an earlier onset of bulk-like water structure.

A dianionic phosphorane intermediate and transition states in an associative A(N)+D(N) mechanism for the ribonucleaseA hydrolysis reaction.

The RNaseA enzyme efficiently cleaves phosphodiester bonds in the RNA backbone. Phosphoryl transfer plays a central role in many biochemical reactions, and this is one of the most studied enzymes. However, there remains considerable controversy about the reaction mechanism. Most of this debate centers around the roles of the conserved residues, structures of the transition state or states, the possibility of a stable intermediate, and the charge and structure of this intermediate. In this communication we report calculations of the mechanism of the hydrolysis step in this reaction using a comprehensive QM/MM theoretical approach that includes a high level calculation of the interactions in the QM region, free energy estimates along an NEB optimized reaction path, and the inclusion of the interaction of the protein surroundings and solvent. Contrary to prior calculations we find a stable pentacoordinated dianionic phosphorane intermediate in the reaction path supporting an A(N)+D(N) reaction mechanism. In the transition state in the path from the reactant to the intermediate state (with barrier of 3.96 kcal/mol and intermediate stability of 2.21 kcal/mol) a proton from the attacking water is partially transferred to the His119 residue and the PO bond only partially formed from the remaining nucleophilic OH(-) species (bond order (BO) 0.11). In passing from the intermediate to the product state (barrier 13.22 kcal/mol) the PO bond on the cyclic phosphorane intermediate is nearly broken (BO 0.28) and the transfer of the proton from the Lys41 is almost complete (Lys41-H BO 0.87). In the product state a proton has been transferred from Lys41 to the O2' position of the sugar. The role of Lys41 as the catalytic acid is a result of the relative positioning of the Lys41 and His12 in the catalytic site. This configuration is supported by calculations and docking studies.

Phosphorylation reaction in cAPK protein kinase-free energy quantum mechanical/molecular mechanics simulations.

We present results of a theoretical analysis of the phosphorylation reaction in cAMP-dependent protein kinase using a combined quantum mechanical and molecular mechanics (QM/MM) approach. Detailed analysis of the reaction pathway is provided using a novel QM/MM implementation of the nudged elastic band method, finite temperature fluctuations of the protein environment are taken into account using free energy calculations, and an analysis of hydrogen bond interactions is performed on the basis of calculated frequency shifts. The late transfer of the substrate proton to the conserved aspartate (D166), the activation free energy of 15 kcal/mol, and the slight exothermic (-3 kcal/mol) character of the reaction are all consistent with the experimental data. The near attack conformation of D166 in the reactant state is maintained by interactions with threonine-201, asparagine-177, and most notably by a conserved water molecule serving as a strong structural link between the primary metal ion and the D166. The secondary Mg ion acts as a Lewis acid, attacking the beta-gamma bridging oxygen of ATP. This interaction, along with a strong hydrogen bond between the D166 and the substrate, contributes to the stabilization of the transition state. Lys-168 maintains a hydrogen bond to a transferring phosphoryl group throughout a reaction process. This interaction increases in the product state and contributes to its stabilization.

Structure and dynamics of the hydration shells of the Al3+ ion.

First principles simulations of the hydration shells surrounding Al3+ ions are reported for temperatures near 300 degrees C. The predicted six water molecules in the octahedral first hydration shell were found to be trigonally coordinated via hydrogen bonds to 12 s shell water molecules in agreement with the putative structure used to analyze the x-ray data, but in disagreement with the results reported from conventional molecular dynamics using two-and three-body potentials. Bond lengths and angles of the water molecules in the first and second hydration shells and the average radii of these shells also agreed very well with the results of the x-ray analysis. Water transfers into and out of the second solvation shell were observed to occur on a picosecond time scale via a dissociative mechanism. Beyond the second shell the bonding pattern substantially returned to the tetrahedral structure of bulk water. Most of the simulations were done with 64 solvating water molecules (20 ps). Limited simulations with 128 water molecules (7 ps) were also carried out. Results agreed as to the general structure of the solvation region and were essentially the same for the first and second shell. However, there were differences in hydrogen bonding and Al-O radial distribution function in the region just beyond the second shell. At the end of the second shell a nearly zero minimum in the Al-O radial distribution was found for the 128 water system. This minimum is less pronounced minimum found for the 64 water system, which may indicate that sizes larger than 64 may be required to reliably predict behavior in this region.

The role of the putative catalytic base in the phosphoryl transfer reaction in a protein kinase: first-principles calculations.

Protein kinases are important enzymes controlling the majority of cellular signaling events via a transfer of the gamma-phosphate of ATP to a target protein. Even after many years of study, the mechanism of this reaction is still poorly understood. Among many factors that may be responsible for the 1011-fold rate enhancement due to this enzyme, the role of the conserved aspartate (Asp166) has been given special consideration. While the essential presence of Asp166 has been established by mutational studies, its function is still debated. The general base catalyst role assigned to Asp166 on the basis of its position in the active site has been brought into question by the pH dependence of the reaction rate, isotope measurements, and pre-steady-state kinetics. Recent semiempirical calculations have added to the controversy surrounding the role of Asp166 in the catalytic mechanism. No major role for Asp166 has been found in these calculations, which have predicted the reaction process consisting of an early transfer of a substrate proton onto the phosphate group. These conclusions were inconsistent with experimental observations. To address these differences between experimental results and theory with a more reliable computational approach and to provide a theoretical platform for understanding catalysis in this important enzyme family, we have carried out first-principles structural and dynamical calculations of the reaction process in cAPK kinase. To preserve the essential features of the reaction, representations of all of the key conserved residues (82 atoms) were included in the calculation. The structural calculations were performed using the local basis density functional (DFT) approach with both hybrid B3LYP and PBE96 generalized gradient approximations. This kind of calculation has been shown to yield highly accurate structural information for a large number of systems. The optimized reactant state structure is in good agreement with X-ray data. In contrast to semiempirical methods, the lowest energy product state places the substrate proton on Asp166. First-principles molecular dynamics simulations provide additional support for the stability of this product state. The latter also demonstrate that the proton transfer to Asp166 occurs at a point in the reaction where bond cleavage at the PO bridging position is already advanced. This mechanism is further supported by the calculated structure of the transition state in which the substrate hydroxyl group is largely intact. A metaphoshate-like structure is present in the transition state, which is consistent with the X-ray structures of transition state mimics. On the basis of the calculated structure of the transition state, it is estimated to be 85% dissociative. Our analysis also indicates an increase in the hydrogen bond strength between Asp166 and substrate hydroxyl and a small decrease in the bond strength of the latter in the transition state. In summary, our calculations demonstrate the importance of Asp166 in the enzymatic mechanism as a proton acceptor. However, the proton abstraction from the substrate occurs late in the reaction process. Thus, in the catalytic mechanism of cAPK protein kinase, Asp166 plays a role of a "proton trap" that locks the transferred phosphoryl group to the substrate. These results resolve prior inconsistencies between theory and experiment and bring new understanding of the role of Asp166 in the protein kinase catalytic mechanism.