Christopher Bayly

Partial atomic charges are a fundamental component underlying classical molecular simulations, but assigning charges remains a computational bottleneck; many common methods rely on quantum mechanical calculations that scale poorly with molecular size and are sensitive to the choice of conformer geometry. We introduce Open Force Field (OpenFF) AshGC, a new graph convolutional neural network charge model, as well as the Sage 2.3.0 small molecule force field for drug-like molecules parameterized to be consistent with AshGC. AshGC is designed to efficiently produce conformer- independent charges of semi-empirical quality at linear cost for molecules of all sizes, from small molecules to macromolecules. AshGC largely generates charges within the range of other accepted AM1-BCC backends such as OpenEye’s oequacpac and AmberTools’ sqm, deviating most in smaller molecules between 4-9 heavy atoms, in negatively charged molecules, and areas of chemistry underrepresented in the training set, such as particular sulfur- and phosphorus-containing functional groups. We further present the development and performance of Sage 2.3.0, which has both Lennard-Jones and valence parameters that are re-trained to be consistent with neural network charges for the first time. Benchmarks spanning gas-phase geometry optimization through protein–ligand binding free energies show Sage 2.3.0 performs comparably to earlier Sage releases, with modest improvements in condensed-phase properties and a slight decrease in non-aqueous solvation free energy accuracy. As with other OpenFF force fields, Sage 2.3.0 was validated in protein-ligand benchmarks to be compatible with Amber protein force fields. All data are publicly available, along with scripts and environments for reproducing the training and benchmarking of AshGC and Sage 2.3.0 at https://github.com/openforcefield/ashgc-v1.0-fit and github.com/openforcefield/ash-sage-rc2 respectively.

Methods for calculating the relative binding free energy (RBFE) between ligands to a target protein are gaining importance in the structure-based drug discovery domain, especially as methodological advances and automation improve accuracy and ease of use. In an RBFE calculation, the difference between the binding affinities of two ligands to a protein is calculated by transforming one ligand into another, in the protein-ligand complex, and in solvent. Alchemical binding free energy calculations are often used for such ligand transformations. Such calculations are not without challenges, however; for example, it can be challenging to handle interfacial waters when these play a crucial role in mediating protein-ligand binding. In some cases, the exchange of the interfacial waters with solvent water might be very infrequent in the course of typical molecular simulations, and such interfacial waters can be considered trapped on the simulation timescale. In these cases, RBFE calculation between two ligands, where one ligand binds with a trapped water while the other ligand displaces it, can result in inaccuracies if the surrounding water structure is not sampled adequately for both ligands. So far, a popular choice for treating the trapped waters in RBFE calculations is to combine free energy calculations with enhanced sampling methods that insert/delete waters in the binding site. Despite recent developments in the enhanced sampling methods, they can result in hysteresis in the RBFE estimate, depending on whether the simulations were started with or without the trapped waters. In this study, we introduce an alternative method, separation of states, to calculate the RBFE between ligand pairs where the ligands bind to the protein with different numbers/positions of trapped waters. The separation of states approach treats the sampling of the trapped waters separately from the free energy calculation of the ligand transformation. In our method, a trapped water in protein's binding site is decoupled from the system first, and the cavity created by its decoupling is stabilized. We then grow a larger ligand into this cavity-- a ligand that is known to displace the trapped water. In this study, we show that our method results in precise and accurate estimates of RBFEs for ligand pairs involving the rearrangement of trapped water via RBFE calculations for five such ligand pairs. We have optimized our simulation protocol to be suited for large distributed computational resources and have automated our RBFE calculation workflow.

A wide range of density functional methods and basis sets are available to derive the electronic structure and properties of molecules. Quantum mechanical calculations are too computationally intensive for routine simulation of molecules in the condensed phase, prompting the development of computationally efficient force fields based on quantum mechanical data. Parametrizing general force fields, which cover a vast chemical space, necessitates generating sizable quantum mechanical datasets with optimized geometries and torsion scans. To achieve this efficiently, it is crucial to choose a quantum mechanical method that balances computational cost and accuracy. In this study we seek to assess the accuracy of quantum mechanical theory for specific properties such as conformer energies, and torsion energetics. To comprehensively evaluate various methods, we focus on a representative set of 59 diverse small molecules, comparing approximately 25 combinations of functional and basis sets against the reference level coupled cluster calculations at complete basis set limit.

We present an efficient polarizable electrostatic model, utilizing typed, atom-centered, polarizabilities and the fast direct approximation, designed for efficient use in molecular dynamics (MD) simulations. The model provides two convenient approaches to assigning partial charges in the context of the atomic polarizabilities. One is a generalization of RESP, called RESP-dPol, and the other, AM1-BCC-dPol, is an adaptation of the widely used AM1-BCC method. Both are designed to accurately replicate gas-phase QM electrostatic potentials. Benchmarks of this polarizable electrostatic model against gas-phase dipole moments, molecular polarizabilities, bulk liquid densities, and static dielectric constants of organic liquids, show good agreement with the reference values. Of note, the model yields markedly more accurate dielectric constants of organic liquids, relative to a matched non-polarizable force field. MD simulations with this method, which is currently parameterized for molecules containing elements C, N, O, and H, run about only 3.6-fold slower than fixed charge force fields, while simulations with the self-consistent mutual polarization average 4.5-fold slower. Our results suggest that RESP-dPol and AM1-BCC-dPol afford improved accuracy, relative to fixed charge force fields, and are good starting points for developing general, affordable, and transferable polarizable force fields. The software implementing these approaches has been designed to utilize the force field fitting frameworks developed and maintained by Open Force Field Initiative, setting the stage for further exploration of this approach to polarizable force field development.

We introduce the Open Force Field (OpenFF) 2.0.0 small molecule force field for drug-like molecules, code-named Sage, which builds upon our previous iteration, Parsley. OpenFF force fields are based on direct chemical perception, which generalizes easily to highly diverse sets of chemistries based on substructure queries. Like the previous OpenFF iterations, the Sage generation of OpenFF force fields was validated in protein-ligand simulations to be compatible with AMBER biopolymer force fields. In this paper we detail the methodology used to develop this force field, as well as the innovations and improvements introduced since the release of Parsley 1.0.0. One particularly significant feature of Sage is a set of improved Lennard-Jones (LJ) parameters retrained against condensed phase mixture data, the first refit of LJ parameters in the OpenFF small molecule force field line. Sage also includes valence parameters refit to a larger database of quantum chemical calculations than previous versions, as well as improvements in how this fitting is performed. Force field benchmarks show improvements in general metrics of performance against quantum chemistry reference data such as root mean square deviations (RMSD) of optimized conformer geometries, torsion fingerprint deviations (TFD), and improved relative conformer energetics (ΔΔ𝐸). We present a variety of benchmarks for these metrics against our previous force fields as well as in some cases other small molecule biomolecular force fields. Sage also demonstrates improved performance in estimating physical properties, including comparison against experimental data from various thermodynamic databases for small molecule properties such as Δ𝐻_𝑚𝑖𝑥, ρ(𝑥), Δ𝐺_𝑠𝑜𝑙𝑣 and Δ𝐺_𝑡𝑟𝑎𝑛𝑠. Additionally, we benchmarked against protein-ligand binding free energies (Δ𝐺_𝑏𝑖𝑛𝑑), where Sage yields results statistically similar to previous force fields. All the data is made publicly available along with complete details on how to reproduce the training results at https://github.com/openforcefield/openff-sage.

Free energy calculations are rapidly becoming indispensable in structure-enabled drug discovery programs. As new methods, force fields, and implementations are developed, assessing their accuracy on real-world systems (benchmarking) becomes critical to provide users with an assessment of the accuracy expected when these methods are applied within their domain of applicability, and developers with a way to assess the expected impact of new methodologies. These assessments require construction of a benchmark—a set of well-prepared, high quality systems with corresponding experimental measurements designed to ensure the resulting calculations provide a realistic assessment of expected performance. To date, the community has not yet adopted a common standardized benchmark, and existing benchmark reports suffer from a myriad of issues, including poor data quality, limited statistical power, and statistically deficient analyses, all of which can conspire to produce benchmarks that are poorly predictive of real-world performance. Here, we address these issues by presenting guidelines for (1) curating experimental data to develop meaningful benchmark sets, (2) preparing benchmark inputs according to best practices to facilitate widespread adoption, and (3) analysis of the resulting predictions to enable statistically meaningful comparisons among methods and force fields. We highlight challenges and open questions that remain to be solved in these areas, as well as recommendations for the collection of new datasets that might optimally serve to measure progress as methods become systematically more reliable. Finally, we provide a curated, versioned, open, standardized benchmark set adherent to these standards (protein-ligand-benchmark) and an open source toolkit for implementing standardized best practices assessments (arsenic) for the community to use as a standardized assessment tool. While our main focus is benchmarking free energy methods based on molecular simulations, these guidelines should prove useful for assessment of the rapidly growing field of machine learning methods for affinity prediction as well.

Accurate small molecule force fields are crucial for predicting thermodynamic and kinetic properties of drug-like molecules in biomolecular systems. Torsion parameters, in particular, are essential for determining conformational distribution of molecules. However, they are usually fit to computationally expensive quantum chemical torsion scans and generalize poorly to different chemical environments. Torsion parameters should ideally capture local through-space non-bonded interactions such as 1-4 steric and electrostatics and non-local through-bond effects such as conjugation and hyperconjugation. Non-local through-bond effects are sensitive to remote substituents and are a contributing factor to torsion parameters poor transferability. Here we show that fractional bond orders such as the Wiberg Bond Order (WBO) are sensitive to remote substituents and correctly captures extent of conjugation and hyperconjugation. We show that the relationship between WBO and torsion barrier heights are linear and can therefore serve as a surrogate to QC torsion barriers, and to interpolate torsion force constants. Using this approach we can reduce the number of computationally expensive QC torsion scans needed while maintaining accurate torsion parameters. We demonstrate this approach to a set of substituted benzene rings.

Accurate molecular mechanics force fields for small molecules are essential for predicting protein-ligand binding affinities in drug discovery and understanding the biophysics of biomolecular systems. Torsion potentials derived from quantum chemical (QC) calculations are critical for determining the conformational distributions of small molecules, but are computationally expensive and scale poorly with molecular size. To reduce computational cost and avoid the complications of distal through-space intramolecular interactions, molecules are generally fragmented into smaller entities to carry out QC torsion scans. However, torsion potentials, particularly for conjugated bonds, can be strongly affected by through-bond chemistry distal to the torsion it-self. Poor fragmentation schemes have the potential to significantly disrupt electronic properties in the region around the torsion by removing important, distal chemistries, leading to poor representation of the parent molecule’s chemical environment and the resulting torsion energy profile. Here we show that a rapidly computable quantity, the fractional Wiberg bond order (WBO), is a sensitive reporter on whether the chemical environment around a torsion has been disrupted. We show that the WBO can be used as a surrogate to assess the robustness of fragmentation schemes and identify conjugated bond sets. We use this concept to construct a validation set by exhaustively fragmenting a set of druglike organic molecules and examine their corresponding WBO distributions derived from accessible conformations that can be used to evaluate fragmentation schemes. To illustrate the utility of the WBO in assessing fragmentation schemes that preserve the chemical environment, we propose a new fragmentation scheme that uses rapidly-computable AM1 WBOs, which are available essentially for free as part of standard AM1-BCC partial charge assignment. This approach can simultaneously maximize the chemical equivalency of the fragment and the substructure in the larger molecule while minimizing fragment size to accelerate QC torsion potential computation for small molecules and reducing undesired through-space steric interactions.

The restrained electrostatic potential (RESP) approach is a highly regarded and widely used method of assigning partial charges to molecules for simulations. RESP uses a quantum-mechanical method that yields fortuitous overpolarization and thereby accounts only approximately for self-polarization of molecules in the condensed phase. Here we present RESP2, a next generation of this approach, where the polarity of the charges is tuned by a parameter, δ, which scales the contributions from gas- and aqueous-phase calculations. When the complete non-bonded force field model, including Lennard-Jones parameters, is optimized to liquid properties, improved accuracy is achieved, even with this reduced set of five Lennard-Jones types. We argue that RESP2 withδ ≈ 0.6 (60% aqueous, 40% gas-phase charges) is an accurate and robust method of generating partial charges, and that a small set of Lennard-Jones types is a good starting point for a systematic re-optimization of this important non-bonded term.

Many molecular simulation force fields represent the charge distributions of molecules with atom-centered partial charges, so simulations with these force fields require that partial charges be assigned to the molecules of interest. The restrained electrostatic potential (RESP) approach is a highly regarded and widely used method of assigning partial charges to varied organic compounds. RESP uses gas-phase HF/6-31G* as the underlying quantum chemical method, intending the resulting overpolarization of molecules to approximate the self-polarization that occurs in the condensed phase setting. However, it is far from clear that this fortuitous overpolarization is optimal or consistent across all compounds. In order to reach a higher level of accuracy, we propose a next generation of this approach, termed RESP2. In RESP2, the charges are derived from higher-level quantum chemical calculations carried out for both gas and aqueous phase, the latter using a continuum solvent model. The polarity of the final charges is tuned by a mixing parameter, δ, which scales the relative contributions of the gas- and aqueous-phase charges. We find that simply substituting RESP2 charges for RESP charges in the context of regular LJ parameters does not lead to clear improvement in liquid-state densities and heats of vaporization but does improve the accuracy of observables expected to depend most strongly on the accuracy of the charge model, i.e., dielectric constants and molecular dipole moments. However, when Lennard-Jones (LJ) parameters are optimized in the context of RESP charges, based on liquid properties, significant improvement in accuracy can be achieved, even with a sharply reduced set of LJ types. We argue that RESP2 with δ≈0.6 (60% aqueous and 40% gas-phase charges) is an accurate and robust method of generating atom-centered partial charges. The present study also highlights the value of optimizing LJ parameters along with the electrostatic model and suggests that a small set of LJ types can be a good starting point for a systematic re-optimization of this important nonbonded term.