Efficient and Precise Force Field Optimization for Biomolecules Using DPA-2
Abstract
Molecular simulations are essential tools in computational chemistry, enabling the prediction and understanding of molecular interactions and thermodynamic properties of biomolecules. However, traditional force fields face significant challenges in accurately representing novel molecules and complex chemical environments due to the labor-intensive process of manually setting optimization parameters and the high computational cost of quantum mechanical calculations. To overcome these difficulties, we fine-tuned a high-accuracy DPA-2 pre-trained model and applied it to optimize force field parameters on-the-fly, significantly reducing computational costs. Our method combines this fine-tuned DPA-2 model with a node-embedding-based similarity metric, allowing seamless augmentation to new chemical species without manual intervention. We applied this process to the TYK2 inhibitor and PTP1B systems and demonstrated its effectiveness through the improvement of free energy perturbation calculation results. This advancement contributes valuable insights and tools for the computational chemistry community.
keywords:
American Chemical Society, LaTeXUnknown University] College of Chemistry and Molecular Engineering, Peking University \alsoaffiliation[DP Tech] DP Technology DP Tech] DP Technology \alsoaffiliation[AISI] AI for Science Institute, Beijing DP Tech] DP Technology DP Tech] DP Technology Unknown University] College of Chemistry and Molecular Engineering, Peking University DP Tech] DP Technology \alsoaffiliation[AISI] AI for Science Institute, Beijing DP Tech] DP Technology DP Tech] DP Technology \abbreviationsIR,NMR,UV
1 Introduction
Molecular simulations play a crucial role in computational chemistry, offering powerful tools for predicting and understanding molecular interactions and the thermodynamic properties of biomolecules. Achieving accurate results in these applications requires extensive sampling of the conformation space over long timescales, ranging from nanoseconds to milliseconds. This extensive sampling is necessary to capture the dynamic behavior of molecules and ensure that the simulations reflect realistic biological conditions.
Traditional methods like ab initio calculations, which are based on quantum mechanics, provide high accuracy but are computationally prohibitive for large systems and long timescales due to their high computational costs. Neural network potentials (NNPs)1, 2, 3, 4, 5, 6 offer similar accuracy with better efficiency but still require substantial computational resources for extensive simulations7, 8.
To balance computational efficiency and accuracy, researchers often use force fields. Molecular force fields allow for faster simulations by approximating the interactions between atoms using predefined parameters. Among the various molecular force fields developed to support these calculations, a few stand out for their widespread use and ongoing development to include a wide range of drug-like molecules: GAFF9, 10, CGenFF11, 12, OpenFF13, 14 and OPLS15, 16, 17, 18. Each force field has been carefully updated and expanded over time to capture better the diversity of drug-like molecules and their interactions with biological targets, reflecting the ongoing commitment within the scientific community to enhance the tools available for drug discovery and development.
Force field optimization can be broadly classified into two categories: optimization of intermolecular terms and intramolecular terms. Optimization of intermolecular terms deals with the interactions between molecules, including electrostatic, repulsion, and dispersion interactions, and there is ongoing progress in this area19, 20, 21, 12, 10, 16, 18. In contrast, optimization of intramolecular terms focuses on forces within a molecule, such as bond stretching, angle bending, and torsion profiles, and lacks a universal solution7.
Optimization of intramolecular terms typically involves acquiring dihedral parameters through torsion scans, requiring extensive ab initio calculations, and constructing a parameter library by manually defining numerous atom types or substructure types. This process is labor-intensive and difficult to generalize to novel molecular species.
To address these limitations, this study introduces a novel on-the-fly force field optimization approach. Firstly, we fine-tuned the DPA-2 pre-trained model 22, for precise potential surface prediction, then applied it to the torsion scan process, achieving significant computational acceleration. Secondly, we developed a node-embedding-based similarity metric to replace hand-crafted atom or substructure types. This allows for seamless augmentation to any new chemical species without manual intervention, thereby ensuring parameter consistency and reducing systematic errors.
By combining these techniques, our on-the-fly force field optimization approach addresses the key challenges of high computational expense and labor-intensive parameter library construction, paving the way for more efficient and broadly applicable force field optimization.
2 Methods
The primary objective of optimizing the intramolecular interaction terms in the force field is to ensure that the calculated energies for bond stretching, angle bending, and torsion profiles closely align with the results obtained from high-accuracy ab initio calculations. Achieving this alignment enhances the precision of molecular simulations, thereby improving the accuracy of predictions related to molecular behavior. Previous studies 13 and our experiments in Fitting in Rotamer Space to Avoid Inconsistencies Between High-Accuracy Potential Surfaces and Molecular Force Fields have shown that quantum mechanical (QM) and molecular mechanics (MM) potential energy surfaces often exhibit different local minima. Fitting in rotamer space can mitigate the impact of these inconsistencies.
The manual design of force field parameters is labor-intensive and challenging to extend to new scenarios. Recent advancements leverage methods such as graph neural networks to identify atoms with similar chemical environments, providing a more efficient and scalable solution23.
Therefore, we designed a force field tuning process, detailed in Figure 1, which involves the following steps:
-
1.
Molecule Fragmentation: Decompose complex organic molecules into fragments containing at least one rotamer using a fragmentation method.
-
2.
Flexible Scan: Perform flexible scans on the high-accuracy potential surface (QM or NNP) around each rotamer for every fragment. This involves fixing the dihedral angle of the rotatable bond and optimizing bond lengths and angles to obtain relaxed conformations.
-
3.
Parameter Optimization: Catalog each fragment’s fingerprint in a fragment library to identify the dihedral parameters that need optimization. Optimize MM dihedral parameters to minimize the relative error between high-accuracy (QM or NNP) and MM potential surfacesc. The fragments’ fingerprints are based on molecular topologies using a node-embedding-based approach, ensuring that similar local environments share identical fingerprints.
-
4.
Parameter Matching: Match the optimized parameters to complex organic molecules using fingerprints.
2.1 Step 1. Molecule Fragmentation
To facilitate the optimization of molecular force fields, we developed a comprehensive approach for the fragmentation of molecular structures. This process involves several key steps to ensure that the chemical integrity and significant features of the molecules are preserved during fragmentation.
Firstly, a detailed analysis is performed to identify important structural elements such as ring systems and the ortho positions relative to specific bonds. This step is crucial for understanding the chemical environment of the molecule and for guiding the subsequent fragmentation process.
To study specific torsion angles, the fragmentation method targets these angles by systematically adding neighboring atoms in layers around the torsion bond. This process encompasses neighboring atoms, functional groups, and connected ring systems, with a keen emphasis on preserving important chemical features such as specific functional groups and ring integrity. Special attention is given to the careful addition of non-carbon and non-hydrogen atoms, carbonyl groups, and specific heteroatoms, ensuring the preservation of chemical significance and functionality in the produced fragments.
Following the determination of the fragments, our method employs RDKit’s functionalities to perform the fragmentation, breaking the selected bonds and optionally capping the resulting fragment ends with methyl groups to maintain valency. This ensures that the fragments are chemically meaningful and suitable for subsequent computational analyses.
The fragmentation method is illustrated in Figure 2. This figure shows the step-by-step addition of atoms around a target torsion angle and the resulting fragments that retain important chemical features.
2.2 Step 2. Torsion Scan with Fine-tuned DPA-2-TB Model
The most time-consuming part of the outlined process is the flexible scan of the high-accuracy potential surface, which can take days using traditional QM methods. To enhance efficiency, we experimented with constructing specific NNPs for organic small molecules to replicate their potential surfaces. The speed of NNPs far exceeds that of QM methods, allowing for molecular force field fine-tuning to be completed on a laptop within minutes.
There are pioneering works on performing torsion scans of biaryl compounds with NNPs24, but their accuracy is not enough for general-purpose torsion scans, and is not easily extensible for building a parameter library.
Large atomic models (LAM) are revolutionizing the way NNPs are produced. The representative LAM, the DPA-222, pre-trained on a vast space of chemical structures and ab initio labels, can easily generalized to downstream NNP fine-tuning tasks with a data efficiency much higher than other NNPs.
We found that directly fitting the DPA-2 model to QM data of organic molecules will result in a relatively large error, as shown in Figure 10. Inspired by previous works of delta-learning, like OrbNet4 and QRNN-TB25, we developed the DPA-2-TB model. These methods utilize semi-empirical methods to describe charge information and long-range interactions that NNPs struggle to learn. In the DPA-2-TB model, delta learning methods are employed to correct GFN2-xTB26 with the DPA-2 model.
original DPA-2 | (1) | |||
(2) |
DPA-2-TB and DPA-2 models predict the atomic energy contribution based on the atomic numbers and the coordinates . represents the descriptor of atom , maps from the atomic numbers and coordinates to a hidden representation that remains invariant under translational, rotational, and permutational (only among atoms with the same atomic number) operations. The structure of the descriptor network and fitting network are detailed in the original reference. DPA-2-TB and DPA-2 share the same descriptor network. Only the delta-learning fitting network should be trained in the fine-tuning stage. The dataset construction and training details of the DPA-2-TB model are displayed in Conclusion.
In practical use, torsion scans are done by geometry optimizations under a single torsion constraint , then scan from . We apply torsion scans on DPA-2-TB and original MM potential energy surfaces separately, then calculate the energy differences as training data, to exclude the effects of other degrees of freedom and 1-4 pair interactions.
(3) |
Where is the energy contribution of a single torsion.
2.3 Step 3. Node-Embedding-Based Rotamer Fingerprint for Consistency of Molecular Force Field Parameters
To maintain robustness and transferability during torsion parameter fitting, structurally similar molecular fragments should share the same parameters as much as possible to minimize potential systematic errors.
During the development of OpenFF13 or other classical forcefields, an SMIRKS-based27, 28 or atom-type-based dihedral parameterization scheme was employed to address this issue, but it required a lot of human effort to design these patterns. We explored a node-embedding-based similarity approach, constructing graphs based on the molecular topologies to assign identical fingerprints to rotamers in the same local environment, thereby achieving parameter consistency and reducing systematic errors.
Given a molecular structure, it is transformed into a graph , where represents the set of vertices (atoms) and the set of edges (bonds). Each vertex is associated with a 131-dimensional feature vector , capturing atomic properties including atomic number, valency, and ring membership. The adjacency matrix , with , denotes the connectivity between atoms, where if a bond exists between atoms and , and otherwise.
The embedding process employs a randomly generated weight matrix , to map atomic number features into a continuous vector space. The iterative update rule for the embedding of atom , after layers, is defined as:
where denotes the set of neighbors of atom , is the feature vector of atom at layer , and is a non-linear activation function such as ReLU.
Torsional angles within the molecule are enumerated by analyzing potential rotational bonds. Each torsion is characterized by a tuple , representing the indices of the atoms involved. The embedding vector for a torsion is computed as the average of the embeddings of its direct and inverse orientations, .
To identify equivalent torsions based on their embeddings, we calculate the Euclidean distance between the embedding vectors of torsion and , . Torsions are considered equivalent if , where is a predefined threshold.
2.4 Software Development and Code Availability
To evaluate the NNP’s effectiveness, we developed a stable and reproducible testing process in NNP benchmark tool (https://github.com/WangXinyan940/schrodinger-nnp-benchmark). Users input an ASE calculator, allowing us to assess the calculator’s accuracy based on Schodinger’s QRNN work, comparing it with QM results, semi-empirical methods, and NNPs published by Schrodinger.
Subsequently, we implemented and open-sourced the dihedral parameter fine-tuning tool, DeePDih (https://github.com/deepmodeling/DeePDih), based on the dihedral optimization process. It allows for structure optimization using user-defined ASE calculators, and saves optimized dihedral parameters in a library, enabling parameter lookup and replacement for new molecules. Its structure optimization tool is JIT accelerated, achieving more than a fivefold performance increase.
3 Results and Discussions
3.1 Evaluation of DPA-2-TB Model
3.1.1 Accuracy of Relative Conformational Energies
We proceed to evaluate the energy prediction capabilities of models. Our objective is to confirm the models’ efficacy in accurately representing the geometries and energies of organic molecule conformers and their ionic counterparts. Our primary focus lies in assessing model transferability by testing on molecules not included in the training set, for the key application of calibrating torsion parameters in conventional force fields.
To this end, we primarily focus on exploring the relative conformational energies of flexible molecules. Instead of examining the heights of torsion barriers as torsion scans do, relative conformational energies delve into the potential energy surface’s minima. To assess the conformational accuracy of our models, we utilize a test set curated by Folmsbee and Hutchison29 then recalculated by Schrodinger25, comprising 576 neutral and 81 charged molecules, all anticipated to exhibit conformational variability. Hutchison et al. conducted a conformational search for these molecules, identifying up to ten low-energy conformers for each.
The test set includes molecules with up to 50 heavy atoms and 23 rotatable bonds, posing a significant challenge to our models’ generalizability, given the smaller size of our training data. We calculate the Mean Absolute Error (MAE) and the square of the Pearson correlation coefficient (R²) across the conformer sets. We employ the median values of these metrics to evaluate our ability to replicate reference energies and the accuracy of conformation ranking. Hutchison’s findings indicate that DFT methods often achieve median R² values above 0.8 and MAE below 0.3 kcal/mol; and as Figure 3(A) shows, the DPA-2-TB model and QRNN-TB model achieve an MAE of 0.2775 kcal/mol and 0.3555 kcal/mol, outperforming the best empirical approaches (GFN2-xTB) and Schrodinger-ANI model.
3.1.2 Accuracy of Optimized Geometries
To evaluate the precision of optimized geometries, we adopted the dataset from QRNN comprising ionic conformers of drug-like molecules. These conformers underwent geometry optimization to evaluate the accuracy of their optimized geometries against B97X-D/6-31G* benchmarks, measuring discrepancies in bond lengths, angles, torsions, and RMSDs across all Cartesian coordinates, as well as in relative energies and median R² values.
Our findings, detailed in Figure 3(B), indicate improvements in errors across bond lengths, angles, torsions, and Cartesian RMSDs from semiempirical methods to DPA-2-TB, achieving remarkably low errors. The results underscore the exceptional congruence of the delta-learned models’ geometries with DFT references, advocating for their potential in future research on conformer energy ranking.
3.1.3 Accuracy of Relative Torsion Energies
We have conducted accuracy assessments on new torsion scans for molecules with drug-like properties selected from QRNN’s benchmark set25, consisting of 388 distinct relaxed torsion scans (optimized using DFT at the reference level) of ”charged” species, and 112 distinct scans of ”neutral” molecules.
3.2 FEP Calculation Results
We demonstrate the effects of the DPA-2-TB corrected force field on free energy perturbation (FEP) calculations with a couple of examples31. We perform the MD simulations with a modified version of Gromacs 2021.232. The AMBER99SB-ILDN33 and general AMBER force field (GAFF) v2 force fields9 are adopted. Both protein-ligand complex and free ligand systems are solvated in an orthogonal box with a water buffer size of 1 nm. The cutoffs for both VDW and Coulomb interactions are 1.1 nm. PME is used to calculate long-range electrostatics. The softcore is applied on both Coulomb and Vdw interactions. Systems are just neutralized by adding Na¸ and Cl ions.
16 replicas with adaptive mesh are set. Bond, vdw, charge are coupled with and uniformly transform from A to B. In each replica, a energy minimization (EM) is first performed until the maximal force is less than 100 kJ/mol/nm. Next, a 20 ps NVT and 20 ps NPT are included for equilibration. Finally, a 2 fs time step MD simulation lasts for 5 ns, and the MBAR34 method is utilized to obtain the free energy difference.
As shown in Figure 4, the RMSE over the Tyk2 inhibitor series improves from 0.91 kcal/mol with GAFF2 to 0.50 kcal/mol with DPA-2-TB corrected GAFF2.
A substantial improvement is also seen in PTP1B, where the RMSE is 1.15 kcal/mol in GAFF2 and 0.91 kcal/mol in DPA-2-TB corrected GAFF2.
More specifically, in a representative edge from molecule 20667 to 23482, as illustrated in Figure 6, the absolute difference of FEP calculated compared to experimental , has lowered from 2.74 kcal/mol with GAFF2 to 0.40 kcal/mol with DPA-2-TB corrected GAFF2. The MM torsion profiles show that the conformation distribution of the sulfonamide group (highlighted gray) in the protein-ligand complex has been corrected from gauche to trans, avoiding the dissociation of the phenyl group to the subpocket, as shown in Figure 7, thus keeping the conformational sampling in a reasonable region.
4 Conclusion
In this study, we developed a novel approach for optimizing molecular force fields by combining the fine-tuned DPA-2 NNP with a node-embedding-based similarity metric. Traditional force fields often struggle with accurately representing novel molecules and complex chemical environments due to the labor-intensive process of manually setting optimization parameters and the high computational cost of QM calculations. Our method addresses these challenges by optimizing force field parameters on-the-fly, ensuring uniform parameter consistency across various molecular fragments, and significantly reducing computational costs.
Our approach involves several key steps to enhance both precision and efficiency. First, we fine-tuned the DPA-2 pre-trained model to predict high-accuracy potential surfaces. This model leverages delta learning techniques to correct semi-empirical QM calculations, allowing it to achieve precision comparable to QM methods but at a fraction of the computational cost. Second, we developed a node-embedding-based similarity metric to replace the traditional, labor-intensive process of manually defining atom types or substructure types. By transforming molecular structures into graph representations, where atoms and bonds are treated as nodes and edges, respectively, our method ensures that structurally similar molecular fragments share identical parameters. This approach enhances parameter consistency across different molecules and reduces systematic errors.
To demonstrate the effectiveness of our method, we applied it to the TYK2 inhibitor and PTP1B systems. The results showed significant improvements in FEP calculations, with the RMSE reduced from 0.91 kcal/mol to 0.50 kcal/mol for the TYK2 inhibitor series and from 1.15 kcal/mol to 0.91 kcal/mol for the PTP1B system. These findings highlight the potential of our approach to enhance the precision of molecular simulations, providing more accurate energy and free energy predictions compared to existing models.
Looking forward, we plan to extend this work in several directions. First, we aim to scan as many drug-like molecules as possible using our method to develop a relatively universal intramolecular force field. This effort will provide a comprehensive and broadly applicable tool for the computational chemistry community. Second, we intend to generalize this method to grid-based energy correction maps (CMAP)35, 36, which will make it possible to calculate the PES of ring flipping and dihedral angle correlations in biomacromolecules accurately. Third, we plan to conduct an independent study focused on improving intermolecular interactions, addressing specific challenges such as halogen bonds and heterocycles in drug molecules, and exploring new schemes for applying optimized force fields in various chemical scenarios.
4.1 Fitting in Rotamer Space to Avoid Inconsistencies Between High-Accuracy Potential Surfaces and Molecular Force Fields
Fig. 8 shows a comparative analysis of the energy differences on MM potential surfaces between QM-optimized structures and MM-optimized structures. The study highlights the potential risks of additional energy uncertainties (about 2 kcal/mol) when directly using QM-optimized conformations.
Meanwhile, the critical role of rotatable bonds in the conformational changes of organic molecules is emphasized, with a detailed analysis using Alanine as an example (Fig. 9). High-dimensional organic molecules can be effectively described in low-dimensional spaces through rotamers, ensuring consistency between traditional force fields and high-accuracy methods in rotamer space.
4.2 Data Set Construction
To construct an NNP suitable for organic small molecules, we specially designed the dataset for this NNP, including the coverage range of the pre-training dataset and model optimization based on a druglike dataset.
The SPICE dataset37, based on published structures of biomolecules, drug molecules, and their complexes generated with classical force fields and MD, serves as an ideal pre-training dataset, increasing the chemical space coverage of the pre-trained model. We selected a representative subset of conformations (approximately 80% of the total dataset) that are consistent with chemical intuition for pre-training. Energies and forces are labeled at B97M-D3BJ/def2-TZVPPD38, 39 level of theory.
The fine-tuning process utilized a druglike dataset, a molecular fragment database built on the ChEMBL dataset40 with 600,000 molecular conformations. All fragments are generated by the method described in section Step 1. Molecule Fragmentation. This dataset is generated by DP-GEN41’s active learning strategy for more continuity in conformation space, making it suitable for fine-tuning potential surfaces. Energies and forces are labeled at B97X-D/def2-SVP42, 39 level of theory. All datasets are publicly available on the OpenLAM project.
4.3 Training and Evaluation of DPA-2-TB Model
The pre-trained DPA-2 model is adopted from the 2024Q1 version of the OpenLAM project, trained on a broad spectrum of chemical and configurational datasets22 including the processed SPICE dataset described in section Data Set Construction.
In the fine-tuning stage, the root mean squared error (RMSE) of energies and forces during fine-tuning is shown in Figure 10(A) and Figure 10(B). DPA-2-TB fine-tuning with the difference between DFT and xTB, can reach an energy RMSE of 0.04 kcal/mol/atom, and force RMSE of 0.8 kcal/mol/Angstrom. Direct DPA2 fine-tuning with DFT energies and forces can only reach an RMSE of 0.1 kcal/mol/atom and 2.0 kcal/mol/Angstrom.
After fine-tuning, the forces predicted by the DPA-2-TB model are much more consistent with DFT calculations, shown in Figure 10(C).
References
- Zhang et al. 2018 Zhang, L.; Han, J.; Wang, H.; Car, R.; E, W. Deep Potential Molecular Dynamics: A Scalable Model with the Accuracy of Quantum Mechanics. Physical Review Letters 2018, 120, 143001.
- Schütt et al. 2017 Schütt, K. T.; Kindermans, P.-J.; Sauceda, H. E.; Chmiela, S.; Tkatchenko, A.; Müller, K.-R. SchNet: A continuous-filter convolutional neural network for modeling quantum interactions. 2017.
- Devereux et al. 2020 Devereux, C.; Smith, J. S.; Huddleston, K. K.; Barros, K.; Zubatyuk, R.; Isayev, O.; Roitberg, A. E. Extending the Applicability of the ANI Deep Learning Molecular Potential to Sulfur and Halogens. Journal of Chemical Theory and Computation 2020, 16, 4192–4202.
- Qiao et al. 2020 Qiao, Z.; Welborn, M.; Anandkumar, A.; Manby, F. R.; Miller, r., T. F. OrbNet: Deep learning for quantum chemistry using symmetry-adapted atomic-orbital features. J Chem Phys 2020, 153, 124111.
- Batzner et al. 2022 Batzner, S.; Musaelian, A.; Sun, L.; Geiger, M.; Mailoa, J. P.; Kornbluth, M.; Molinari, N.; Smidt, T. E.; Kozinsky, B. E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials. Nature communications 2022, 13, 2453–NA.
- Kovács et al. 2023 Kovács, D. P.; Moore, J. H.; Browning, N. J.; Batatia, I.; Horton, J. T.; Kapil, V.; Witt, W. C.; Magdău, I.-B.; Cole, D. J.; Csányi, G. MACE-OFF23: Transferable Machine Learning Force Fields for Organic Molecules. 2023.
- Lahey and Rowley 2020 Lahey, S. J.; Rowley, C. N. Simulating protein-ligand binding with neural network potentials. Chem Sci 2020, 11, 2362–2368.
- Sabanes Zariquiey et al. 2024 Sabanes Zariquiey, F.; Galvelis, R.; Gallicchio, E.; Chodera, J. D.; Markland, T. E.; De Fabritiis, G. Enhancing Protein-Ligand Binding Affinity Predictions Using Neural Network Potentials. J Chem Inf Model 2024, 64, 1481–1485.
- Wang et al. 2004 Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general amber force field. Journal of computational chemistry 2004, 25, 1157–1174.
- He et al. 2020 He, X.; Man, V. H.; Yang, W.; Lee, T.-S.; Wang, J. A fast and high-quality charge model for the next generation general AMBER force field. The Journal of chemical physics 2020, 153, 114502–NA.
- Vanommeslaeghe and MacKerell 2012 Vanommeslaeghe, K.; MacKerell, A. D. Automation of the CHARMM General Force Field (CGenFF) I: bond perception and atom typing. Journal of chemical information and modeling 2012, 52, 3144–3154.
- Soteras Gutierrez et al. 2016 Soteras Gutierrez, I.; Lin, F. Y.; Vanommeslaeghe, K.; Lemkul, J. A.; Armacost, K. A.; Brooks, r., C. L.; MacKerell, J., A. D. Parametrization of halogen bonds in the CHARMM general force field: Improved treatment of ligand-protein interactions. Bioorg Med Chem 2016, 24, 4812–4825.
- Qiu et al. 2021 Qiu, Y. et al. Development and Benchmarking of Open Force Field v1.0.0-the Parsley Small-Molecule Force Field. Journal of chemical theory and computation 2021, 17, 6262–6280.
- Boothroyd et al. 2023 Boothroyd, S. et al. Development and Benchmarking of Open Force Field 2.0.0: The Sage Small Molecule Force Field. J Chem Theory Comput 2023, 19, 3251–3275.
- Jorgensen et al. 1996 Jorgensen, W. L.; Maxwell, D.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. Journal of the American Chemical Society 1996, 118, 11225–11236.
- Harder et al. 2015 Harder, E. et al. OPLS3: A Force Field Providing Broad Coverage of Drug-like Small Molecules and Proteins. Journal of chemical theory and computation 2015, 12, 281–296.
- Roos et al. 2019 Roos, K.; Wu, C.; Damm, W.; Reboul, M.; Stevenson, J. M.; Lu, C.; Dahlgren, M. K.; Mondal, S.; Chen, W.; Wang, L.; Abel, R.; Friesner, R. A.; Harder, E. D. OPLS3e: Extending Force Field Coverage for Drug-Like Small Molecules. J Chem Theory Comput 2019, 15, 1863–1874.
- Lu et al. 2021 Lu, C.; Wu, C.; Ghoreishi, D.; Chen, W.; Wang, L.; Damm, W.; Ross, G. A.; Dahlgren, M. K.; Russell, E.; Von Bargen, C. D.; Abel, R.; Friesner, R. A.; Harder, E. OPLS4: Improving Force Field Accuracy on Challenging Regimes of Chemical Space. Journal of chemical theory and computation 2021, 17, 4291–4300.
- Jakalian et al. 2000 Jakalian, A.; Bush, B. L.; Jack, D. B.; Bayly, C. I. Fast, efficient generation of high-quality atomic charges. AM1-BCC model: I. Method. Journal of Computational Chemistry 2000, 21, 132–146.
- Jakalian et al. 2002 Jakalian, A.; Jack, D. B.; Bayly, C. I. Fast, efficient generation of high-quality atomic charges. AM1-BCC model: II. Parameterization and validation. Journal of Computational Chemistry 2002, 23, 1623–1641.
- Bayly et al. 1993 Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A. A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model. The Journal of Physical Chemistry 1993, 97, 10269–10280.
- Zhang et al. 2023 Zhang, D. et al. DPA-2: Towards a universal large atomic model for molecular and material simulation. 2023.
- Wang et al. 2022 Wang, Y.; Fass, J.; Kaminow, B.; Herr, J. E.; Rufa, D.; Zhang, I.; Pulido, I.; Henry, M.; Bruce Macdonald, H. E.; Takaba, K.; Chodera, J. D. End-to-end differentiable construction of molecular mechanics force fields. Chem Sci 2022, 13, 12016–12033.
- Lahey et al. 2020 Lahey, S. J.; Thien Phuc, T. N.; Rowley, C. N. Benchmarking Force Field and the ANI Neural Network Potentials for the Torsional Potential Energy Surface of Biaryl Drug Fragments. J Chem Inf Model 2020, 60, 6258–6268.
- Jacobson et al. 2022 Jacobson, L. D.; Stevenson, J. M.; Ramezanghorbani, F.; Ghoreishi, D.; Leswing, K.; Harder, E. D.; Abel, R. Transferable Neural Network Potential Energy Surfaces for Closed-Shell Organic Molecules: Extension to Ions. J Chem Theory Comput 2022, 18, 2354–2366.
- Bannwarth et al. 2019 Bannwarth, C.; Ehlert, S.; Grimme, S. GFN2-xTB-An Accurate and Broadly Parametrized Self-Consistent Tight-Binding Quantum Chemical Method with Multipole Electrostatics and Density-Dependent Dispersion Contributions. J Chem Theory Comput 2019, 15, 1652–1671.
- 27 SMIRKS. https://www.daylight.com/dayhtml/doc/theory/theory.smirks.html.
- Mobley et al. 2018 Mobley, D. L.; Bannan, C. C.; Rizzi, A.; Bayly, C. I.; Chodera, J. D.; Lim, V. T.; Lim, N. M.; Beauchamp, K. A.; Slochower, D. R.; Shirts, M. R.; Gilson, M. K.; Eastman, P. K. Escaping Atom Types in Force Fields Using Direct Chemical Perception. Journal of Chemical Theory and Computation 2018, 14, 6076–6092, PMID: 30351006.
- Folmsbee and Hutchison 2020 Folmsbee, D.; Hutchison, G. Assessing conformer energies using electronic structure and machine learning methods. International Journal of Quantum Chemistry 2020, 121.
- Stevenson et al. 2019 Stevenson, J.; Jacobson, L. D.; Zhao, Y.; Wu, C.; Maple, J.; Leswing, K.; Harder, E.; Abel, R. Schrodinger-ANI: An Eight-Element Neural Network Interaction Potential with Greatly Expanded Coverage of Druglike Chemical Space. 2019,
- Wang et al. 2015 Wang, L. et al. Accurate and Reliable Prediction of Relative Ligand Binding Potency in Prospective Drug Discovery by Way of a Modern Free-Energy Calculation Protocol and Force Field. Journal of the American Chemical Society 2015, 137, 2695–2703.
- Abraham et al. 2015 Abraham, M.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015, 1, 19–25.
- Lindorff-Larsen et al. 2010 Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved side-chain torsion potentials for the Amber ff99SB protein force field. Proteins: Structure, Function, and Bioinformatics 2010, 78, 1950–1958.
- Shirts and Chodera 2008 Shirts, M. R.; Chodera, J. D. Statistically optimal analysis of samples from multiple equilibrium states. The Journal of Chemical Physics 2008, 129.
- MacKerell Jr et al. 2004 MacKerell Jr, A. D.; Feig, M.; Brooks, C. L. Improved treatment of the protein backbone in empirical force fields. Journal of the American Chemical Society 2004, 126, 698–699.
- Mackerell Jr. et al. 2004 Mackerell Jr., A. D.; Feig, M.; Brooks III, C. L. Extending the treatment of backbone energetics in protein force fields: Limitations of gas-phase quantum mechanics in reproducing protein conformational distributions in molecular dynamics simulations. Journal of Computational Chemistry 2004, 25, 1400–1415.
- Eastman et al. 2023 Eastman, P.; Behara, P. K.; Dotson, D. L.; Galvelis, R.; Herr, J. E.; Horton, J. T.; Mao, Y.; Chodera, J. D.; Pritchard, B. P.; Wang, Y.; De Fabritiis, G.; Markland, T. E. SPICE, A Dataset of Drug-like Molecules and Peptides for Training Machine Learning Potentials. Sci Data 2023, 10, 11.
- Mardirossian and Head-Gordon 2016 Mardirossian, N.; Head-Gordon, M. omegaB97M-V: A combinatorially optimized, range-separated hybrid, meta-GGA density functional with VV10 nonlocal correlation. J Chem Phys 2016, 144, 214110.
- Rappoport and Furche 2010 Rappoport, D.; Furche, F. Property-optimized gaussian basis sets for molecular response calculations. J Chem Phys 2010, 133, 134105.
- Gaulton et al. 2012 Gaulton, A.; Bellis, L. J.; Bento, A. P.; Chambers, J.; Davies, M.; Hersey, A.; Light, Y.; McGlinchey, S.; Michalovich, D.; Al-Lazikani, B.; Overington, J. P. ChEMBL: a large-scale bioactivity database for drug discovery. Nucleic Acids Res 2012, 40, D1100–7.
- Zhang et al. 2020 Zhang, Y.; Wang, H.; Chen, W.; Zeng, J.; Zhang, L.; Wang, H.; E, W. DP-GEN: A concurrent learning platform for the generation of reliable deep learning based potential energy models. Computer Physics Communications 2020, 253, 107206–NA.
- Chai and Head-Gordon 2008 Chai, J. D.; Head-Gordon, M. Long-range corrected hybrid density functionals with damped atom-atom dispersion corrections. Phys Chem Chem Phys 2008, 10, 6615–20.