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.

The formation of protein-ligand complexes involves displacement of water molecules that were previously occupying the protein's binding site. In some cases, however, some water molecules may not be displaced by the ligand's binding, and they can stabilize the complex by mediating the interactions between the ligand and the protein. A relative binding free energy (RBFE) calculation between two ligands, one of which binds to the protein with an intermediate water while the other displaces the water, can yield wrong results if the water fails to rearrange itself within the simulation timescale. Enhanced sampling methods have previously been used to address the sampling of such "trapped" waters, inserting or deleting waters in the protein's binding site during ligand transformation. While sometimes effective, the enhanced sampling methods typically require long simulation times to converge and may lead to differences in RBFE estimates (i.e., hysteresis) based on initial water placement. In this study, we present a non-equilibrium switching (NES) method to calculate RBFEs in systems with trapped waters. Our approach requires the knowledge of the positions of the trapped waters prior to performing the free energy calculation for ligand transformation and then uses this information to efficiently calculate the RBFE between the ligands. In our simulation protocol, we perform ligand transformation in the binding site of the target protein by using three consecutive NES switches. The three NES switches implement restraints, transform the ligand, and then remove the restraints. We demonstrate that our NES simulation-based method results in RBFE estimates within 1.1 kcal mol-1 of experimental RBFEs, with associated statistical errors under 0.4 kcal mol-1, for eight systems involving trapped water displacement. Our method provides a computationally inexpensive alternative for estimating RBFEs for systems involving trapped waters by leveraging distributed computational resources.

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, electrostatic properties, and torsion energetics. To comprehensively evaluate various methods, we focus on a representative set of 59 diverse small molecules, comparing approximately 24 combinations of functional and basis sets against the reference level coupled cluster calculations at complete basis set limit.

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.

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.

In this chapter, we introduce Orion, OpenEye's cloud-native molecular design platform. We demonstrate the scientific advances that we have made in several areas of the drug discovery process, from early hit finding, through lead-optimization to drug formulation. These advances are due, in part, to the seemingly infinite compute resources at Amazon Web Services (AWS) coupled with cutting-edge scientific algorithms. We show that virtual screening of up to eight billion virtual compounds (and beyond), is not just achievable, but can be done robustly and routinely. Moreover, Orion is an excellent resource for computationally intensive applications. We show how efficient relative binding free energy calculations, getting mechanistic insights from membrane permeability simulations, and predicting small molecule crystal structures can all assist in lead selection and optimization for the next generation of medicines. Orion is more than just an efficient compute engine; it is a platform, an instant data center in the cloud, with support for complex and efficient calculations, data storage, analysis, visualization, and team communication.

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 expected 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 when these methods are deployed within their domains of applicability. 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 (PLBenchmarks) 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 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.

We present a methodology for defining and optimizing a general force field for classical molecular simulations, and we describe its use to derive the Open Force Field 1.0.0 small molecule force field, code-named Parsley. Rather than traditional atom-typing, our approach builds on the SMIRKS-native Open Force Field (SMIRNOFF) parameter assignment formalism, which handles increases in the diversity and specificity of the force field definition without needlessly increasing the complexity of the specification. Parameters are optimized with the ForceBalance tool, based on reference quantum chemical data that include torsion potential energy profiles, optimized gas-phase structures, and vibrational frequencies. These quantum reference data are computed and are maintained with QCArchive, an open-source and freely available distributed computing and database software ecosystem. In this initial application of the method, we present essentially a full optimization of all valence parameters and report tests of the resulting force field against compounds and data types outside the training set. These tests show improvements in optimized geometries and conformational energetics and demonstrate that Parsley's accuracy for liquid properties is similar to that of other general force fields, as is accuracy on binding free energies. We find that this initial Parsley force field affords accuracy similar to that of other general force fields when used to calculate relative binding free energies spanning 199 protein-ligand systems. Additionally, the resulting infrastructure allows us to rapidly optimize an entire new force field with minimal human intervention.

1Computational Chemistry, Janssen Research & Development, Turnhoutseweg 30, Beerse B-2340, Belgium; 2OpenEye Scientific Software, 9 Bisbee Court, Suite D, Santa Fe, NM 87508 USA; 3Computational and Systems Biology Program, Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center, New York, NY 10065 USA; 4MSD, The Francis Crick Institute, 1 Midland Road, London, NW1 1AT, United Kingdom; 5EaStCHEM School of Chemistry, David Brewster Road, Joseph Black Building, The King’s Buildings, Edinburgh, EH9 3FJ, UK; 6Departments of Pharmaceutical Sciences and Chemistry, University of California, Irvine, CA USA; 7Computational Chemistry & Biology, Merck KGaA, Frankfurter Str. 250, 64289 Darmstadt, Germany; 8DeepCure, 131 Dartmouth St, Boston, MA 02116 USA

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.

Force fields are used in a wide variety of contexts for classical molecular simulation, including studies on protein-ligand binding, membrane permeation, and thermophysical property prediction. The quality of these studies relies on the quality of the force fields used to represent the systems. Focusing on small molecules of fewer than 50 heavy atoms, this data compares nine force fields: GAFF, GAFF2, MMFF94, MMFF94S, OPLS3e, SMIRNOFF99Frosst, and the Open Force Field Parsley, versions 1.0, 1.1, and 1.2. On a dataset comprising 22,675 molecular structures of 3,271 molecules, we analyzed force field-optimized geometries and conformer energies compared to reference quantum mechanical (QM) data. The data was created using scripts of the benchmarkff github repository. A corresponding manuscript is submitted, a preprint is available on ChemRxiv: Lim, Victoria T.; Hahn, David F.; Tresadern, Gary; Bayly, Christopher I.; Mobley, David (2020): Benchmark Assessment of Molecular Geometries and Energies from Small Molecule Force Fields. ChemRxiv. Preprint Read below or the file README.md for further information and description of the content: # README Version: 04 Nov 2020 For Python scripts that are NOT found in these directories, please check the [BenchmarkFF Github repo](https://github.com/MobleyLab/benchmarkff/tree/master/tools). ## Procedure 1. Prep OPLS3e file for analysis: standardize format by OpenEye in case of differences and convert from kJ/mol to kcal/mol. ``` cd prep python convert_extension.py -i opls3e_minimized.sd -o opls3e.sdf ``` 2. Remove mols that couldn't parameterize by ALL FFs. ``` python get_by_tag.py -i opls3e.sdf -s "SMILES QCArchive" -list trim3.txt -o trim3_full_opls3e.sdf ``` 3. Run analysis. ``` conda activate parsley # calc ddE, RMSD, and TFD distributions python compare_ffs.py -i match.in -t 'SMILES QCArchive' --plot > metrics.out # match_minima, only in 01_analysis_all and 02_analysis_all_smaller_cutoff python match_minima.py -i match.in --plot --cutoff 1.0 --readpickle # look at specific subsets, only in 01_analysis_all python color_by_moiety.py -i match.in -p metrics.pickle -s N-N.dat azetidine.dat octahydrotetracene.dat -o scatter_tfd_3_ # look at outliers,only in 01_analysis_all and 02_analysis_all_smaller_cutoff python tailed_parameters.py -i refdata_trim_overlap_full_openff_unconstrained-1.2.0.sdf -f --metric 'TFD' --cutoff 0.12 --tag "TFD to trim_overlap_full_qcarchive.sdf" --tag_smiles "SMILES QCArchive" > output_tfd.dat ``` ## Brief description of contents * High level: ``` . ├── 00_prep │ ├── convert_extension.py │ ├── opls3e_minimized.sd OPLS3e minimized structures from Schrodinger Maestro │ ├── opls3e.sdf standardized through OpenEye tools │ ├── opt_openff*.sdf OpenFF minimized conformations ├── 01_analysis_all compare all ffs (qm, GAFF(2), MMFF94(S), Smirnoff, OpenFF-X.X, OPLS3e) ├── 02_analysis_all_smaller_cutoff compare all ffs (qm, GAFF(2), MMFF94(S), Smirnoff, OpenFF-X.X, OPLS3e) with a smaller cutoff of .3 for match_minima ├── 03_analysis_latest_ffs compare only the latest versions of ffs (qm, GAFF2, MMFF94S, OpenFF-1.2, OPLS3e) ├── 04_analysis_openff_only compare only OpenFF ffs (qm, Smirnoff, OpenFF-X.X) └── README.md ``` * Inside an output directory: ``` YY_analysis_* various output files of above mentioned scripts, some are listed and described below: ├── bar*.png parameter coverage bar plots ├── ddE.dat relative energies data ├── fig_density_*.png scatter plots of ddE vs (RMSD or TFD) for each force field ├── match.in input file for compare_ffs.py ├── metrics.out output file for compare_ffs.py ├── metrics.pickle pickle file for compare_ffs.py -- you can read this into compare_ffs instead of rerunning the full analysis ├── refdata_*.sdf output SDF files with stored RMSD / TFD scores with reference to QM for each structure ├── relene_*.dat relative energies of matched conformers ├── ridge_dde.png compared energies plot ├── ridge_rmsd.svg compared rmsds plot ├── ridge_tfd.svg compared tfds plot ├── fig_scatter_*.png scatter plots of ddE vs (RMSD or TFD). these are noisy; I don't use these ├── trim3_*.sdf input SDF files for compare_ffs.py listed in match.in file ├── violin*.* violin plot showing ddE distributions ```

Force fields are used in a wide variety of contexts for classical molecular simulation, including studies on protein-ligand binding, membrane permeation, and thermophysical property prediction. The quality of these studies relies on the quality of the force fields used to represent the systems. Focusing on small molecules of fewer than 50 heavy atoms, our aim in this work is to compare nine force fields: GAFF, GAFF2, MMFF94, MMFF94S, OPLS3e, SMIRNOFF99Frosst, and the Open Force Field Parsley, versions 1.0, 1.1 and 1.2. On a dataset comprising 22,675 molecular structures of 3,271 molecules, we analyzed force field-optimized geometries and conformer energies compared these to reference quantum mechanical (QM) data. We show that while OPLS3e performs best, the latest Open Force Field Parsley release is approaching a comparable level of accuracy in reproducing QM geometries and energetics for this set of molecules. Meanwhile, the performance of established force fields such as MMFF94s and GAFF2 is generally somewhat worse. We also find that the series of recent Open Force Field versions provide significant increases in accuracy. Our molecule set and results are available for other researchers to use in testing.

Accurate hydrogen placement in molecular modeling is crucial for studying the interactions and dynamics of biomolecular systems. It is difficult to locate hydrogen atoms from many experimental structural characterization approaches, such as due to the weak scattering of x-ray radiation. Hydrogen atoms are usually added and positioned in silico when preparing experimental structures for modeling and simulation. The carboxyl functional group is a prototypical example of a functional group that requires protonation during structure preparation. To our knowledge, when in their neutral form, carboxylic acids are typically protonated in the syn conformation by default in classical molecular modeling packages, with no consideration of alternative conformations, though we are not aware of any careful examination of this topic. Here, we investigate the general belief that carboxylic acids should always be protonated in the syn conformation. We calculate and compare the relative energetic stabilities of syn and anti acetic acid using ab initio quantum mechanical calculations and atomistic molecular dynamics simulations. We show that while the syn conformation is the preferred state, the anti state may in some cases also be present under normal NPT conditions in solution.