Custom force fields for Molecular Dynamics

Accurate force fields are essential to reliably predict the thermodynamic properties of small molecules. However, general transferable force fields like AMBER GAFF/GAFF2 and OpenFF, have relatively limited number of different torsion parameters and, therefore, might not correctly represent the conformational distribution of an exotic or less well-parameterized molecule during dynamics and FEP calculations. Computational chemists can address this by generating custom force fields in a fully automated manner using Flare™.  

This new custom parameterization workflow in Flare V6 is significantly more robust than previous versions, where the user had to manually select the torsion of interest and edit the molecule to remove unimportant atoms (a process needing good chemistry knowledge) to cut the computational cost associated with the underlying calculations. In Flare V6, all the torsions are set to be parameterized by default and the molecule is automatically fragmented around each torsion of interest without affecting the chemistry around the bond. 

Figure 1. Custom force field generation for OpenFF force fields is activated by simply ticking the box next to ‘Custom parameters’ in the Flare Dynamics Calculation set-up. The number of rotamers to be scanned can also be changed as per need.

The underlying components behind the generation of custom parameters in Flare V6 are described below.

Fragmentation

Performing rotation scan on a specific dihedral in a molecule can sometimes become complicated and time consuming if the molecule is particularly flexible and difficult to optimize. To address this issue, we fragment the molecule around the central rotatable bond based on Wiberg bond Orders (WBO). Changes in WBO are indicative of changes in the chemistry around the bond of interest, and so the parent (ligand) is fragmented to the smallest possible entity till the WBO of the bond of interest is within a specific threshold. If no such fragment can be identified then the parent molecule is used as the fragment for that specific bond. This is done by our implementation of openff-fragmenter [1] with AmberTools sqm for WBOs calculation.

Figure 2. Example TYK2 inhibitor parent molecule and its fragments as generated in our workflow. The corresponding torsion of interest in the parent and fragment are highlighted with the same color.

Torsion scan

Traditionally, custom torsion potentials are generated from quantum mechanical (QM) rotation scans, but this comes with a significant computational cost depending on the flexibility (number of rotatable bonds) of the molecule. However, we take advantage of ANI-2X [2], a machine-learned approximation to QM, that has been specifically refined to better predict torsion profiles. With this, we get the accuracy of its reference density functional theory (ωB97X/6-31G(d)) at a fractional computational cost.

Figure 3. Torsion energy profile (right) generated using ANI-2X for a fragment (left) with the bond of interest highlighted in magenta. It took less than 4 minutes to generate this profile.

Parameter optimization

Finally, the ANI generated torsion potentials are used as reference data to optimize the torsion parameters and obtain a custom force field for the molecule using ForceBalance [3-4].

The custom force fields thus generated show a significant improvement in the accuracy of binding free energy calculations (Figure 4). An added advantage with our workflow is that the results are cached, and so this calculation would only be done once for common cores in congeneric series, thus saving a significant amount of time and resources that can be used for running FEP and dynamics. Further, as the final custom parameters are stored locally, they can be shared with colleagues working on the same molecule(s) and reused in any subsequent experiment.

Figure 4. Example showing the improvement in accuracy of binding free energy calculations for a congeneric series of TYK2 inhibitors with the use of custom OpenFF force fields (right) Vs. standard GAFF2 (left). The correlation (R2) between experimental and predicted binding free energies improved from 0.4 to almost 1, while the Mean Unsigned Error (MUE) reduced from 0.83 to just 0.31 kcal/mol.

Sneak peek

The upcoming release of Flare will see a significant speed improvement over the current version along with some user options in the custom force fields workflow. For instance, the AmberTools sqm based WBO calculation in the fragmentation task has been replaced with the much faster semi-empirical program, xTB, offering a speed up by at least 4-12 times compared to the current version.

Another noteworthy improvement is in the way torsion scans are performed. In the current version, we optimize the geometry of pre-generated conformers individually at different torsion angles using ANI with geomeTRIC [5]. Though this makes the optimization tasks highly parallelizable and fast, the final torsion profiles may not look smooth and the low energy conformations are not always optimal. We, therefore, have replaced this with our own implementation of TorsionDrive [6] with wavefront propagation to get more accurate and smoother profiles. Further, the user can now choose the level of theory at which they would like to run the torsion scan: machine learned ANI-2X potentials or semi-empirical xTB method (>30x faster).

Try Flare on your project

Request an evaluation to try Flare on your project.

References

1. Stern, Chaya D., et al. 'Capturing non-local through-bond effects when fragmenting molecules for quantum chemical torsion scans.' bioRxiv, doi: 10.1101/2020.08.27.270934
2. Devereux, Christian, et al. "Extending the applicability of the ANI deep learning molecular potential to sulfur and halogens." Journal of Chemical Theory and Computation 16.7 (2020): 4192-4202
3. Wang, Lee-Ping, Jiahao Chen, and Troy Van Voorhis. 'Systematic parametrization of polarizable force fields from quantum chemistry data.' Journal of chemical theory and computation 9.1 (2013): 452-460
4. Wang, Lee-Ping, Todd J. Martinez, and Vijay S. Pande. 'Building force fields: An automatic, systematic, and reproducible approach.' The journal of physical chemistry letters 5.11 (2014): 1885-1891
5. Wang, L.-P.; Song, C.C. (2016) 'Geometry optimization made simple with translation and rotation coordinates', J. Chem, Phys. 144, 214108
6. Qiu, Yudong, et al. 'Driving torsion scans with wavefront propagation.' The Journal of Chemical Physics 152.24 (2020): 244116