Force Field Simulations Research Articles

The mobility and solvation structure of a hydroxyl radical in a water nanodroplet: a Born-Oppenheimer molecular dynamics study.

Hydroxyl radicals (OH*) play a crucial role in atmospheric chemistry and biological processes. In this study, Born-Oppenheimer molecular dynamics simulations are performed under ambient conditions for a hydroxyl radical in a water nanodroplet containing 191 water molecules. Density functional theory calculations are performed at the BLYP-D3 level with some test calculations at the B3LYP-D3 level. In two 150 ps trajectories, either with OH* initially located in the interior region or at the surface of the water nanodroplet, the OH* radical ends up in the subsurface layer of the nanodroplet, which is different from the "surface preference" predicted from previous empirical force field simulations. The solvation structure of OH* contains fluctuating hydrogen bonds, as well as a two-center three-electron hemibond in some cases. The mobility of OH* is enhanced by hydrogen transfer, which has a free energy barrier of ∼4.6 kcal mol-1. The results presented in this study deepen our understanding of the structure and dynamics of OH* in aqueous solutions, especially around the air-water interface.

Effect of protein dimerization on ion conductivity of gramicidin a channel studied using polarizable force field

In studies of ion channel systems, due to the huge computational cost of polarizable force fields, classical force fields remain the most widely used for a long time. In this work, we used the AMOEBA polarizable atomic multipole force field in enhanced sampling simulations of single-channel gramicidin A (gA) and double-channel gA systems and investigated its reliability in characterizing ion-transport properties of the gA ion channel under dimerization. The influence of gA dimerization on the permeation of potassium and sodium ions through the channel was described in terms of conductance, diffusion coefficient, and free energy profile. Results from the polarizable force field simulations show that the conductance of potassium and sodium ions passing through the single- and double-channel agrees well with experimental values. Further data analysis reveals that the molecular mechanism of protein dimerization affects the ion-transport properties of gA channels, i.e., protein dimerization accelerates the permeation of potassium and sodium ions passing through the double-channel by adjusting the environment around gA protein (the distribution of phospholipid head groups, ions outside the channel, and bulk water), rather than directly adjusting the conformation of gA protein.

A Bayes-inspired theory for optimally building an efficient coarse-grained folding force field.

Because of their potential utility in predicting conformational changes and assessing folding dynamics, coarse-grained (CG) RNA folding models are appealing for rapid characterization of RNA molecules. Previously, we reported the iterative simulated RNA reference state (IsRNA) method for parameterizing a CG force field for RNA folding, which consecutively updates the simulation force field to reflect marginal distributions of folding coordinates in the structure database and extract various energy terms. While the IsRNA model was validated by showing close agreement between the IsRNA-simulated and experimentally observed distributions, here, we expand our theoretical understanding of the model and, in doing so, improve the parameterization process to optimize the subset of included folding coordinates, which leads to accelerated simulations. Using statistical mechanical theory, we analyze the underlying, Bayesian concept that drives parameterization of the energy function, providing a general method for developing predictive, knowledge-based, polymer force fields on the basis of limited data. Furthermore, we propose an optimal parameterization procedure, based on the principal of maximum entropy.

Lysine-arginine advanced glycation end-product cross-links and the effect on collagen structure: A molecular dynamics study.

The accumulation of advanced glycation end‐products is a fundamental process that is central to age‐related decline in musculoskeletal tissues and locomotor system function and other collagen‐rich tissues. However, although computational studies of advanced glycation end‐product cross‐links could be immensely valuable, this area remains largely unexplored given the limited availability of structural parameters for the derivation of force fields for Molecular Dynamics simulations. In this article, we present the bonded force constants, atomic partial charges and geometry of the arginine‐lysine cross‐links DOGDIC, GODIC, and MODIC. We have performed in vacuo Molecular Dynamics simulations to validate their implementation against quantum mechanical frequency calculations. A DOGDIC advanced glycation end‐product cross‐link was then inserted into a model collagen fibril to explore structural changes of collagen and dynamics in interstitial water. Unlike our previous studies of glucosepane, our findings suggest that intra‐collagen DOGDIC cross‐links furthers intra‐collagen peptide hydrogen‐bonding and does not promote the diffusion of water through the collagen triple helices.

Modification of lipid membrane compressibility induced by an electric field.

Changes in membrane deformation and compressibility, induced by an external electric field, are investigated using coarse-grained martini force field simulations in a salt-free environment. We observe changes in the area of the membrane above a critical electric field. Below this value, the membrane compressibility modulus is found to decrease monotonically. For higher electric fields, the membrane projected area remains constant while the net interfacial area increases, with the corresponding compressibility moduli, show the opposite behavior. We find that the mechanical parameters, surface tension and bending modulus, of a freely floating membrane in the absence of explicit ions, are unaffected by the presence of the electric field. We believe these results have a bearing on our understanding of the electroformation of uncharged lipids in a salt-free environment.

Conduction and Gating Properties of the TRAAK Channel from Molecular Dynamics Simulations with Different Force Fields.

In recent years, the K2P family of potassium channels has been the subject of intense research activity. Owing to the complex function and regulation of this family of ion channels, it is common practice to complement experimental findings with the atomistic description provided by computational approaches such as molecular dynamics (MD) simulations, especially, in light of the unprecedented timescales accessible at present. However, despite recent substantial improvements, the accuracy of MD simulations is still undermined by the intrinsic limitations of force fields. Here, we systematically assessed the performance of the most popular force fields employed to study ion channels at timescales that are orders of magnitude greater than the ones accessible when these energy functions were first developed. Using 32 μs of trajectories, we investigated the dynamics of a member of the K2P ion channel family, the TRAAK channel, using two established force fields in simulations of biological systems: AMBER and CHARMM. We found that while results are comparable on the nanosecond timescales, significant inconsistencies arise at microsecond timescales.

Computational and experimental studies on the influence of theTrp-cage inter-residue interactions and other external factors on its structural and folding stability

<p indent=0mm>The major driving force of protein folding depends on their inter-residue interactions, especially for sidechain-sidechain interactions. The Trp-cage is a fast-folding mini-protein, which has already been acknowledged as an ideal model for tackling complicated folding-related problems of proteins. Although this protein model is merely composed of 20 amino acid residues, it exhibits an extremely complex interaction network between residues with different types of sidechain orientation. Because current force fields exhibit different degrees of accuracy inadequacies in the descriptions of atomic interactions, it remains a matter of debate at what stage of folding the helix, hydrophobic cluster and salt-bridge is being conclusively formed in its possible folding mechanisms. The structures and mutations of all the Trp-cage system reported since 2002 are summarized in this paper, and its two typical structures (TC5b and TC10b) are taken as examples to summarize the composition and types of their inter-residue interactions as well as the sidechain features and orientations. Meanwhile, the recent computational and experimental studies of several possible folding mechanisms are also generalized and discussed, and their relationships with simulation temperature and force field as well as the switching from one folding pathway to the other are also discussed in detail. Finally, the effects of numerous factors on structural stability and folding dynamics have been elaborately discussed, including different types of inter-residue interactions (such as salt-bridge, hydrophobic and hydrogen-bonding interactions) and environmental factors (such as temperature, pressure, pH and other potential factors). We also propose that the clarification of the above factors and their synergistic effects are beneficial not only to reveal the driving force of protein stabilization and rapid formation, but also to provide theoretical guidance for further structural regulation and optimal design.

Millimeter-Scale and Billion-Atom Reactive Force Field Simulation on Sunway Taihulight

Large-scale molecular dynamics (MD) simulations on supercomputers play an increasingly important role in many research areas. With the capability of simulating charge equilibration (QEq), bonds and so on, Reactive force field (ReaxFF) enables the precise simulation of chemical reactions. Compared to the first principle molecular dynamics (FPMD), ReaxFF has far lower requirements on computational resources so that it can achieve higher efficiencies for large-scale simulations. In this article, we present our efforts on scaling ReaxFF on the Sunway TaihuLight Supercomputer (TaihuLight). We have carefully redesigned the force analysis and neighbor list building steps. By applying fine-grained optimizations we gain better single process performance. For the many-body interactions, we propose an isolated computation and update strategy and implement inverse trigonometric functions. For QEq, we implement a pipelined conjugate gradient (CG) approach to achieving better scalability. Furthermore, we reorganize the data layout and implement the update operation based on data locality in ReaxFF. Our experiments show that this approach can simulate chemical reactions with 1,358,954,496 atoms using 4,259,840 cores with a performance of 0.015 ns/day. To our best knowledge, this is the first realization of chemical reaction simulation with a millimeter-scale force field.

Different Force Fields Give Rise to Different Amyloid Aggregation Pathways in Molecular Dynamics Simulations.

The progress toward understanding the molecular basis of Alzheimers's disease is strongly connected to elucidating the early aggregation events of the amyloid-β (Aβ) peptide. Molecular dynamics (MD) simulations provide a viable technique to study the aggregation of Aβ into oligomers with high spatial and temporal resolution. However, the results of an MD simulation can only be as good as the underlying force field. A recent study by our group showed that none of the common force fields can distinguish between aggregation-prone and nonaggregating peptide sequences, producing a similar and in most cases too fast aggregation kinetics for all peptides. Since then, new force fields specially designed for intrinsically disordered proteins such as Aβ were developed. Here, we assess the applicability of these new force fields to studying peptide aggregation using the Aβ16-22 peptide and mutations of it as test case. We investigate their performance in modeling the monomeric state, the aggregation into oligomers, and the stability of the aggregation end product, i.e., the fibrillar state. A main finding is that changing the force field has a stronger effect on the simulated aggregation pathway than changing the peptide sequence. Also the new force fields are not able to reproduce the experimental aggregation propensity order of the peptides. Dissecting the various energy contributions shows that AMBER99SB-disp overestimates the interactions between the peptides and water, thereby inhibiting peptide aggregation. More promising results are obtained with CHARMM36m and especially its version with increased protein-water interactions. It is thus recommended to use this force field for peptide aggregation simulations and base future reparameterizations on it.

Multipolar Description of Atom-Atom Electrostatic Interaction Energies in Single/Double-Stranded DNAs.

Molecular force field simulation is an effective method to explore the properties of DNA molecules in depth. Almost all current popular force fields calculate atom-atom electrostatic interaction energies for DNAs based on the atomic charge and dipole or quadrupole moments, without considering high-rank atomic multipole moments for more accurate electrostatics. Actually, the distribution of electrons around atomic nuclei is not spherically symmetric but is geometry dependent. In this work, a multipole expansion method that allows us to combine polarizability and anisotropy was applied. One single-stranded DNA and one double-stranded DNA were selected as pilot systems. Deoxynucleotides were cut out from pilot systems and capped by mimicking the original DNA environment. Atomic multipole moments were integrated instead of fixed-point charges to calculate atom-atom electrostatic energies to improve the accuracy of force fields for DNA simulations. Also, the applicability of modeling the behavior of both single-stranded and double-stranded DNAs was investigated. The calculation results indicated that the models can be transferred from pilot systems to test systems, which is of great significance for the development of future DNA force fields.

Polarizabilities of neutral atoms and atomic ions with a noble gas electron configuration.

Atomic polarizabilities play an important role in the development of force fields for molecular simulations, as well as for the development of qualitative concepts of atomic and molecular behavior. Coupled cluster theory at the coupled cluster singles doubles triples level with very large correlation-consistent basis sets with extended diffuse functions has been used to predict the polarizabilities of the atomic neutrals, mono-cations and mono-anions with a noble gas configuration. Additional corrections for scalar relativistic and spin-orbit effects were also included for the electron configurations of Kr, Xe, and Rn. The results are in excellent agreement with experiment or with other high level calculations where available. The current results for most of these species represent the best available values for the polarizabilities. The results show that the polarizability of H- is very difficult to calculate without extremely diffuse functions. The polarizability of H- is the largest value, 34.05 Å3, calculated for all species in the current study. The polarizabilities of the remaining halogen anions are also the best available values. The polarizabilities of the halogen anions (excluding F-) and H- have a linear correlation with the electron affinity of the neutral atom. Spin-orbit effects, even for closed shell species, cannot be ignored for quantitative accuracy, and the inclusion of spin-orbit effects for Fr+, Rn, and At- increases the polarizability by 4%, 6%, and 15%, respectively.

Critical Assessment of Common Force Fields for Molecular Dynamics Simulations of Potassium Channels.

For the last two decades, the KcsA K+ channel has served as a case study to understand how potassium channels operate at the atomic scale, and molecular dynamics simulations have contributed significantly to the current knowledge of the atomic mechanisms of conduction, inactivation, and gating in this family of membrane proteins. Currently, microsecond trajectories are becoming the new standard in the field, and consequently, it is critical to assess and compare the performance of the classical force fields ordinarily used in simulations of biological systems as well as to quantitatively assess the agreement with experimental data for trajectories of this order of magnitude. To that extent, we performed classical molecular dynamics simulations with CHARMM36 and AMBER-ff14sb force fields using atomic models based on the experimental structure of the KcsA channel in the open/conductive state, at conditions of ionic concentrations and membrane potentials resembling the ones adopted in experiments. In simulations using the CHARMM force field, the experimental open/conductive structure of the channel exhibited high conformational plasticity and fast collapse toward an occluded state. In contrast, in an identical set of simulations using the AMBER force field, no major deviations from the experimental structure were recorded. Force field development is a complex process in which many approximations are typically required and adopted. The results presented here provide additional motivation to further improve the existing models and, crucially, alert practitioners about limitations.

Simplified force field for molecular dynamics simulations of amorphous SiO2 for solar applications

Mixtures of molten-salts (e.g. NaNO3, KNO3) with α-SiO2 nanoparticles have potential applications in concentrated solar power generation both as a thermal energy storage material and heat transfer fluid. In this study, a set of Lennard-Jones parameters for amorphous silica is proposed, aimed at understanding the molten-salt/amorphous-silica interface at the nanoscale through molecular dynamics simulations. Current force fields used for the simulation of α-SiO2 do not have validated rules for cross-term interactions or are computationally expensive. In this work, a novel set of Lennard-Jones parameters capable of accurately simulating the structure of amorphous SiO2 is presented. Its advantage is that it has precise validated mixing rules, which can be used to model cross term (different species) interactions, which are dominant in simulations involving interfaces. The proposed potential is validated by experimental results and data from previous studies: all the calculated properties (i.e. radial distribution, angular distribution, mean bond length, coordination number, density, thermal expansion coefficient and thermal conductivity) are in good agreement with the data in the literature. Additionally, the molten salt-amorphous silica interface is investigated. The existence of a compressed liquid layer at the liquid-nanoparticle interface is observed. The density of the molten salt mixture is found to be significantly higher and the self-diffusivity lower, in the proximity of the interface. Finally, the Na+ and K+ ions appear to be ordered in the proximity of the interface, indicating the existence of an ordered double-layer.

Force-Field Simulations of a Hydrated Lanreotide-Based Derivative: Hydration, Dynamics, and Numerical Evidence of Self-Assembly in Dimers.

Recently, self-organization of the cyclic octapeptide lanreotide and lanreotide-based derivatives in a nanotube to from a dimer structure has been experimentally evidenced. While the nature of the interactions between both monomers has been strongly investigated no molecular details of the hydration of the monomer and the formation of the dimer have been provided. Using molecular dynamics simulations, this work focuses on the structure, hydration, and dynamics of water and an analog of lanreotide. To do so, several models of monomers based on different schemes of partial charges and electrostatic interaction calculations are considered. By comparison with the experiments, we show that the model based on the combination of the AMBER force-field, CHELPG charge calculation, Ewald sum is the most relevant. Additionally, by mapping the interfacial hydration of the lanreotide monomer we evidence a heterogeneous surface in terms of hydrophilicity involving heterogeneous hydration. Furthermore, we show a slowdown in the translational dynamics of water molecules located close to the lanreotide surface. We also provide the molecular details of the self-assembly in the dimer in terms of structure, hydration, and energy.

Molecular dynamics and density functional theory studies of γ-butyrolactone (GBL) + ethanol and γ-valerolactone (GVL) + ethanol liquid mixtures

γ-Butyrolactone (GBL) and γ-valerolactone (GVL) are two of the lactones with wide application in industry. The interactions of these two lactones with alcohols have been a subject of experimental research. This article reports the theoretical investigation of GBL/ethanol and GVL/ethanol mixtures throughout the entire composition range using molecular dynamics and ab initio calculations. Molecular dynamics simulation using the optimized potentials for liquid simulations force field reproduces the experimental thermodynamics results quite well. Investigation of the liquid structure show that ethanol interacts with GBL and GVL molecules at specific sites, namely, almost collinearly with the lactone CO bond, below and above the lactone ring plane, and between the CC ring bonds of the lactone. While ab initio calculations largely support the interaction of ethanol with the lactones via O-H---O hydrogen bonds, hydrogen boding via lactone CH and ethanol O clearly have a significant role as well. The O-H---O interactions are dominated by the carbonyl oxygen. Dispersion interactions seem to play a more important role in stabilizing those interactions in which the hydrogen bonding is via the lactone ring oxygen. Generally, the two lactones behave very similarly in both liquid and gas phase studies, with minor differences attributed to more steric hindrance due the methyl group in γ-valerolactone.

A potential for molecular simulation of compounds with linear moieties

The harmonic angle bending potential is used in many force fields for (bio)molecular simulation. The force associated with this potential is discontinuous at angles close to 180°, which can lead to numeric instabilities. Angle bending of linear groups, such as alkynes or nitriles, or linear molecules, such as carbon dioxide, can be treated by a simple harmonic potential if we describe the fluctuations as a deviation from a reference position of the central atom, the position of which is determined by the flanking atoms. The force constant for the linear angle potential can be derived analytically from the corresponding force constant in the traditional potential. The new potential is tested on the properties of alkynes, nitriles, and carbon dioxide. We find that the angles of the linear groups remain about 2° closer to 180° using the new potential. The bond and angle force constants for carbon dioxide were tuned to reproduce the experimentally determined frequencies. An interesting finding was that simulations of liquid carbon dioxide under pressure with the new flexible model were stable only when explicitly modeling the long-range Lennard-Jones (LJ) interactions due to the very long-range nature of the LJ interactions (>1.7 nm). In the other tested liquids, we find that a Lennard-Jones cutoff of 1.1 nm yields similar results as the particle mesh Ewald algorithm for LJ interactions. Algorithmic factors influencing the stability of liquid simulations are discussed as well. Finally, we demonstrate that the linear angle potential can be used in free energy perturbation calculations.

Optimization of the Reax force field for the lithium-oxygen system using a high fidelity charge model.

Studies using molecular dynamics (MD) have long struggled to simulate the failure modes of materials, predicting unrealistically high ductility and failing to capture brittle fracture. The primary cause of this shortcoming is an inadequate description of bond breaking. While reactive force fields such as ReaxFF show improvements compared to traditional force fields, the charge models used yield unphysical partial charges, especially during dissociation of ionic bonds. This flaw may be remedied by using the atom-condensed Kohn-Sham density functional theory (DFT) approximated to a second order (ACKS2) charge model for determining partial charges. In this work, we present a new ACKS2-enabled Reax force field for fracture simulations of lithium oxide systems, which was obtained by training against an extensive set of DFT, multireference configuration interaction (MRCI), and MRCI+Q reference data using genetic optimization techniques. This new force field significantly improves the bond breaking behavior, but still cannot fully capture the brittle fracture in MD simulations, suggesting more research is needed to improve simulation of brittle fracture.

A Simple Method for Including Polarization Effects in Solvation Free Energy Calculations When Using Fixed-Charge Force Fields: Alchemically Polarized Charges.

The incorporation of polarizability in classical force-field molecular simulations is an ongoing area of research. We focus here on its application to hydration free energy simulations of organic molecules. In contrast to computationally complex approaches involving the development of explicitly polarizable force fields, we present herein a simple methodology for incorporating polarization into such simulations using standard fixed-charge force fields, which we call the alchemically polarized charges (APolQ) method. APolQ employs a standard classical alchemical free energy change simulation to calculate the free energy difference between a fully polarized solute particle in a condensed phase and its unpolarized state in a vacuum. APolQ can in principle be applied to any microscopically homogeneous system (e.g., pure or mixed solvents). We applied APolQ to hydration free energy data for a test set of 45 neutral solute molecules in the FreeSolv database and compared results obtained using three different water models (SPC/E, TIP3P, and OPC3) and using minimal basis iterative Stockholder (MBIS) and restrained electrostatic potential (RESP) partial charge methodologies. In comparison with AM1-BCC, we found that APolQ outperforms it for the test set. Despite our method using default GAFF parameters, the MBIS partial charges yield absolute average deviations 1.5–1.9 kJ mol–1 lower than using AM1 bond charge correction (AM1-BCC). We conjecture that this method can be further improved by fitting the Lennard-Jones and torsional parameters to partial charges derived using MBIS or RESP methodologies.

The Protection of Lithium Metal Enabled by LiNO3 for Lithium-Sulfur Batteries

An insight into the role of LiNO3 additive in an ether-based electrolyte for lithium sulfur battery has been presented. Herein, we proposed a formation mechanism of solid-electrolyte interphase on the surface of lithium metal anode by using the theoretical reactive force field (ReaxFF) simulation method. The interaction between the reactive lithium metal and nitrate ions results in the formation of LixNOy clusters, distributed homogeneously in both the SEI layer and electrolyte phase. Not only can these clusters be an efficient surface protection layer of lithium metal anode but also, they can suppress the lithium polysulfides shuttle effect by adsorbing liquid lithium polysulfide intermediates.

Hybrid Combination of Quantum Mechanics with Quantum-Based Polarizable Reactive Force Field for Large Scale Full Solvent Simulations of Electrocatalysis

To develop new generations of electrocatalysts required for energy and environmental sustainability, we need the accuracy of full solvent quantum mechanics (QM) (free energy barriers to 0.05 eV, onset potentials to 0.05 V) but for practical sized nanoparticles and catalysts (1000’s to millions of atoms). We report here a solution to this problem. We start with the RexPoN reactive force field that provides higher accuracy than density functional theory (DFT) and combine it with QM to accurately include long-range interactions and polarization effects to enable reactive simulations with QM accuracy in the presence of solvent including 1000’s to millions of waters. Here we apply this RexPoN embedded QM (ReQM) to reactive simulations of electrocatalysis demonstrating that ReQM accurately replaces DFT water for computing the Raman frequencies of reaction intermediates during CO2 reduction to ethylene, with comparisons to operando electrocatalysis experiments and to full solvent QM calculations.