Polarizable Molecular Mechanics Force Field Research Articles

On the Accuracy of QM/MM Models: A Systematic Study of Intramolecular Proton Transfer Reactions of Amino Acids in Water.

This work presents a systematic assessment of QM/QM' and QM/MM models with respect to direct QM calculations for the tautomerization (neutral to zwitterion) reactions of amino acids (glycine, alanine, valine, aspartate, and neutral and protonated histidine) solvated in a 160 water cluster. The effect of varying QM region size and choice of embedding potentials, including fixed-charge and polarizable molecular mechanics force fields (TIP3P and EFP) and various semiempirical QM methods (PM7, GFN2-xTB, DFTBA, DFTB3, HF-3c, and PBEh-3c), on the accuracy of the models was examined. A surprising finding was that molecular mechanics force fields outperformed many of the semiempirical methods. Generally, the errors in the QM/QM' and QM/MM models converge slowly with respect to the QM region size, requiring 50 or more waters to be included in the QM region before the error in the model falls below 1 kcal mol-1 of its pure QM result. Different QM region selection schemes were also compared, and it was found that selection based on Natural Population Analysis (NPA) atomic charges significantly reduced the error in the QM/QM' and QM/MM models particularly if a low-quality embedding potential was used. It is envisaged that these results will be useful for the development of future hybrid QM models.

Exciton properties and optical spectra of light harvesting complex II from a fully atomistic description.

We present a fully atomistic simulation of linear optical spectra (absorption, fluorescence and circular dichroism) of the Light Harvesting Complex II (LHCII) trimer using a hybrid approach, which couples a quantum chemical description of the chlorophylls with a classical model for the protein and the external environment (membrane and water). The classical model uses a polarizable Molecular Mechanics force field, thus allowing mutual polarization effects in the calculations of the excitonic properties. The investigation is performed both on the crystal structure and on structures generated by a μs long classical molecular dynamics simulation of the complex within a solvated membrane. The results show that this integrated approach not only provides a good description of the excitonic properties and optical spectra without the need for additional refinements of the excitonic parameters, but it also allows an atomistic investigation of the relative importance of electronic, structural and environment effects in determining the optical spectra.

Coupled Cluster Theory with Induced Dipole Polarizable Embedding for Ground and Excited States.

In this work, we present the theory and implementation of the coupled cluster single and double excitations (CCSD) method combined with a classical polarizable molecular mechanics force field (MMPol) based on the induced dipole model. The method is developed to compute electronic excitation energies within the state specific (SS) and linear response (LR) formalisms for the interaction of the quantum mechanical and classical regions. Furthermore, we consider an approximate expression of the correlation energy, originally developed for CCSD with implicit solvation models, where the interaction term is linear in the coupled cluster density. This approximation allows us to include the explicit contribution of the environment to the CC equations without increasing the computational effort. The test calculations on microsolvated systems, where the CCSD/MMPol method is compared to full CCSD calculations, demonstrates the reliability of this computational protocol for all interaction schemes (errors < 2%). We also show that it is important to include induced dipoles on all atom centers of the classical region and that too diffuse functions in the basis set may be problematic due to too strong interaction with the environment.

Retardation in electron dynamics simulations based on time-dependent density functional theory

When a molecular system is subjected to an external electric perturbation, originating from an electromagnetic wave or from a charged moving particle (in molecule-ion collisions), the electrons are the first particles to respond, on the attosecond timescale. Electron dynamics (ED) can be simulated by so-called real-time time-dependent density functional theory (RT-TDDFT). Within this framework, ED is driven by the electrostatic potential and by exchange-correlation potentials that fluctuate on the attosecond time scale. In vacuum the speed of light approaches 3 A as−1. Therefore, when simulating ED in extended molecular systems the question of retardation in the propagation of the potentials has to be posed. In this contribution we investigate two types of retardation; the first one deals with retardation in the potential created by a collision with a charged projectile. This is done through the Lienard-Wiechert potential (LWP). The second one deals with retardation in the electrostatic interaction between the time-dependent electron density and its environment, here in the context of hybrid schemes coupling RT-TDDFT to polarizable Molecular Mechanics force fields (MMpol). We found that the latter retardation effects can be safely neglected because of the rapid damping vs. distance of the electric fields created by electrostatic dipole moments. This conclusion is also relevant for methodologies, coupling RT-TDDFT to implicit polarizable continuum models. On the other hand, our results recommend the use of the LWP for modelling molecule-ion collisions by first-principles simulations. Remarkably, ionization takes place on faster time scales when relativistic corrections are introduced even for incident kinetic energies of 0.1 MeV.

Polarizable Molecular Simulations Reveal How Silicon-Containing Functional Groups Govern the Desalination Mechanism in Nanoporous Graphene.

We report a molecular dynamics (MD) simulation study of reverse osmosis desalination using nanoporous monolayer graphene passivated by SiH2 and Si(OH)2 functional groups. A highly accurate and detailed polarizable molecular mechanics force field model was developed for simulating graphene nanopores of various sizes and geometries. The simulated water fluxes and ion rejection percentages are explained using detailed atomistic mechanisms derived from analysis of the simulation trajectories. Our main findings are as follows: (1) The Si(OH)2 pores possess superior ion rejection rates due to selective electrostatic repulsion of Cl- ions, but Na+ ions are attracted to the pore and block water transfer. (2) By contrast, the SiH2 pores operate via a steric mechanism that excludes ions based on the size and flexibility of their hydration layers. (3) In the absence of ions, water flux is directly proportional to the solvent accessible area within the pore; however, simulated fluxes are lower than those inferred from recent experimental work. We also provide some hypotheses that could resolve the differences between simulation and experiment.

Efficient treatment of induced dipoles.

Most existing treatments of induced dipoles in polarizable molecular mechanics force field calculations use either the self-consistent variational method, which is solved iteratively, or the "direct" approximation that is non-iterative as a result of neglecting coupling between induced dipoles. The variational method is usually implemented using assumptions that are only strictly valid under tight convergence of the induced dipoles, which can be computationally demanding to enforce. In this work, we discuss the nature of the errors that result from insufficient convergence and suggest a strategy that avoids such problems. Using perturbation theory to reintroduce the mutual coupling into the direct algorithm, we present a computationally efficient method that combines the precision of the direct approach with the accuracy of the variational approach. By analyzing the convergence of this perturbation series, we derive a simple extrapolation formula that delivers a very accurate approximation to the infinite order solution at the cost of only a few iterations. We refer to the new method as extrapolated perturbation theory. Finally, we draw connections to our previously published permanent multipole algorithm to develop an efficient implementation of the electric field and Thole terms and also derive some necessary, but not sufficient, criteria that force field parameters must obey.

Three-Dimensional RISM Integral Equation Theory for Polarizable Solute Models.

Modeling solute polarizability is a key ingredient for improving the description of solvation phenomena. In recent years, polarizable molecular mechanics force fields have emerged that circumvent the limitations of classical fixed charge force fields by the ability to adapt their electrostatic potential distribution to a polarizing environment. Solvation phenomena are characterized by the solute's excess chemical potential, which can be computed by expensive fully atomistic free energy simulations. The alternative is to employ an implicit solvent model, which poses a challenge to the formulation of the solute-solvent interaction term within a polarizable framework. Here, we adapt the three-dimensional reference interaction site model (3D RISM) integral equation theory as a solvent model, which analytically yields the chemical potential, to the polarizable AMOEBA force field using an embedding cluster (EC-RISM) strategy. The methodology is analogous to our earlier approach to the coupling of a quantum-chemical solute description with a classical 3D RISM solvent. We describe the conceptual physical and algorithmic basis as well as the performance for several benchmark cases as a proof of principle. The results consistently show reasonable agreement between AMOEBA and quantum-chemical free energies in solution in general and allow for separate assessment of energetic and solvation-related contributions. We find that, depending on the parametrization, AMOEBA reproduces the chemical potential in better agreement with reference quantum-chemical calculations than the intramolecular energies, which suggests possible routes toward systematic improvement of polarizable force fields.

Ab Initio Extension of the AMOEBA Polarizable Force Field to Fe(2.).

We extend the AMOEBA polarizable molecular mechanics force field to the Fe(2+) cation in its singlet, triplet, and quintet spin states. Required parameters are obtained either directly from first principles calculations or optimized so as to reproduce corresponding interaction energy components in a hexaaquo environment derived from quantum mechanical energy decomposition analyses. We assess the importance of the damping of point-dipole polarization at short distance as well as the influence of charge-transfer for metal-water interactions in hydrated Fe(2+); this analysis informs the selection of model systems employed for parametrization. We validate our final Fe(2+) model through comparison of molecular dynamics (MD) simulations to available experimental data for aqueous ferrous ion in its quintet electronic ground state.

Computation of Induced Dipoles in Molecular Mechanics Simulations Using Graphics Processors

In this work, we present a tentative step toward the efficient implementation of polarizable molecular mechanics force fields with GPU acceleration. The computational bottleneck of such applications is found in the treatment of electrostatics, where higher-order multipoles and a self-consistent treatment of polarization effects are needed. We have implemented a GPU accelerated code, based on the Tinker program suite, for the computation of induced dipoles. The largest test system used shows a speedup factor of over 20 for a single precision GPU implementation, when comparing to the serial CPU version. A discussion of the optimization and parametrization steps is included. Comparison between different graphic cards and CPU-GPU embedding is also given. The current work demonstrates the potential usefulness of GPU programming in accelerating this field of applications.

Computational Studies of the Molecular Mechanisms Responsible for Ca2+ Permeation and Mg2+ Block of NMDA Receptors

Glutamate receptors are membrane proteins activated by the neurotransmitter glutamate that mediate fast synaptic excitation in the mammalian brain. NMDA receptors constitute a glutamate receptor subfamily specifically activated by N-methyl-d-aspartate. Due to their high Ca2+ permeability and voltage-dependent channel block by Mg2+, NMDA receptors play a central role in development through stabilization of synaptic connections, as well as in learning and memory by mediating many forms of synaptic plasticity. The mechanisms by which the ion channel of NMDA receptors selects Ca2+ for permeation over all other physiological ions, while binding Mg2+ and restricting its permeation, are not well understood. We hypothesize that the slightly different radii and electronic properties of Mg2+ and Ca2+ ions result in drastically different free energy barriers for transition of the ions from a binding site in the selectivity filter toward the intracellular solution. We are applying quantitative theoretical “bottom up” approaches to this complex system by combining methods of computational chemistry, molecular mechanics (MM), and bioinformatics. The structure of the NMDA receptor channel is constructed and refined using experimental information, homology modeling and extensive molecular dynamics simulations. We are performing quantum chemical calculations to determine the energy of the transition state of model ligand exchange reactions that mimic divalent ion transition from the selectivity filter of NMDA receptors to water. Quantum calculations are used to parameterize a polarizable molecular mechanics force field for these divalent ion interactions with organic ligands with the further goal to perform MM simulations. Umbrella sampling and thermodynamic integration simulations are used to compute free energies of transfer of divalent ions between water and a model of the NMDA ion selectivity filter as well as free energy barriers for these transitions.

Molecular Mechanisms of Ca 2+ Selectivity and Mg 2+ Block of NMDA Receptors

Glutamate receptors are membrane proteins activated by the neurotransmitter glutamate that mediate synaptic excitation in the mammalian brain. NMDA receptors constitute a glutamate receptor subfamily specifically activated by N-methyl-d-aspartate. Due to their high Ca2+ permeability and voltage-dependent channel block by Mg2+, NMDA receptors play a central role in development through stabilization of synaptic connections, as well as in learning and memory by mediating many forms of synaptic plasticity. The mechanisms by which the ion channel of NMDA receptors selects Ca2+ for permeation over all other physiological ions, while binding and restricting Mg2+ permeation, are not well understood. We hypothesize that the slightly different radii of Mg2+ and Ca2+ ions result in drastically different free energy barriers for movement of the ions from a binding site in the selectivity filter to water. We applied quantitative theoretical “bottom up” approaches to this complex system by combining methods of computational chemistry, molecular mechanics, and bioinformatics. We performed high accuracy quantum chemical (QC) calculations to determine the energy of the transition state of the ligand exchange reaction that occurs when divalent ions transition from the selectivity filter of NMDA receptors to water. The results of QC calculations are used to parameterize accurately the polarizable molecular mechanics force fields for divalent ion-ligand interactions. Molecular dynamics simulations with polarizable and non-polarizable force fields are used to predict quantitatively ion interactions with the heterogeneous environment of the pore. The structure of the NMDA receptor channel is constructed and refined using a combination of experimental studies, homology modeling, molecular dynamics, and an optimized structure of a Mg2+ - Asn ligand complex. A homology model of NMDA receptors is refined through extensive equilibration and all-atom molecular dynamics simulations on the multi-nanosecond time-scale to acquire a stable model structure.

Spectroscopic Polarizable Force Field for Amide Groups in Polypeptides

The infrared spectra of polypeptides are dominated by the so-called amide bands. These bands originate from the electrostatically coupled vibrations of the strongly polar amide groups (AGs) making up the polypeptide backbone. Because the AGs are highly polarizable, external electric fields can shift the frequencies of the amide normal modes over wide spectral ranges. The sensitivity to external fields and the strong polarity are the reasons why the shapes of the amide bands can code the structure of the polypeptide backbone. Aiming at a decoding of these band shapes, Schultheis et al. (J. Phys. Chem. B 2008, 112, 12217) have recently suggested a polarizable molecular mechanics (PMM) force field for AGs, which employs field dependent force constants and enables the computation of the amide bands from molecular dynamics simulations. Here we extend and refine this first suggestion of such a PMM force field. The extension rests on the choice of suitable internal coordinates for the AGs and on the inclusion of the complete AG Hessian and of its field dependence. The force field parameters are calculated from density functional theory. The improved quality of the resulting PMM descriptions is demonstrated using very simple examples and an outlook is given.

A Polarizable Force Field for Computing the Infrared Spectra of the Polypeptide Backbone

The shapes of the amide bands in the infrared (IR) spectra of proteins and peptides are caused by electrostatically coupled vibrations within the polypeptide backbone and code the structures of these biopolymers. A structural decoding of the amide bands has to resort to simplified models because the huge size of these macromolecules prevents the application of accurate quantum mechanical methods such as density functional theory (DFT). Previous models employed transition-dipole coupling methods that are of limited accuracy. Here we propose a concept for the computation of protein IR spectra, which describes the molecular mechanics (MM) of polypeptide backbones by a polarizable force field of "type II". By extending the concepts of conventional polarizable MM force fields, such a PMM/II approach employs field-dependent parameters not only for the electrostatic signatures of the molecular components but also for the local potentials modeling the stiffness of chemical bonds with respect to elongations, angle deformations, and torsions. Using a PMM/II force field, the IR spectra of the polypeptide backbone can be efficiently calculated from the time dependence of the backbone's dipole moment during a short (e.g., 100 ps) MD simulation by Fourier transformation. PMM/II parameters are derived for harmonic bonding potentials of amide groups in polypeptides from a series of DFT calculations on the model molecule N-methylacetamide (NMA) exposed to homogeneous external electric fields. The amide force constants are shown to vary by as much as 20% for relevant field strengths. As a proof of principle, it is shown that the large solvatochromic effects observed in the IR spectra of NMA upon transfer from the gas phase into aqueous solution are not only excellently reproduced by DFT/MM simulations but are also nicely modeled by the PMM/II approach. The tasks remaining for a proof of practice are specified.

Development of an Accurate and Robust Polarizable Molecular Mechanics Force Field from ab Initio Quantum Chemistry

We have produced a polarizable force field for a series of small molecules, representative of functional groups found in organic and biochemical systems. We have used high-level ab initio results for fitting values of all the parameters except for the dispersion-term coefficient B in the −B/r6 energy term, which, although obtained from comparison with experimental condensed-phase data, depended only on atomic number of the site in hand. Heats of vaporization and densities of the pure liquids, computed with molecular dynamics, agreed with experiment within ca. 0.5 kcal/mol and 5%, respectively.