Standard Force Field Research Articles

Insight into uranyl binding by cyclic peptides from molecular dynamics and density functional theory

It is a challenging task to develop uranyl-chelating agents based on peptide chemistry. A recently developed cationic dummy atom model of uranyl in conjunction with the classical molecular dynamics simulation presents a helpful utility to study the chelation of uranyl by peptides with a low computational cost. In the present study, it was used to describe the chelation of uranyl by the cyclic decapeptide with 4 Glu residues cyc-GluArgGluProGlyGluTrpGluProGly and its derivatives containing two phosphorylated serines in place of two Glu, termed pS16, pS18, pS38, and pS68. The obtained structures were further studied by density functional theory (DFT) and subsequent density analysis. We show that a combination of steered molecular dynamics and simulated annealing, using standard forcefields for peptide with the cationic dummy atom model of uranyl, can quickly and reliably obtain binding modes of uranyl-peptide complexes. Classical molecular dynamics simulation in explicit water produces geometry very close to the DFT-optimized structure. The presence of uranyl completely changes the conformation of these cyclic peptides from unstructured to organised. The simulation of a peptide with two uranyl units explained why only the 1: 1 ratio of peptide and chelated-uranyl is observed experimentally in most cases, by the insufficiency of the anionic residues for the chelation of two UO22+ units, but that pS16 can accommodate two such units.

NP-Collabs: Investigating Counterion-Mediated Bridges in the Multiply Phosphorylated Tau-R2 Repeat.

Tau is an intrinsically disordered (IDP) microtubule-associated protein (MAP) that plays a key part in microtubule assembly and organization. The function of tau can be regulated by multiple phosphorylation sites. These post-translational modifications are known to decrease the binding affinity of tau for microtubules, and abnormal tau phosphorylation patterns are involved in Alzheimer's disease. Using all-atom molecular dynamics simulations, we compared the conformational landscapes explored by the tau R2 repeat domain (which comprises a strong tubulin binding site) in its native state and with multiple phosphorylations on the S285, S289, and S293 residues, with four different standard force field (FF)/water model combinations. We find that the different parameters used for the phosphate groups (which can be more or less flexible) in these FFs and the specific interactions between bulk cations and water lead to the formation of a specific type of counterion bridge, termed nP-collab (for nphosphate collaboration, with n being an integer), where counterions form stable structures binding with two or three phosphate groups simultaneously. The resulting effect of nP-collabs on the tau-R2 conformational space differs when using sodium or potassium cations and is likely to impact the peptide overall dynamics and how this MAP interacts with tubulins. We also investigated the effect of phosphoresidue spacing and ionic concentration by modeling polyalanine peptides containing two phosphoserines located one-six residues apart. Three new metrics specifically tailored for IDPs (proteic Menger curvature, local curvature, and local flexibility) were introduced, which allow us to fully characterize the impact of nP-collabs on the dynamics of disordered peptides at the residue level.

Adjusting the Energy Profile for CH-O Interactions Leads to Improved Stability of RNA Stem-Loop Structures in MD Simulations.

The role of ribonucleic acid (RNA) in biology continues to grow, but insight into important aspects of RNA behavior is lacking, such as dynamic structural ensembles in different environments, how flexibility is coupled to function, and how function might be modulated by small molecule binding. In the case of proteins, much progress in these areas has been made by complementing experiments with atomistic simulations, but RNA simulation methods and force fields are less mature. It remains challenging to generate stable RNA simulations, even for small systems where well-defined, thermostable structures have been established by experiments. Many different aspects of RNA energetics have been adjusted in force fields, seeking improvements that are transferable across a variety of RNA structural motifs. In this work, the role of weak CH···O interactions is explored, which are ubiquitous in RNA structure but have received less attention in RNA force field development. By comparing data extracted from high-resolution RNA crystal structures to energy profiles from quantum mechanics and force field calculations, it is shown that CH···O interactions are overly repulsive in the widely used Amber RNA force fields. A simple, targeted adjustment of CH···O repulsion that leaves the remainder of the force field unchanged was developed. Then, the standard and modified force fields were tested using molecular dynamics (MD) simulations with explicit water and salt, amassing over 300 μs of data for multiple RNA systems containing important features such as the presence of loops, base stacking interactions as well as canonical and noncanonical base pairing. In this work and others, standard force fields lead to reproducible unfolding of the NMR-based structures. Including a targeted CH···O adjustment in an otherwise identical protocol dramatically improves the outcome, leading to stable simulations for all RNA systems tested.

FeNNol: An efficient and flexible library for building force-field-enhanced neural network potentials.

Neural network interatomic potentials (NNPs) have recently proven to be powerful tools to accurately model complex molecular systems while bypassing the high numerical cost of abinitio molecular dynamics simulations. In recent years, numerous advances in model architectures as well as the development of hybrid models combining machine-learning (ML) with more traditional, physically motivated, force-field interactions have considerably increased the design space of ML potentials. In this paper, we present FeNNol, a new library for building, training, and running force-field-enhanced neural network potentials. It provides a flexible and modular system for building hybrid models, allowing us to easily combine state-of-the-art embeddings with ML-parameterized physical interaction terms without the need for explicit programming. Furthermore, FeNNol leverages the automatic differentiation and just-in-time compilation features of the Jax Python library to enable fast evaluation of NNPs, shrinking the performance gap between ML potentials and standard force-fields. This is demonstrated with the popular ANI-2x model reaching simulation speeds nearly on par with the AMOEBA polarizable force-field on commodity GPUs (graphics processing units). We hope that FeNNol will facilitate the development and application of new hybrid NNP architectures for a wide range of molecular simulation problems.

Quantifying the Uncertainty of Force Field Selection on Adsorption Predictions in MOFs.

Comparisons between simulated and experimental adsorption isotherms in MOFs are fraught with challenges. On the experimental side, there is significant variation between isotherms measured on the same system, with a significant percentage (∼20%) of published data being considered outliers. On the simulation side, force fields are often chosen "off-the-shelf" with little or no validation. The effect of this choice on the reliability of simulated adsorption predictions has not yet been rigorously quantified. In this work, we fill this gap by systematically quantifying the uncertainty arising from force field selection on adsorption isotherm predictions. We choose methane adsorption, where electrostatic interactions are negligible, to independently study the effect of the framework Lennard-Jones parameters on a series of prototypical materials that represent the most widely studied MOF "families". Using this information, we compute an adsorption "consensus isotherm" from simulations, including a quantification of uncertainty, and compare it against a manually curated set of experimental data from the literature. By considering many experimental isotherms measured by different groups and eliminating outliers in the data using statistical analysis, we conduct a rigorous comparison that avoids the pitfalls of the standard approach of comparing simulation predictions to a single experimental data set. Our results show that (1) the uncertainty in simulated isotherms can be as large as 15% and (2) standard force fields can provide reliable predictions for some systems but can fail dramatically for others, highlighting systematic shortcomings in those models. Based on this, we offer recommendations for future simulation studies of adsorption, including high-throughput computational screening of MOFs.

Non-adiabatic direct quantum dynamics using force fields: Toward solvation.

Quantum dynamics simulations are becoming a powerful tool for understanding photo-excited molecules. Their poor scaling, however, means that it is hard to study molecules with more than a few atoms accurately, and a major challenge at the moment is the inclusion of the molecular environment. Here, we present a proof of principle for a way to break the two bottlenecks preventing large but accurate simulations. First, the problem of providing the potential energy surfaces for a general system is addressed by parameterizing a standard force field to reproduce the potential surfaces of the molecule's excited-states, including the all-important vibronic coupling. While not shown here, this would trivially enable the use of an explicit solvent. Second, to help the scaling of the nuclear dynamics propagation, a hierarchy of approximations is introduced to the variational multi-configurational Gaussian method that retains the variational quantum wavepacket description of the key quantum degrees of freedom and uses classical trajectories for the remaining in a quantum mechanics/molecular mechanics like approach. The method is referred to as force field quantum dynamics (FF-QD), and a two-state ππ*/nπ* model of uracil, excited to its lowest bright ππ* state, is used as a test case.

Potential energy landscape of a flexible water model: Equation of state, configurational entropy, and Adam-Gibbs relationship.

The potential energy landscape (PEL) formalism is a tool within statistical mechanics that has been used in the past to calculate the equation of states (EOS) of classical rigid model liquids at low temperatures, where computer simulations may be challenging. In this work, we use classical molecular dynamics (MD) simulations and the PEL formalism to calculate the EOS of the flexible q-TIP4P/F water model. This model exhibits a liquid-liquid critical point (LLCP) in the supercooled regime, at (Pc = 150 MPa, Tc = 190 K, and ρc = 1.04 g/cm3) [using the reaction field technique]. The PEL-EOS of q-TIP4P/F water and the corresponding location of the LLCP are in very good agreement with the MD simulations. We show that the PEL of q-TIP4P/F water is Gaussian, which allows us to calculate the configurational entropy of the system, Sconf. The Sconf of q-TIP4P/F water is surprisingly similar to that reported previously for rigid water models, suggesting that intramolecular flexibility does not necessarily add roughness to the PEL. We also show that the Adam-Gibbs relation, which relates the diffusion coefficient D with Sconf, holds for the flexible q-TIP4P/F water model. Overall, our results indicate that the PEL formalism can be used to study molecular systems that include molecular flexibility, the common case in standard force fields. This is not trivial since the introduction of large bending/stretching mode frequencies is problematic in classical statistical mechanics. For example, as shown previously, we find that such high frequencies lead to unphysical (negative) entropy for q-TIP4P/F water when using classical statistical mechanics (yet, the PEL formalism can be applied successfully).

A Simple Electron-Density Based Force Field Model for High-Energy Interactions between Atoms and Molecules.

In high-energy molecular dynamics or Monte Carlo simulations, standard force fields optimized for simulations at ambient temperatures are inadequate. This is largely because their repulsive parts have been regarded as not very significant, even well below zero interaction energies. It is, therefore, not obvious which force fields to resort to for simulating hot gases or plasmas. A force field model that uses the electronic densities of noninteracting atoms or molecules within the pair approximation is introduced. We start by deriving a naïve model that neglects any exchange and correlation effects between the electronic clouds and then correct this model by adding a term calibrated from ab initio calculations using the CCSD(T)/cc-pVTZ level of theory. The resulting expression for this term can be regarded as a simple exchange-correlation function. We compare the results for the repulsive part of the potential energy hypersurfaces with the force fields commonly used on some dimers of small molecules.

Interatomic force fields for zirconium based on the embedded atom method and the tabulated Gaussian Approximation Potential

The accuracy of interatomic interaction potentials – also known as force fields – is the main factor determining the physical soundness of classical molecular dynamics (MD) simulations. Here, we present a multi-objective framework to generate embedded-atom-method (EAM) force fields using ab initio data. The EAM force fields were tuned via particle swarm optimization to capture the non-linear association between atomic structures and system energies. Using this framework, 95 standard EAM force fields for zirconium were developed and 45 physical features for each developed force field were tracked. Principal component analysis (PCA) was performed to provide insights into the compromises that must be made when generating EAM force fields. Of note, by assigning large fitting weights to generalized stacking fault energy (GSFE) surfaces, there exist EAM force fields with properly positioned minima on prismatic GSFE surfaces and containing no spurious minima in basal GSFE surfaces. However, while standard EAM force fields achieved this without explicitly taking the angular dependence of atomic interactions into account, they led to a severe mismatch between other important physical properties and benchmarks. Hence, we also constructed two machine-learned tabulated Gaussian approximation potentials (tabGAP) with an additional three-body term, which successfully tackled the aforementioned issue and exhibit acceptable prediction accuracy across many physical properties (lattice parameters, elastic properties, dimer potential energies, melting temperatures, phase stability, point defect formation energies, point defect migration energies, and free surface energies) of Zr. Remarkably, its computational efficiency is only 6 times slower than standard EAM force fields.

OPLS/2020 Force Field for Unsaturated Hydrocarbons, Alcohols, and Ethers.

The OPLS all-atom force field was updated and applied to modeling unsaturated hydrocarbons, alcohols, and ethers. Testing has included gas-phase conformational energetics, properties of pure liquids, and free energies of hydration. Monte Carlo statistical mechanics (MC) calculations were used to model 60 liquids. In addition, a robust, automated procedure was devised to compute the free energies of hydration with high precision via free-energy perturbation (FEP) calculations using double annihilation. Testing has included larger molecules than in the past, and parameters are reported for the first time for some less common groups including alkynes, allenes, dienes, and acetals. The average errors in comparison with experimental data for the computed properties of the pure liquids were improved with the modified force field (OPLS/2020). For liquid densities and heats of vaporization, the average unsigned errors are 0.01 g/cm3 and 0.2 kcal/mol. The average error and signed error for free energies of hydration are both 1.2 kcal/mol. As noted before, this reflects a systematic overestimate of the hydrophobicity of organic molecules when the parametrization is done to minimize the errors for properties of pure liquids. Implications for the modeling of biomolecular systems with standard force fields are considered.

Gaussian attractive potential for carboxylate/cobalt surface interactions.

Ligand-decorated metal surfaces play a pivotal role in various areas of chemistry, particularly in selective catalysis. Molecular dynamics simulations at the molecular mechanics level of theory are best adapted to gain complementary insights to experiments regarding the structure and dynamics of such organic films. However, standard force fields tend to capture only weak physisorption interactions. This is inadequate for ligands that are strongly adsorbed such as carboxylates on metal surfaces. To address this limitation, we employ the Gaussian Lennard-Jones (GLJ) potential, which incorporates an attractive Gaussian potential between the surface and ligand atoms. Here, we develop this approach for the interaction between cobalt surfaces and carboxylate ligands. The accuracy of the GLJ approach is validated through the analysis of the interaction of oxygen with two distinct cobalt surfaces. The accuracy of this method reaches a root mean square deviation (RMSD) of about 3 kcal/mol across all probed configurations, which corresponds to a percentage error of roughly 4%. Application of the GLJ force field to the dynamics of the organic layer on these surfaces reveals how the ligand concentration influences the film order, and highlights differing mobility in the x and y directions, attributable to surface corrugation on Co(112̄0). GLJ is versatile, suitable for a broad range of metal/ligand systems, and can, subsequently, be utilized to study the organic film on the adsorption/desorption of reactants and products during a catalytic process.

Parameter Dependence of the Solubility Limit for Disodium Phosphate.

The solubility limit was calculated for supersaturated solutions of disodium phosphate in water as a function of the sodium-oxygen Lennard-Jones radius parameter Rmin. We found that changes in the sodium-oxygen Rmin were clearly exponentially related to the concentration of the solubility limit. Starting from standard force fields more suited to nucleic acids and phospholipids, only relatively small changes were required to achieve the experimentally known solubility limit. Simultaneously, we found that it was possible to achieve the solubility limit and the osmotic pressure with the same model parameters. Based on transferability, the adjusted Rmin parameter can be used to more accurately model phosphorylated proteins.

AutoParams: An Automated Web-Based Tool To Generate Force Field Parameters for Molecular Dynamics Simulations.

Many research questions benefit from molecular dynamics simulations to observe the motions and conformations of molecules over time, which rely on force fields that describe sets of common molecules by category. With the increase of importance for large data sets used in machine learning and growing computational efficiency, the ability to rapidly create large numbers of force field inputs is of high importance. Unusual molecules, such as nucleotide analogues, functionalized carbohydrates, and modified amino acids, are difficult to describe consistently using standard force fields, requiring the development of custom parameters for each unique molecule. While these parameters may be created by individual users, the process can become time-consuming or may introduce errors that may not be immediately apparent. We present an open-source automated parameter generation service, AutoParams, which requires minimal input from the user and creates useful Amber force field parameter sets for most molecules, particularly those that combine molecular types (e.g., a carbohydrate functionalized with a benzene). We include hierarchical atom-typing logic that makes it straightforward to expand with additional force fields and settings, and options for creating monomers in polymers, such as functionalized amino acids. It can be straightforwardly linked to any charge generation program and currently has interfaces to Psi4, PsiRESP, and TeraChem. It is open source and is available via GitHub. It includes error checking and testing protocols to ensure the parameters will be sufficient for subsequent molecular dynamics simulations and streamlines the creation of force field databases.

Upgrading of the general AMBER force field 2 for fluorinated alcohol biosolvents: A validation for water solutions and melittin solvation.

The standard GAFF2 force field parameterization has been refined for the fluorinated alcohols 2,2,2-trifluoroethanol (TFE), 1,1,1,3,3,3-hexafluoro-2-propanol (HFIP), and 1,1,1,3,3,3-hexafluoropropan-2-one (HFA), which are commonly used to study proteins and peptides in biomimetic media. The structural and dynamic properties of both proteins and peptides are significantly influenced by the biomimetic environment created by the presence of these cosolvents in aqueous solutions. Quantum mechanical calculations on stable conformers were used to parameterize the atomic charges. Different systems, such as pure liquids, aqueous solutions, and systems formed by melittin protein and cosolvent/water solutions, have been used to validate the new models. The calculated macroscopic and structural properties are in agreement with experimental findings, supporting the validity of the newly proposed models.

Simulating water adsorption in metal–organic frameworks with open metal sites using the 12-6-4 Lennard–Jones potential

ABSTRACT Metal–organic frameworks (MOFs) with coordinatively unsaturated open metal sites (OMSs) are promising sorbent materials in atmospheric water harvesting (AWH) systems at low relative humidity (RH) due to their strong interactions with water molecules. However, accurate computational modelling of water adsorption in those materials are challenging, as standard force fields (FFs) based on the 12-6 Lennard–Jones (L–J) potential cannot properly describe the water–OMS interactions. In biomolecular simulations, the 12-6-4 L–J potential has been successfully used to model metal ion–water interactions in which the extra 1/r 4-term is for charge-induced dipole electrostatics. In this work, we adopted this strategy and used the 12-6-4 L–J potential to model water–OMS interactions in MOFs for low RH AWH applications. Without any modifications, the parameters used for aqueous metal ions were able to greatly improve the accuracy in predicting water adsorption isotherms in MOF-74, which highlights the simplicity and transferability of this method.

An implementation of the Martini coarse-grained force field in OpenMM

We describe a complete implementation of Martini 2 and Martini 3 in the OpenMM molecular dynamics software package. Martini is a widely used coarse-grained force field with applications in biomolecular simulation, materials, and broader areas of chemistry. It is implemented as a force field but makes extensive use of facilities unique to the GROMACS software, including virtual sites and bonded terms that are not commonly used in standard atomistic force fields. OpenMM is a flexible molecular dynamics package widely used for methods development and is competitive in speed on GPUs with other commonly used packages. OpenMM has facilities to easily implement new force field terms, external forces and fields, and other nonstandard features, which we use to implement all force field terms used in Martini 2 and Martini 3. This allows Martini simulations, starting with GROMACS topology files that are processed by custom scripts, with all the added flexibility of OpenMM. We provide a GitHub repository with test cases, compare accuracy and performance between GROMACS and OpenMM, and discuss the limitations of our implementation in terms of direct comparison with GROMACS. We describe a use case that implements the Modeling Employing Limited Data method to apply experimental constraints in a Martini simulation to efficiently determine the structure of a protein complex. We also discuss issues and a potential solution with the Martini 2 topology for cholesterol.

IDP Force Fields Applied to Model PPII-Rich 33-mer Gliadin Peptides.

The 33-mer gliadin peptide and its deamidated metabolite, 33-mer DGP, are the immunodominant peptides responsible for the adaptive immune response in celiac disease (CD). CD is a complex autoimmune chronic disorder triggered by gluten ingestion that affects the small intestine and affects ∼1% of the global population. The 33-mers are polyproline II-rich (PPII) and intrinsically disordered peptides (IDPs), whose structures remain elusive. We sampled the conformational ensembles of both 33-mer peptides via molecular dynamics simulations employing two force fields (FFs) (Amber ff03ws and Amber ff99SB-disp) specifically validated for other IDPs. Our results show that both FFs allow the extensive exploration of the conformational landscape, which was not possible with the standard FF GROMOS53A6 reported before. Clustering analysis of the trajectories showed that the five largest clusters (78-88% of the total structures) present elongated, semielongated, and curved conformations in both FFs. Large average radius of gyration and solvent-exposed surfaces characterized these structures. While the structures sampled are similar, the Amber ff99SB-disp trajectories explored folded conformations with a higher probability. In addition, PPII secondary structure was preserved throughout the trajectories (58-73%) together with a non-negligible content of β structures (11-23%), in agreement with previous experimental results. This work represents the initial step in studying further the interaction of these peptides with other biologically relevant molecules, which could lead to finally disclose the molecular events that lead to CD.

Quantum Mechanical-Cluster Approach to Solve the Bioisosteric Replacement Problem in Drug Design.

Bioisosteres are molecules that differ in substituents but still have very similar shapes. Bioisosteric replacements are ubiquitous in modern drug design, where they are used to alter metabolism, change bioavailability, or modify activity of the lead compound. Prediction of relative affinities of bioisosteres with computational methods is a long-standing task; however, the very shape closeness makes bioisosteric substitutions almost intractable for computational methods, which use standard force fields. Here, we design a quantum mechanical (QM)-cluster approach based on the GFN2-xTB semi-empirical quantum-chemical method and apply it to a set of H → F bioisosteric replacements. The proposed methodology enables advanced prediction of biological activity change upon bioisosteric substitution of -H with -F, with the standard deviation of 0.60 kcal/mol, surpassing the ChemPLP scoring function (0.83 kcal/mol), and making QM-based ΔΔG estimation comparable to ∼0.42 kcal/mol standard deviation of in vitro experiment. The speed of the method and lack of tunable parameters makes it affordable in current drug research.

Mechanistic understanding of the correlation between structure and dynamics of liquid carbonate electrolytes: impact of polarization.

Liquid electrolyte design and modelling is an essential part of the development of improved lithium ion batteries. For mixed organic carbonates (ethylene carbonate (EC) and ethyl-methyl carbonate (EMC) mixtures)-based electrolytes with LiPF6 as salt, we have compared a polarizable force field with the standard non-polarizable force field with and without charge rescaling to model the structural and dynamic properties. The result of our molecular dynamics simulations shows that both polarizable and non-polarizable force fields have similar structural factors, which are also in agreement with X-ray diffraction experimental results. In contrast, structural differences are observed for the lithium neighborhood, while the lithium-anion neighbourhood is much more pronounced for the polarizable force field. Comparison of EC/EMC coordination statistics with Fourier transformed infrared spectroscopy (FTIR) shows the best agreement for the polarizable force field. Also for transport quantities such as ionic conductivities, transference numbers, and viscosities, the agreement with the polarizable force field is by far better for a large range of salt concentrations and EC : EMC ratios. In contrast, for the non-polarizable variants, the dynamics are largely underestimated. The excellent performance of the polarizable force field is explored in different ways to pave the way to a realistic description of the structure-dynamics relationships for a wide range of salt and solvent compositions for this standard electrolyte. In particular, we can characterize the distinct correlation terms between like and unlike ions, relate them to structural properties, and explore to which degree the transport in this electrolyte is mass or charge limited.

Charge-transfer-to-solvent states provide a sensitive spectroscopic probe of the local solvent structure around anions

ABSTRACT This computational study characterises charge-transfer-to-solvent (CTTS) states of aqueous thiocyanate anion using equation-of-motion coupled-cluster methods combined with electrostatic embedding quantum mechanics/molecular mechanics (QM/MM) scheme. Equilibrium sampling was carried out using classical molecular dynamics (MD) with standard force-fields and QM/MM ab initio molecular dynamics (AIMD) using density functional theory. The two calculations yield significantly different local structure around solvated SCN − . Because of the diffuse character of CTTS states, they are very sensitive to the local structure of solvent around the solute and its dynamic fluctuations. Owing to this sensitivity, the spectra computed using MD and AIMD based snapshots differ considerably. This sensitivity suggests that the spectroscopy exploiting CTTS transitions can provide an experimental handle for assessing the quality of force-fields and density functionals. By combining CTTS-based spectroscopies with reliable theoretical modeling, detailed microscopic information of the solvent structure can be obtained. We present a robust computational protocol for modeling spectra of solvated anions and emphasise the use of an ab initio characterization of individual electronic transitions as CTTS or local excitations.