Ab Initio Extension of the AMOEBA Polarizable Force Field to Fe2+

Solvation energies of the ferrous ion in water and in ammonia at various temperatures

CONTEXT

The solvation of metal ions is crucial to understanding relevant properties in physics, chemistry, or biology. Therefore, we present solvation enthalpies and solvation free energies of the ferrous ion in water and ammonia. Our results agree well with the experimental reports for the hydration free energy and hydration enthalpy. We obtained [Formula: see text] kJ mol[Formula: see text] for the hydration free energy and [Formula: see text] kJ mol[Formula: see text] for the hydration enthalpy of ferrous ion in water at room temperature. At ambient temperature, we obtained [Formula: see text] kJ mol[Formula: see text] as the [Formula: see text] ammoniation free energy and [Formula: see text] kJ mol[Formula: see text] for the ammoniation enthalpy. In addition, the free energy of solvation is deeply affected when the temperature increases. This pattern can be attributed to the rise of entropy when the temperature rises. Besides, the temperature does not affect the ammoniation enthalpies and the hydration enthalpy of the [Formula: see text] ion.

METHOD

All the geometry optimizations are performed at the MP2 methods associated with the 6-31++g(d,p) basis set of Pople. solvated phase structures of [Formula: see text] ion in water or in ammonia are performed using the PCM model. The [Formula: see text] program suite was used to perform all the calculations. The program TEMPO was also used to evaluate the temperature sensitivity of the different obtained geometries.

Thermodynamics of Metal-Acetate Interactions

Metal ions play crucial roles in protein- and ligand-mediated interactions. They not only act as catalysts to facilitate biological processes but are also important as protein structural elements. Accurately predicting metal ion interactions in computational studies has always been a challenge, and various methods have been suggested to improve these interactions. One such method is the 12-6-4 Lennard-Jones (LJ)-type nonbonded model. Using this model, it has been possible to successfully reproduce the experimental properties of metal ions in aqueous solution. The model includes induced dipole interactions typically ignored in the standard 12-6 LJ nonbonded model. In this we expand the applicability of this model to metal ion-carboxylate interactions. Using 12-6-4 parameters that reproduce the solvation free energies of the metal ions leads to an overestimation of metal ion-acetate interactions, thus, prompting us to fine-tune the model to specifically handle the latter. We also show that the standard 12-6 LJ model significantly falls short in reproducing the experimental binding free energy between acetate and 11 metal ions (Ni(II), Mg(II), Cu(II), Zn(II), Co(II), Cu(I), Fe(II), Mn(II), Cd(II), Ca(II), and Ag(I)). In this study, we describe optimized C4 parameters for the 12-6-4 LJ nonbonded model to be used with three widely employed water models (Transferable Intermolecular Potential with 3 Points (TIP3P), Simple Point Charge Extended (SPC/E), and Optimal Point Charge (OPC) water models). These parameters can accurately match the experimental binding free energy between 11 metal ions and acetate. These parameters can be applied to the study of metalloproteins and transition metal ion channels and transporters, as acetate serves as a representative of the negatively charged amino acid side chains from aspartate and glutamate.

Hybrid Functional and Plane Waves based Ab Initio Molecular Dynamics Study of the Aqueous Fe2+ /Fe3+ Redox Reaction

Kohn-Sham density functional theory and plane wave basis set based ab initio molecular dynamics (AIMD) simulation is a powerful tool for studying complex reactions in solutions, such as electron transfer (ET) reactions involving Fe²⁺ /Fe³⁺ ions in water. In most cases, such simulations are performed using density functionals at the level of Generalized Gradient Approximation (GGA). The challenge in modelling ET reactions is the poor quality of GGA functionals in predicting properties of such open-shell systems due to the inevitable self-interaction error (SIE). While hybrid functionals can minimize SIE, standard plane-wave based AIMD at that level of theory is typically 150 times slower than GGA for systems containing ∼100 atoms. Among several approaches reported to speed-up AIMD simulations with hybrid functionals, the noise-stabilized MD (NSMD) procedure, together with the use of localized orbitals to compute the required exchange integrals, is an attractive option. In this work, we demonstrate the application of the NSMD approach for studying the Fe²⁺ /Fe³⁺ redox reaction in water. It is shown here that long AIMD trajectories at the level of hybrid density functionals can be obtained using this approach. Redox properties of the aqueous Fe²⁺ /Fe³⁺ system computed from these simulations are compared with the available experimental data for validation.

Manifolds of low energy structures for a magic number of hydrated sulfate: SO42-(H2O)24

Low energy structures of SO₄^2-(H₂O)₂₄ have been obtained using a combination of classical molecular dynamics simulations and refinement of structures and energies by quantum chemical calculations. Extensive exploration of the potential energy surface led to a number of low-energy structures, confirmed by accurate calibration calculations. An overall analysis of this large set was made after devising appropriate structural descriptors such as the numbers of cycles and their combinations. Low energy structures bear common motifs, the most prominent being fused cycles involving alternatively four and six water molecules. The latter adopt specific conformations which ensure the appropriate surface curvature to form a closed cage without dangling O-H bonds and at the same time provide 12-coordination of the sulfate ion. A prominent feature to take into account is isomerism via inversion of hydrogen bond orientations along cycles. This generates large families of ca. 100 isomers for this cluster size, spanning energy windows of 10-30 kJ mol^-1. This relatively ignored isomerism must be taken into account to identify reliably the lowest energy minima. The overall picture is that the magic number cluster SO₄^2-(H₂O)₂₄ does not correspond to formation of a single, remarkable structure, but rather to a manifold of structural families with similar stabilities. Extensive calculations on isomerization mechanisms within a family indicate that large barriers are associated to direct inversion of hydrogen bond networks. Possible implications of these results for magic number clusters of other anions are discussed.

Evaluation of the AMOEBA force field for simulating metal halide perovskites in the solid state and in solution

In this work, we compare the existing nonpolarizable force fields developed to study the solid or solution phases of hybrid organic-inorganic halide perovskites with the AMOEBA polarizable force field. The aim is to test whether more computationally expensive polarizable force fields like AMOEBA offer better transferability between solution and solid phases, with the ultimate goal being the study of crystal nucleation, growth, and other interfacial phenomena involving these ionic compounds. In the context of hybrid perovskites, AMOEBA force field parameters already exist for several elements in solution, and we decided to leave them unchanged and to only parameterize the missing ones (Pb²⁺ and CH₃NH₃ ⁺ ions) in order to maximize transferability and avoid overfitting to the specific examples studied here. Overall, we find that AMOEBA yields accurate hydration free energies (within 5%) for typical ionic species while showing the correct ordering of stability for the different crystal polymorphs of CsPbI₃ and CH₃NH₃PbI₃. Although the existing parameters do not accurately reproduce all transition temperatures and lattice parameters, AMOEBA offers better transferability between solution and solid states than existing nonpolarizable force fields.

Evidence for an egg-box-like structure in iron(ii)–polygalacturonate hydrogels: a combined EXAFS and molecular dynamics simulation study

The coordination of Fe(ii) with polygalacturonic acid (polyGalA) in Fe(ii)–polyGalA hydrogels exhibits an octahedral geometry that follows the “egg-box model”.

Polarizable Force Fields for Biomolecular Simulations: Recent Advances and Applications

Realistic modeling of biomolecular systems requires an accurate treatment of electrostatics, including electronic polarization. Due to recent advances in physical models, simulation algorithms, and computing hardware, biomolecular simulations with advanced force fields at biologically relevant timescales are becoming increasingly promising. These advancements have not only led to new biophysical insights but also afforded opportunities to advance our understanding of fundamental intermolecular forces. This article describes the recent advances and applications, as well as future directions, of polarizable force fields in biomolecular simulations.

Infrared Spectra of Deprotonated Dicarboxylic Acids: IRMPD Spectroscopy and Empirical Valence‐Bond Modeling

Experimental infrared multiple-photon dissociation (IRMPD) spectra recorded for a series of deprotonated dicarboxylic acids, HO₂ (CH₂ )_n CO 2 - (n=2-4), are interpreted using a variety of computational methods. The broad bands centered near 1600 cm^-1 can be reproduced neither by static vibrational calculations based on quantum chemistry nor by a dynamical description of individual structures using the many-body polarizable AMOEBA force field, strongly suggesting that these molecules experience dynamical proton sharing between the two carboxylic ends. To confirm this assumption, AMOEBA was combined with a two-state empirical valence-bond (EVB) model to allow for proton transfer in classical molecular dynamics simulations. Upon suitable parametrization based on ab initio reference data, the EVB-AMOEBA model satisfactorily reproduces the experimental infrared spectra, and the finite temperature dynamics reveals a significant amount of proton sharing in such systems.

Classical Drude Polarizable Force Field Model for Methyl Phosphate and Its Interactions with Mg2

Phosphate groups are essential components of nucleic acids and proteins, whose interactions with solvent, metal ions, and ionic side chains help control folding and binding. Methyl phosphate (MP) represents a simple analog of phosphate moieties that are post-translation modifications in proteins and present at the termini of nucleic acids, among other environments. In the present study, we optimized parameters for use in polarizable molecular dynamics simulations of MP in its mono- and dianionic forms, MP- ≡ CH3HPO4- and MP2- ≡ CH3PO42-, along with P i2- ≡ HPO42-, in the context of the classical Drude oscillator model. Parameter optimization was done in a manner consistent with the remainder of the Drude molecular mechanics force field, choosing atomic charges and polarizabilities to reproduce molecular properties from quantum mechanics as well as experimental hydration free energies. Optimized parameters were similar to existing dimethyl phosphate parameters, with a few significant differences. The developed parameters were then used to compute magnesium binding affinities in aqueous solution, using alchemical molecular dynamics free energy simulations. Good agreement with experiment was obtained, and outer sphere binding was shown to be predominant for MP- and MP2-.

Chemical shift extremum of 129Xe(aq) reveals details of hydrophobic solvation

The ¹²⁹Xe chemical shift in an aqueous solution exhibits a non-monotonic temperature dependence, featuring a maximum at 311 K. This is in contrast to most liquids, where the monotonic decrease of the shift follows that of liquid density. In particular, the shift maximum in water occurs at a higher temperature than that of the maximum density. We replicate this behaviour qualitatively via a molecular dynamics simulation and computing the ¹²⁹Xe chemical shift for snapshots of the simulation trajectory. We also construct a semianalytical model, in which the Xe atom occupies a cavity constituted by a spherical water shell, consisting of an even distribution of solvent molecules. The temperature dependence of the shift is seen to result from a product of the decreasing local water density and an increasing term corresponding to the energetics of the Xe-H₂O collisions. The latter moves the chemical shift maximum up in temperature, as compared to the density maximum. In water, the computed temperature of the shift maximum is found to be sensitive to both the details of the binary chemical shift function and the coordination number. This work suggests that, material parameters allowing, the maximum should be exhibited by other liquids, too.

AMOEBA Polarizable Force Field Parameters of the Heme Cofactor in Its Ferrous and Ferric Forms

We report the first parameters of the heme redox cofactors for the polarizable AMOEBA force field in both the ferric and ferrous forms. We consider two types of complexes, one with two histidine side chains as axial ligands and one with a histidine and a methionine side chain as ligands. We have derived permanent multipoles from second-order Møller-Plesset perturbation theory (MP2). The sets of parameters have been validated in a first step by comparison of AMOEBA interaction energies of heme and a collection of biologically relevant molecules with MP2 and Density Functional Theory (DFT) calculations. In a second validation step, we consider interaction energies with large aggregates comprising around 80 H₂O molecules. These calculations are repeated for 30 structures extracted from semiempirical PM7 DM simulations. Very encouraging agreement is found between DFT and the AMOEBA force field, which results from an accurate treatment of electrostatic interactions. We finally report long (10 ns) MD simulations of cytochromes in two redox states with AMOEBA testing both the 2003 and 2014 AMOEBA water models. These simulations have been carried out with the TINKER-HP (High Performance) program. In conclusion, owing to their ubiquity in biology, we think the present work opens a wide array of applications of the polarizable AMOEBA force field on hemeproteins.

Strategy for Modeling the Infrared Spectra of Ion-Containing Water Drops

Hydrated ions are ubiquitous in environmental and biological media. Understanding the perturbation exerted by an ion on the water hydrogen bond network is possible in the nanodrop regime by recording vibrational spectra in the O-H bond stretching region. This has been achieved experimentally in recent years by forming gaseous ions containing tens to hundreds of water molecules and recording their infrared photodissociation spectra. In this paper, we demonstrate the capabilities of a modeling strategy based on an extension of the AMOEBA polarizable force field to implement water atomic charge fluctuations along with those of intramolecular structure along the dynamics. This supplementary flexibility of nonbonded interactions improves the description of the hydrogen bond network and, therefore, the spectroscopic response. Finite temperature IR spectra are obtained from molecular dynamics simulations by computing the Fourier transform of the dipole moment autocorrelation function. Simulations of 1-2 ns are required for extensive sampling in order to reproduce the experimental spectra. Furthermore, bands are assigned with the driven molecular dynamics approach. This method package is shown to compare successfully with experimental spectra for 11 ions in water drops containing 36-100 water molecules. In particular, band frequency shifts of the free O-H stretching modes at the cluster surface are well reproduced as a function of both ion charge and drop size.

Extending the Nonbonded Cationic Dummy Model to Account for Ion-Induced Dipole Interactions

Modeling metalloproteins often requires classical molecular dynamics (MD) simulations in order to capture their relevant motions, which in turn necessitates reliable descriptions of the metal centers involved. One of the most successful approaches to date is provided by the "cationic dummy model", where the positive charge of the metal ion is transferred toward dummy particles that are bonded to the central metal ion in a predefined coordination geometry. While this approach allows for ligand exchange, and captures the correct electrostatics as demonstrated for different divalent metal ions, current dummy models neglect ion-induced dipole interactions. In the present work, we resolve this weakness by taking advantage of the recently introduced 12-6-4 type Lennard-Jones potential to include ion-induced dipole interactions. We revise our previous dummy model for Mg²⁺ and demonstrate that the resulting model can simultaneously reproduce the experimental solvation free energy and metal-ligand distances without the need for artificial restraints or bonds. As ion-induced dipole interactions become particularly important for highly charged metal ions, we develop dummy models for the biologically relevant ions Al³⁺, Fe³⁺, and Cr³⁺. Finally, the effectiveness of our new models is demonstrated in MD simulations of several diverse (and highly challenging to simulate) metalloproteins.

Computational investigation of O2 diffusion through an intra-molecular tunnel in AlkB; influence of polarization on O2 transport

E. Coli AlkB catalyzes the direct dealkylation of various alkylated bases in damaged DNA. The diffusion of molecular oxygen to the active site in AlkB is an essential step for the oxidative dealkylation activity. Despite detailed studies on the stepwise oxidation mechanism of AlkB, there is no conclusive picture of how O2 molecules reach the active site of the protein. Yu et al. (Nature, 439, 879) proposed the existence of an intra-molecular tunnel based on their initial crystal structures of AlkB. We have employed computational simulations to investigate possible migration pathways inside AlkB for O2 molecules. Extensive molecular dynamics (MD) simulations, including explicit ligand sampling and potential of mean force (PMF) calculations, have been performed to provide a microscopic description of the O2 delivery pathway in AlkB. Analysis of intra-molecular tunnels using the CAVER software indicates two possible pathways for O2 to diffuse into the AlkB active site. Explicit ligand sampling simulations suggests that only one of these tunnels provides a viable route. The free energy path for an oxygen molecule to travel along each of these tunnels has been determined with AMBER and AMOEBA. Both PMFs indicate passive transport of O2 from the surface of the protein. However, the inclusion of explicit polarization shows a very large barrier for diffusion of the co-substrate out of the active site, compared with the non-polarizable potential. In addition, our results suggest that the mutation of a conserved residue along the tunnel, Y178, has dramatic effects on the dynamics of AlkB and on the transport of O2 along the tunnel.

Polarizable and Non-Polarizable Force Field Representations of Ferric Cation and Validations

The AMOEBA polarizable force field of ferric ion was optimized and applied to study the hydration of ferric ion and its complexation with porphine in the aqueous phase. The nonpolarizable force field was also optimized for comparison. The AMOEBA force field was found to give a more accurate hydration free energy than the nonpolarizable force field with respect to experimental data, and correctly predict the most stable electronic state of hydrated Fe³⁺, which is the sextet state, and of the Fe^(III)-Por complex, which is the quartet state, consistent with the literature that was carried out using the DFT method. The explicit inclusion of charge transfer between Fe³⁺ and ligand was found to be important in order to obtain a precise picture of polarization energy and van der Waals energy, which otherwise deviate from the corresponding energy components derived from ab initio calculations. The successful application of the AMOEBA force field in the characterization of aquo Fe^(III)-Por complexes suggests that its use may be extended to the study of the dynamics of metalloenzyme containing highly charged metal ions in the condensed phase with reliable treatment of the interactions between metal atom and protein.

Metal Ion Modeling Using Classical Mechanics

Metal ions play significant roles in numerous fields including chemistry, geochemistry, biochemistry, and materials science. With computational tools increasingly becoming important in chemical research, methods have emerged to effectively face the challenge of modeling metal ions in the gas, aqueous, and solid phases. Herein, we review both quantum and classical modeling strategies for metal ion-containing systems that have been developed over the past few decades. This Review focuses on classical metal ion modeling based on unpolarized models (including the nonbonded, bonded, cationic dummy atom, and combined models), polarizable models (e.g., the fluctuating charge, Drude oscillator, and the induced dipole models), the angular overlap model, and valence bond-based models. Quantum mechanical studies of metal ion-containing systems at the semiempirical, ab initio, and density functional levels of theory are reviewed as well with a particular focus on how these methods inform classical modeling efforts. Finally, conclusions and future prospects and directions are offered that will further enhance the classical modeling of metal ion-containing systems.

Molecular dynamics study of response of liquid N,N-dimethylformamide to externally applied electric field using a polarizable force field

The behavior of Liquid N,N-dimethylformamide subjected to a wide range of externally applied electric fields (from 0.001 V/nm to 1 V/nm) has been investigated through molecular dynamics simulation. To approach the objective the AMOEBA polarizable force field was extended to include the interaction of the external electric field with atomic partial charges and the contribution to the atomic polarization. The simulation results were evaluated with quantum mechanical calculations. The results from the present force field for the liquid at normal conditions were compared with the experimental and molecular dynamics results with non-polarizable and other polarizable force fields. The uniform external electric fields of higher than 0.01 V/nm have a significant effect on the structure of the liquid, which exhibits a variation in numerous properties, including molecular polarization, local cluster structure, rotation, alignment, energetics, and bulk thermodynamic and structural properties.

Long-range electron transfer with myoglobin immobilized at Au/mixed-SAM junctions: mechanistic impact of the strong protein confinement

Horse muscle myoglobin (Mb) was tightly immobilized at Au-deposited ~15-Å-thick mixed-type (1:1) alkanethiol SAMs, HS-(CH₂)₁₁-COOH/HS-(CH₂)₁₁-OH, and placed in contact with buffered H₂O or D₂O solutions. Fast-scan cyclic voltammetry (CV) and a Marcus-equation-based analysis were applied to determine unimolecular standard rate constants and reorganization free energies for electron transfer (ET), under variable-temperature (15-55 °C) and -pressure (0.01-150 MPa) conditions. The CV signal was surprisingly stable and reproducible even after multiple temperature and pressure cycles. The data analysis revealed the following values: standard rate constant, 33 s⁻¹ (25 °C, 0.01 MPa, H₂O); reorganization free energy, 0.5 ± 0.1 eV (throughout); activation enthalpy, 12 ± 3 kJ mol⁻¹; activation volume, -3.1 ± 0.2 cm³ mol⁻¹; and pH-dependent solvent kinetic isotope effect (k(H)⁰/k(D)⁰), 0.7-1.4. Furthermore, the values for the rate constant and reorganization free energy are very similar to those previously found for cytochrome c electrostatically immobilized at the monocomponent Au/HS-(CH₂)₁₁-COOH junction. In vivo, Mb apparently forms a natural electrostatic complex with cytochrome b₅ (cyt-b₅) through the "dynamic" (loose) docking pattern, allowing for a slow ET that is intrinsically coupled to the water's removal from the "defective" heme iron (altogether shaping the biological repair mechanism for Mb's "met" form). In contrary, our experiments rather mimic the case of a "simple" (tight) docking of the redesigned (mutant) Mb with cyt-b₅ (Nocek et al. J. Am. Chem. Soc. 2010, 132, 6165-6175). According to our analysis, in this configuration, Mb's distal pocket (linked to the "ligand channel") seems to be arrested within the restricted configuration, allowing the rate-determining reversible ET process to be coupled only to the inner-sphere reorganization (minimal elongation/shortening of an Fe-OH₂ bond) rather than the pronounced detachment (rebinding) of water and, hence, to be much faster.