GAFF Force Field Research Articles

Comparative assessment of physics-based in silico methods to calculate relative solubilities.

Relative solubilities, i.e. whether a given molecule is more soluble in one solvent compared to others, is a critical parameter for pharmaceutical and agricultural formulation development and chemical synthesis, material science, and environmental chemistry. In silico predictions of this crucial variable can help reducing experiments, waste of solvents and synthesis optimization. In this study, we evaluate the performance of different physics-based methods for predicting relative solubilities. Our assessment involves quantum mechanics-based COSMO-RS and molecular dynamics-based free energy methods using OPLS4, the open-source OpenFF Sage, and GAFF force fields, spanning over 200 solvent-solute combinations. Our investigation highlights the important role of compound multimerization, an effect which must be accounted for to obtain accurate relative solubility predictions. The performance landscape of these methods is varied, with significant differences in precision depending on both the method used and the solute considered, thereby offering an improved understanding of the predictive power of physics-based methods in chemical research.

Inclusion of Water Multipoles into the Implicit Solvation Framework Leads to Accuracy Gains.

The current practical "workhorses" of the atomistic implicit solvation─the Poisson-Boltzmann (PB) and generalized Born (GB) models─face fundamental accuracy limitations. Here, we propose a computationally efficient implicit solvation framework, the Implicit Water Multipole GB (IWM-GB) model, that systematically incorporates the effects of multipole moments of water molecules in the first hydration shell of a solute, beyond the dipole water polarization already present at the PB/GB level. The framework explicitly accounts for coupling between polar and nonpolar contributions to the total solvation energy, which is missing from many implicit solvation models. An implementation of the framework, utilizing the GAFF force field and AM1-BCC atomic partial charges model, is parametrized and tested against the experimental hydration free energies of small molecules from the FreeSolv database. The resulting accuracy on the test set (RMSE ∼ 0.9 kcal/mol) is 12% better than that of the explicit solvation (TIP3P) treatment, which is orders of magnitude slower. We also find that the coupling between polar and nonpolar parts of the solvation free energy is essential to ensuring that several features of the IWM-GB model are physically meaningful, including the sign of the nonpolar contributions.

GAFF-Based Polarizable Force Field Development and Validation for Ionic Liquids.

Ionic liquids (ILs) have been used in many applications, including gas separations, electrochemistry, lubrication, and catalysis. Understanding how the different properties of ILs are related to their chemical structure and composition is crucial for these applications. Experimental investigations often provide limited insights and can be tedious in exploring a range of state points. Therefore, molecular simulations have emerged as a powerful tool that not only offers a microscopic perspective but also enables rapid screening and prediction of physical properties. The accuracy of these predictions, however, depends on the quality of the intermolecular potentials (force fields) used. The widely used classical fixed charge models, such as GAFF, OPLS, and CL&P, are popular due to their simplicity and computational efficiency. However, it has been shown that the use of integer charges with these classical models leads to sluggish dynamics. The use of scaled charge models can improve the dynamics, but these mean-field approaches are unable to account for polarization effects explicitly. Several different approaches have been proposed to include polarizability in IL force fields. In this work, we follow the protocol of the CL&Pol model to develop a Drude oscillator model based on the GAFF force field (Goloviznina, K., et al. J. Chem. Theory Comput. 2019, 15, 5858). We compare the performance of the model for eight imidazolium- and pyrrolidinium-based ILs against that of other models. We find that the new model provides reasonable estimations of density, self-diffusivity, and structural properties for these ILs and suggests a relatively simple way of extending the general GAFF model to more ILs.

Assessment of GAFF and OPLS Force Fields for Urea: Crystal and Aqueous Solution Properties.

Molecular simulations such as Monte Carlo, molecular dynamics, and metadynamics have been used to provide insight into crystallization phenomena, including nucleation and crystal growth. However, these simulations depend on the force field used, which models the atomic and molecular interactions, to adequately reproduce relevant material properties for the phases involved. Two widely used force fields, the General AMBER Force Field (GAFF) and the Optimized Potential for Liquid Simulations (OPLS), including several variants, have previously been used for studying urea crystallization. In this work, we investigated how well four different versions of the GAFF force field and five different versions of the OPLS force field reproduced known urea crystal and aqueous solution properties. Two force fields were found to have the best overall performance: a specific urea charge-optimized GAFF force field and the original all-atom OPLS force field. It is recommended that a suitable testing protocol involving both solution and solid properties, such as that used in this work, is adopted for the validation of force fields used for simulations of crystallization phenomena.

Structure and interactions in reline (1:2 Choline chloride – Urea mixture) according to quantum chemical calculations and molecular dynamics simulation

Quantum chemical calculations using the B3LYP + GD3/6–311++g(2d,2p) method were carried out to study the structure and interactions in a number of complexes containing a choline cation, a chloride anion and urea molecules. The contributions from Coulomb and van der Waals interactions were calculated. The QTAIM method was used to evaluate the strength of various bonds in these complexes. A molecular dynamics simulation of reline was also carried out in the NPT ensemble at a pressure of 1 atm and at temperatures of 298 and 323 K, with the GAFF force field and charges taken from our quantum chemical calculations. The particle radial distribution functions of centers of mass, local mole fraction functions, atom–atom radial distribution functions, and spatial distribution functions were calculated. Particles of different types were shown to form bonds of different types competing with each other.

Exploring potential SARS-CoV-2 Mpro non-covalent inhibitors through docking, pharmacophore profile matching, molecular dynamic simulation, and MM-GBSA.

In the replication of SARS-CoV-2, the main protease (Mpro/3CLpro) is significant. It is conserved in a number of novel coronavirus variations, and no known human proteases share its cleavage sites. Therefore, 3CLpro is an ideal target. In the report, we screened five potential inhibitors (1543, 2308, 3717, 5606, and 9000) of SARS-CoV-2 Mpro through a workflow. The calculation of MM-GBSA binding free energy showed that three of the five potential inhibitors (1543, 2308, 5606) had similar inhibitor effects to X77 against Mpro of SARS-CoV-2. In conclusion, the manuscript lays the groundwork for the design of Mpro inhibitors. In the virtual screening phase, we used structure-based virtual screening (Qvina2.1) and ligand-based virtual screening (AncPhore). In the molecular dynamic simulation part, we used the Amber14SB + GAFF force field to perform molecular dynamic simulation of the complex for 100ns (Gromacs2021.5) and performed MM-GBSA binding free energy calculation according to the simulation trajectory.

Role of Non-Covalent Interactions in Carbonic Anhydrase I-Topiramate Complex Based on QM/MM Approach.

Carbonic anhydrase (CA) I with a Topiramate (TPM) complex was investigated on the basis of a Quantum Mechanics/Molecular Mechanics (QM/MM) approach. The QM part was treated using Density Functional Theory (DFT) while the MM was simulated using Amberff14SB and GAFF force fields. In addition, the TIP3P model was applied to reproduce the polar environment influence on the studied complex. Next, three snapshots (after 5 ps, 10 ps, and 15 ps of the simulation time) were taken from the obtained trajectory to provide an insight into the non-covalent interactions present between the ligand and binding pocket of the protein. Our special attention was devoted to the binding site rearrangement, which is known in the literature concerning the complex. This part of the computations was performed using ωB97X functional with Grimme D3 dispersion corrections as well as a Becke-Johnson damping function (D3-BJ). Two basis sets were applied: def2-SVP (for larger models) and def2-TZVPD (for smaller models), respectively. In order to detect and describe non-covalent interactions between amino acids of the binding pocket and the ligand, Independent Gradient Model based on Hirshfeld partitioning (IGMH), Interaction Region Indicator (IRI), Quantum Theory of Atoms in Molecules (QTAIM) and Natural Bond Orbitals (NBO) methods were employed. Finally, Symmetry-Adapted Perturbation Theory (SAPT) was applied for energy decomposition between the ligand and protein. It was found that during the simulation time, the ligand position in the binding site was preserved. Nonetheless, amino acids interacting with TPM were exchanging during the simulation, thus showing the binding site reorganization. The energy partitioning revealed that dispersion and electrostatics are decisive factors that are responsible for the complex stability.

The reconciliation between the experimental and calculated octanol-water partition coefficient of 1,2-dipalmitoyl-sn-glycero-3-phosphatidylcholine using atomistic molecular dynamics: an open question

The octanol-water partition coefficient of 1,2-dipalmitoyl-sn-glycero-3-phosphatidylcholine (DPPC) was investigated using atomistic molecular dynamics simulations via thermodynamic integration and multistate Bennett acceptance ratio methods. The GAFF and CHARMM36 force fields were used with six water models widely used in molecular dynamics simulations. The OPC4 water model provided the best agreement with the experimental octanol-water partition coefficient of DPPC using the two force fields. However, there is still plenty of room for improvement in water models with correct estimation of surface tension that uses better and suitable non-bonded interaction parameters between water-water and water-DPPC. The Gibbs free energy of transferring DPPC from octanol to water phase was calculated to be 19.8 ± 0.3 and 20.2 ± 0.3 kcal mol−1, giving a partition coefficient of 14.5 ± 0.4 and 14.8 ± 0.3 for the GAFF and CHARMM36 force fields, respectively. This study reinforces the importance of developing new water models that reproduce experimental surface tensions to reconcile the water-water and water-DPPC non-bonded interactions and the existing discrepancy between experimental measurements of amphiphilic molecules that are important in many areas of scientific applications and industry such as biophysics, surfactant, colloids, membranes, medicine, nanotechnology, and food and pharmaceutical industries, and so on. It raises two important open questions: Is the experimental octanol-water partition coefficient of DPPC reliable? Or is its calculation accurate using the OPC4 water model? With respect to the experimental measurements, there may be non-treated aspects such as the formation of aggregates in aqueous phase and limit of detection of the applied method. And, in the calculation, some effects are not possible to be considered in a correct way or viable time such as calculating quantum effects, sampling all conformations, considering phase transitions, and correctly evaluating the intermolecular forces to estimate an accurate surface tension. Communicated by Ramaswamy H. Sarma

Parameterization and Application of the General Amber Force Field to Model Fluoro Substituted Furanose Moieties and Nucleosides.

Molecular mechanics force field calculations have historically shown significant limitations in modeling the energetic and conformational interconversions of highly substituted furanose rings. This is primarily due to the gauche effect that is not easily captured using pairwise energy potentials. In this study, we present a refinement to the set of torsional parameters in the General Amber Force Field (gaff) used to calculate the potential energy of mono, di-, and gem-fluorinated nucleosides. The parameters were optimized to reproduce the pseudorotation phase angle and relative energies of a diverse set of mono- and difluoro substituted furanose ring systems using quantum mechanics umbrella sampling techniques available in the IpolQ engine in the Amber suite of programs. The parameters were developed to be internally consistent with the gaff force field and the TIP3P water model. The new set of angle and dihedral parameters and partial charges were validated by comparing the calculated phase angle probability to those obtained from experimental nuclear magnetic resonance experiments.

Influence of Temperature and Salt Concentration on the Hydrophobic Interactions of Adamantane and Hexane

One of the definitions of hydrophobic interactions is the aggregation of nonpolar particles in a polar solvent, such as water. While this phenomenon appears to be very simple, it is crucial for many complex processes, such as protein folding, to take place. In this work, the hydrophobic association of adamantane and hexane at various temperatures and ionic strengths was studied using molecular dynamics simulations with the AMBER 16.0 program and the GAFF force field. The potentials of mean force of hydrophobic dimer formation, as well as the excess free energy, excess energy, excess entropy, and excess heat capacity corresponding to the formation of the contact minimum, were determined and analyzed. For both systems, the depth of the contact minimum in the potential of mean force was found to increase with both temperature and ionic strength. The excess heat capacity of the association at the contact minimum and T = 298 K was found to be negative and to decrease, while the excess entropy and energy were found to be positive and to increase for both systems, the changes being more pronounced for the hexane dimer. The excess heat capacity is also greater in absolute value for the hexane dimer.

Binding free energies for the SAMPL8 CB8 \u201cDrugs of Abuse\u201d challenge from umbrella sampling combined with Hamiltonian replica exchange

Umbrella sampling along a one-dimensional order parameter in combination with Hamiltonian replica exchange was employed to calculate the binding free energy of five guest molecules with known affinity to cucurbit[8]uril. A simple empirical approach correcting for the overestimation of the affinity by the GAFF force field was proposed and subsequently applied to the seven guest molecules of the “Drugs of Abuse” SAMPL8 challenge. Compared to the uncorrected binding free energies, the systematic error decreased but quantitative agreement with experiment was only reached for a few compounds. From a retrospective analysis a weak point of the correction term was identified.

Comparison of the GAFF, OPLSAA and CHARMM27 force field for the reproduction of the thermodynamics properties of furfural, 2-methylfuran, 2,5-dimethylfuran and 5-hydroxymethylfurfural

The biofuels 2-methylfuran and 2,5-dimethylfuran, derived from furfural and 5-hydroxymethylfurfural, have a high potential to be used as clean and sustainable energy resources. As a consequence, investigating the properties of these compounds using molecular simulations is getting more important. In this work, we evaluated the widely used force fields OPLS-AA, GAFF and CHARMM27 regarding their ability to correctly reproduce the thermodynamic properties of these bio-masses and bio-fuels. We therefore performed Gibbs Ensemble Monte Carlo and molecular dynamics simulations to compute liquid density and vapor liquid equilibria properties of these bio-compounds. Subsequently, the obtained results are compared with experimental data and correlations by the PC-SAFT equation of state. Notably, the reported data create a basis for selecting an appropriate model or modifying the available force fields for further investigations.

Global Optimization of the Lennard-Jones Parameters for the Drude Polarizable Force Field.

Molecular dynamics (MD) simulations based on atomic models play an important role in the drug-discovery process to screen molecules, estimate binding free energies, and optimize lead compounds in chemical space. Accurate computations of thermodynamic and kinetic properties using MD simulations are highly dependent on the accuracy of the underlying atomic force field. In this context, going beyond the nonpolarizable fixed-charge model by accounting explicitly for induced polarization is highly desirable. The CHARMM polarizable force field based on classical Drude oscillators, in which an auxiliary charged particle is attached via a harmonic spring to its parent nucleus, offers both a computationally convenient and rigorous framework to model explicitly induced electronic polarization in MD simulations. For any molecule of interest, electrostatic partial charges, atomic polarizabilities, and Thole shielding factors, as well as bonded parameters can either be determined from ab initio calculations or ascribed from the knowledge-based library of the CHARMM Generalized force field. While this approach is fairly reliable in general, it is well understood that the overall accuracy of the models with respect to thermodynamic properties such as bulk density, enthalpies, and solvation free energies is particularly sensitive to the nonbonded Lennard-Jones (LJ) parameters. In the present study, we systematically refined the set of LJ parameters for the atom types available in the Drude force field to best match the experimental thermodynamic properties for 416 small druglike organic molecules. To further test the transferability of the optimized parameters, the hydration free energy of 372 molecules was computed. The calculations resulted in a small average error of 0.46 kcal/mol and a Pearson R of 0.9, representing a significant improvement over the additive GAFF force field in our previous study, where an average error of ∼2 kcal/mol was obtained. Such an improvement is consistent with the ability of the polarizable Drude model to more accurately model interactions in different environments. The effort provides a roadmap for the global optimization of force field parameters using experimental data. It is hoped that the present effort will further the application of the Drude polarizable force field in molecular simulations including drug design and discovery.

Evaluation of thermal conductivity of organic phase-change materials from equilibrium and non-equilibrium computer simulations: Paraffin as a test case

An accurate in silico evaluation of the thermal conductivity is critical for improving the thermal properties of organic phase-change materials on a rational basis. To explore the impact of a theoretical model on the computed thermal conductivity, here we employed the equilibrium and the non-equilibrium molecular dynamics (MD) simulations to study paraffin (n-eicosane) bulk samples, in both crystalline and liquid states, with the use of 10 atomistic force fields, both all-atom and united-atom ones. Overall, we found that the equilibrium MD method is preferable for computing the thermal conductivity of n-eicosane samples (at least for a 10-nm-size simulation box). For the n-eicosane crystals, the all-atom models provide larger thermal conductivity coefficients than their united-atom counterparts and, correspondingly, a better match with the experimental data. This is most likely because the crystalline lattice of the models with explicit hydrogen atoms is additionally stabilized by the electrostatic interactions. In contrast, in the liquid state, most all-atom models overestimate the experimental data for n-eicosane, providing thereby worse performance as compared to the united-atom force fields. However, when it comes to the experimentally observed increase in the thermal conductivity of n-eicosane samples upon crystallization, only all-atom models are able to reproduce quantitatively the experimental data. Each force field of n-eicosane was also characterized by an overall score which accumulated the deviations of the computed thermal conductivity coefficients from the experimental values, for both crystalline and liquid samples. It turns out that the best performance among 10 atomistic models of n-eicosane is observed for the all-atom GAFF force field. All in all, our study clearly demonstrates that a proper choice of the model for computing the thermal conductivity is a non-trivial task: even for such relatively simple compounds as paraffins (n-alkanes), different models perform quite differently, in equilibrium and in non-equilibrium MD simulations, as well as in crystalline and liquid phases.

Computation of the Dissociation Temperature of TBAB Semiclathrate in an Aqueous Solution Using Molecular Simulations.

Quaternary ammonium salts such as tetra-n-butyl ammonium bromide (TBAB) are known to form semiclathrate hydrates. Since they form at much milder conditions compared to gas hydrates, they have evoked much interest in development of new technologies for gas storage and gas separations. In this work, we present a method to compute the phase equilibrium of TBAB semiclathrate. The TBAB molecule is modeled using OPLS-AA and GAFF force fields and the results from simulations are compared with experimental data to determine the ability of the force fields to accurately predict the semiclathrate phase equilibria. It was observed that the OPLS-AA force field outperforms the GAFF force field in predicting the experimental phase equilibrium data.

Computation of gas solubilities in choline chloride urea and choline chloride ethylene glycol deep eutectic solvents using Monte Carlo simulations

Deep eutectic solvents (DESs) are considered as green alternatives to room temperature ionic liquids (RTILs), due to their lower-cost synthesis and more environmentally friendly nature. In this work, Monte Carlo (MC) simulations have been used to compute the solubilities of CO2, H2S, CH4, CO, H2, and N2 in choline chloride urea (ChClU) and choline chloride ethylene glycol (ChClEg) DESs. Due to the strong intermolecular interactions of DESs, leading to high viscosities, MC simulations present significant challenges with respect to system equilibration and solute molecule insertions. The Continuous Fractional Component Monte Carlo (CFCMC) method has been used with our open-source code, Brick-CFCMC, to improve molecule insertions and equilibration of the system, and directly compute the excess chemical potential and solubility (in terms of the Henry constant) of the gas molecules in the DESs. Pure DES properties, such as density and radial distribution functions (RDFs), were well reproduced by MC simulations. The solubilities of gases were, however, underestimated by the CFCMC simulations compared to available experimental data from literature. The order of solubilities of the different gases in ChClU at 328 K was obtained as H2S > CO2 > CH4 > H2 > CO > N2, which reasonably agrees with experimental data from literature. The OPLS force field resulted in larger average Henry constants in ChClEg, compared to the GAFF force field, implying the better suitability of the GAFF force field for the calculations. Smaller ionic charge scaling factors were shown to increase the solubilities of the gases in the DESs, but result in lower densities. The differences between the computed Henry constants from MC simulations and experimental data from literature may be caused by the unsuitability of the used force field parameters of the DESs in combination with those of the solute gases. Nonetheless, experimental data from literature is scarce (except for CO2) and in some cases contradictory, which makes the comparison with the computational results difficult.

A theoretical investigation of different point charges combined with GAFF and OPLS-AA for acetic anhydride

Four different charge methods such as RESP2, AM1-BCC, CM1 and CM5 were employed to investigate the properties of acetic anhydride in a bulk phase simulation with GAFF and OPLS force field potentials. These charges were obtained from quantum mechanical calculations by fitting electrostatic potential for RESP2, AM1-BCC and Hirshfeld population for CM5 charges. The calculated density and viscosity values from molecular dynamics simulation were in good agreement with experimental values; however, GAFF with RESP2/AM1-BCC models exhibited higher error (2.1–3.07%) than OPLS with CM1/CM5, which showed the error around 0.1–0.5%.

SAMPL6 Octanol-water partition coefficients from alchemical free energy calculations with MBIS atomic charges.

In molecular modeling the description of the interactions between molecules forms the basis for a correct prediction of macroscopic observables. Here, we derive atomic charges from the implicitly polarized electron density of 11 molecules in the SAMPL6 challenge using the Hirshfeld-I and Minimal Basis Set Iterative Stockholder (MBIS) partitioning method. These atomic charges combined with other parameters in the GAFF force field and different water/octanol models were then used in alchemical free energy calculations to obtain hydration and solvation free energies, which after correction for the polarization cost, result in the blind prediction of the partition coefficient. From the tested partitioning methods and water models the S-MBIS atomic charges with the TIP3P water model presented the smallest deviation from the experiment. Conformational dependence of the free energies and the energetic cost associated with the polarization of the electron density are discussed.

Cooperativity of hydrogen bonding network in microsolvated biotin, the ligand of avidin class proteins.

Biotin is well known to be bound with exceptional strength by the avidin class of proteins. This ability comes from a match between the biotin-binding pocket of the protein and the structural elements of biotin, including its ureido and thiolane rings. Here we investigate the solvation shell of biotin in water as revealed by classical force field molecular dynamics with GAFF force field. Snapshots from the classical molecular dynamics were then used to generate microsolvated structures. Details of hydrogen bonding patterns present in these microsolvated structures were studied by symmetry-adapted perturbation theory (SAPT). Interaction energy values for small models of biotin hydrated by 5 or 6 water molecules show that the cooperativity constitutes 15-22% of the total interaction energy and corresponds roughly to formation of one additional hydrogen bond to biotin. The SAPT analysis shows the differences underlying hydrogen bonds of similar strength (with oxygen or sulfur atoms as the hydrogen bond acceptors, and with nitrogen atom playing a dual role of the donor and acceptor).

Python implementation of the restrained electrostatic potential charge model

AbstractThe restrained electrostatic potential (RESP) charge model is widely used in molecular dynamics simulations, especially for the AMBER and GAFF force fields. We have implemented the RESP scheme using the accessible and widely used Python language and the NumPy numerical library. This article provides a programming‐oriented introduction to the RESP scheme and highlights some of the features of NumPy that are useful in scientific computing.