Simon Boothroyd

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 parametrized to be consistent with AshGC. AshGC is designed to efficiently produce conformer-independent charges of semiempirical 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 and 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 retrained 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 nonaqueous 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 https://github.com/openforcefield/ash-sage-rc2 respectively.

Common methods for assigning atom-centered partial charges in computational chemistry, such as RESP and AM1-BCC, rely on quantum mechanical or semiempirical calculations of the molecule of interest, which are expensive to compute and dependent on the choice of input molecular conformer(s). Graph neural network (GNN) based continuous atom embeddings have been shown to be a fast and flexible solution for partial charge assignment, but those developed so far for condensed phase modeling have usually been trained to reproduce AM1-BCC charges, which themselves seek to reproduce the HF/6-31G(d) molecular electrostatic potential. Here, we investigate the suitability of various common charge assignment schemes, including ESP and atoms-in-molecule (AIM) based approaches, as training targets for new GNN-based charge models. We show that the strengths of both approaches can be combined by cotraining GNN models to AIM charges and molecular dipoles and electrostatic potentials. We collect a data set of quantum mechanical AIM properties computed at a high level of theory (ωB97X-D/def2-tzvpp), in both vacuum and implicit solvent, and train new GNN charge models to each. Charges can be scaled between the vacuum and solvated charge sets, and combined with Lennard-Jones parameters optimized using the Open Force Field infrastructure, to yield force fields that are suitably polarized for condensed phase modeling. We further demonstrate that the charge models may be applied to explore electrostatics-driven structure-activity relationships in medicinal chemistry. The charge models are freely available at: https://github.com/cole-group/nagl-mbis/.

The Structure and TOpology Replica Molecular Mechanics (STORMM) code is a next-generation molecular simulation engine and associated libraries optimized for performance on fast, vectorized central processor units and graphics processing units (GPUs) with independent memory and tens of thousands of threads. STORMM is built to run thousands of independent molecular mechanical calculations on a single GPU with novel implementations that tune numerical precision, mathematical operations, and scarce on-chip memory resources to optimize throughput. The libraries are built around accessible classes with detailed documentation, supporting fine-grained parallelism and algorithm development as well as copying or swapping groups of systems on and off of the GPU. A primary intention of the STORMM libraries is to provide developers of atomic simulation methods with access to a high-performance molecular mechanics engine with extensive facilities to prototype and develop bespoke tools aimed toward drug discovery applications. In its present state, STORMM delivers molecular dynamics simulations of small molecules and small proteins in implicit solvent with tens to hundreds of times the throughput of conventional codes. The engineering paradigm transforms two of the most memory bandwidth-intensive aspects of condensed-phase dynamics, particle–mesh mapping, and valence interactions, into compute-bound problems for several times the scalability of existing programs. Numerical methods for compressing and streamlining the information present in stored coordinates and lookup tables are also presented, delivering improved accuracy over methods implemented in other molecular dynamics engines. The open-source code is released under the MIT license.

The Structure and TOpology Replica Molecular Mechanics (STORMM) code is a next-generation molecular simulation engine and associated libraries optimized for performance on fast, multicore central processor units (CPUs) and graphics processing units (GPUs) with independent memory and tens of thousands of threads. STORMM is built to run thousands of independent molecular mechanical calculations on a single GPU with novel implementations that optimize numerical precision, mathematical operations, throughput, and resource management. The libraries are built around accessible classes with detailed documentation, supporting fine-grained parallelism and algorithm development as well as macroscopic manipulations of groups of systems on and off of the GPU. A primary intention of the STORMM libraries is to provide developers of atomic simulation methods with access to a high-performance molecular mechanics engine with extensive facilities to prototype and develop bespoke tools aimed toward drug discovery applications. In its present state, STORMM delivers molecular dynamics simulations of small molecules and small proteins in implicit solvent with tens to hundreds of times the throughput of conventional codes. The engineering paradigm also transforms two of the most memory bandwidth-intensive aspects of condensed-phase dynamics, particle-mesh mapping and valence interactions, into compute-bound problems for several times the scalability of existing programs. Numerical methods for getting the most out of each bit of information present in stored coordinates and lookup tables are also presented, delivering improved accuracy over methods implemented in other molecular dynamics engines. The open-source code is released under the MIT license.

Software to more rapidly and accurately predict protein-ligand binding affinities is of high interest for early-stage drug discovery, and physics-based methods are among the most widely used technologies for this purpose. The accuracy of these methods depends critically on the accuracy of the potential functions that they use. Potential functions are typically trained against a combination of quantum chemical and experimental data. However, although binding affinities are among the most important quantities to predict, experimental binding affinities have not to date been integrated into the experimental data set used to train potential functions. In recent years, the use of host-guest complexes as simple and tractable models of binding thermodynamics has gained popularity due to their small size and simplicity, relative to protein-ligand systems. Host-guest complexes can also avoid ambiguities that arise in protein-ligand systems such as uncertain protonation states. Thus, experimental host-guest binding data are an appealing additional data type to integrate into the experimental data set used to optimize potential functions. Here, we report the extension of the Open Force Field Evaluator framework to enable the systematic calculation of host-guest binding free energies and their gradients with respect to force field parameters, coupled with the curation of 126 host-guest complexes with available experimental binding free energies. As an initial application of this novel infrastructure, we optimized generalized Born (GB) cavity radii for the OBC2 GB implicit solvent model against experimental data for 36 host-guest systems. This refitting led to a dramatic improvement in accuracy for both the training set and a separate test set with 90 additional host-guest systems. The optimized radii also showed encouraging transferability from host-guest systems to 59 protein-ligand systems. However, the new radii are significantly smaller than the baseline radii and lead to excessively favorable hydration free energies (HFEs). Thus, users of the OBC2 GB model currently may choose between GB cavity radii that yield more accurate binding affinities and GB cavity radii that yield more accurate HFEs. We suspect that achieving good accuracy on both will require more far-reaching adjustments to the GB model. We note that binding free-energy calculations using the OBC2 model in OpenMM gain about a 10x speedup relative to corresponding explicit solvent calculations, suggesting a future role for implicit solvent absolute binding free-energy (ABFE) calculations in virtual compound screening. This study proves the principle of using host-guest systems to train potential functions that are transferrable to protein-ligand systems and provides an infrastructure that enables a range of applications.

The Open Force Field software stack is employed to automatically train a transferable, small molecule force field, based on the double exponential functional form, on over 1000 experimental condensed phase physical properties.

Developing accurate classical force field representations of molecules is key to realizing the full potential of molecular simulations, both as a powerful route to gaining fundamental insight into a broad spectrum of chemical and biological phenomena, and for predicting physicochemical and mechanical properties of substances. The Open Force Field Consortium is an industry-funded open science effort to this end, developing open source tools for rapidly generating new, high-quality small molecule force fields. An integral aspect of this is the parameterization and assessment of force fields against high-quality, condensed-phase physical property data, curated from open data sources such the NIST ThermoML Archive, alongside quantum chemical data. The quantity of such experimental data in open data archives alone would require an onerous amount of human and compute resources to both curate and estimate manually, especially when estimations must be made for numerous sets of force field parameters. Here we present an entirely automated, highly scalable framework for evaluating physical properties and their gradients in terms of force field parameters. It is written as a modular and extensible Python framework, which employs an intelligent multiscale estimation approach that allows for the automated estimation of properties from simulation and cached simulation data, and a pluggable API for estimation of new properties. In this study we demonstrate the utility of the framework by benchmarking the OpenFF 1.0.0 small molecule force field, GAFF 1.8 and GAFF 2.1 force fields against a test set of binary density and enthalpy of mixing measurements curated using the frameworks utilities. Further, we demonstrate the framework’s utility as part of force field optimization by using it alongside ForceBalance, a framework for systematic force field optimization, to retrain a set of non-bonded van der Waals parameters against a training set of density and enthalpy of vaporization measurements.

A high level of physical detail in a molecular model improves its ability to perform high accuracy simulations but can also significantly affect its complexity and computational cost. In some situations, it is worthwhile to add complexity to a model to capture properties of interest; in others, additional complexity is unnecessary and can make simulations computationally infeasible. In this work, we demonstrate the use of Bayesian inference for molecular model selection, using Monte Carlo sampling techniques accelerated with surrogate modeling to evaluate the Bayes factor evidence for different levels of complexity in the two-centered Lennard-Jones + quadrupole (2CLJQ) fluid model. Examining three nested levels of model complexity, we demonstrate that the use of variable quadrupole and bond length parameters in this model framework is justified only for some chemistries. Through this process, we also get detailed information about the distributions and correlation of parameter values, enabling improved parametrization and parameter analysis. We also show how the choice of parameter priors, which encode previous model knowledge, can have substantial effects on the selection of models, penalizing careless introduction of additional complexity. We detail the computational techniques used in this analysis, providing a roadmap for future applications of molecular model selection via Bayesian inference and surrogate modeling.

While the solubility of a substance is a fundamental property of widespread significance, its prediction from first principles (starting from only the knowledge of the molecular structure of the solute and solvent) remains a challenge. Recently, we proposed a robust and efficient method to predict the solubility from the density of states of a solute-solvent system using classical molecular simulation. The efficiency, and indeed the generality, of the method has now been enhanced by extending it to calculate solution chemical potentials (rather than probability distributions as done previously), from which solubility may be accessed. The method has been employed to predict the chemical potential of Form 1 of urea in both water and methanol for a range of concentrations at ambient conditions and for two charge models. The chemical potential calculations were validated by thermodynamic integration with the two sets of values being in excellent agreement. The solubility determined from the chemical potentials for urea in water ranged from 0.46 to 0.50 mol kg⁻¹, while that for urea in methanol ranged from 0.62 to 0.85 mol kg⁻¹, over the temperature range 298–328 K. In common with other recent studies of solubility prediction from molecular simulation, the predicted solubilities differ markedly from experimental values, reflecting limitations of current forcefields.

Solubility is a fundamental property of widespread significance. Despite its importance, its efficient and accurate prediction from first principles remains a major challenge. Here we propose a novel method to predict the solubility of molecules using a density of states (DOS) approach from classical molecular simulation. The method offers a potential route to solubility prediction for large (including drug-like) molecules over a range of temperatures and pressures, all from a modest number of simulations. The method was employed to predict the solubility of sodium chloride in water at ambient conditions, yielding a value of 3.77(5) mol kg⁻¹. This is in close agreement with other approaches based on molecular simulation, the consensus literature value being 3.71(25) mol kg⁻¹. The predicted solubility is about half of the experimental value, the disparity being attributed to the known limitation of the Joung–Cheatham force field model employed for NaCl. The proposed method also accurately predicted the NaCl model's solubility over the temperature range 298–373 K directly from the density of states data used to predict the ambient solubility.