 Good evening, everyone. I'm Bhavan Bairam. I was a open-affect postdoc in David Mobley's lab at UC Irvine. And today, I'm going to talk about small molecule force field development, and specifically about balance parameter optimization and share a few quick stories. Force field development is pretty exciting. And I think this is how every force field developer feels like after they make a new iteration of force field. And I hope our force fields only burn your computer and not your molecules. And the reason being, developing a force field, general force field broadly involves taking care of these force things, chemical typing, training data, optimization, and benchmarking. So chemical typing are the parameters you put into your force fields. Most of this is driven by chemistry wizards. And there are significant efforts by Trevor and Tobias to automate the chemical perception tests, as they have shown this morning. And next comes the training data you choose and whether it covers diverse chemistries and then the quality of the data set. And next would be the way you perform your optimization, the realization skills you use, and the initial point you're putting your optimization to and the definition of objective function, weights you are assigning to different targets. And finally, another most essential part that decides whether the force field is good or not is benchmarking the output from optimization. And if you feel something is not right, you go back and fix it and repeat the cycle. So how do we do this? Here is how OpenFF is successful in culminating all the steps involved in building a force field by developing a software and data infrastructure, bringing together all the bits and pieces needed. In this workflow, we start with an initial set of parameters and we use OpenFF toolkit for processing the molecule graphs and assigning parameters. And it has wrappers around ADKT and OpenI toolkit backends for most of the heavy chem informatics work. And then we have force balance or force field outmasters that can ingest both quantum mechanical and physical property data and outmesters the force field. After that, we perform orthogonal benchmarks on the output and make a decision on whether to make a release or further improvements are warranted. All of this infrastructure is permissively licensed so that you can use it for your commercial endeavors as well. And I think all the individual parts are industry standard and most of it is community driven. And I think this picture summarizes it like any lab that can go using our infrastructure. So today I will share some of the details of building a force field starting with our conduct chemistry data generation in the context of training and testing our force fields and many other auxiliary projects that make use of these datasets. So this is our workflow of QC data set generation where anyone can come with a set of 2D or 3D representation of molecules and make a pull request on this repository which will then process the submission and send calculations to HPC clusters at all the different universities and sends back the calculated data to QC archive which is a public repository of data hosted by Molesi. And here is the team that helped keep the infrastructure resilient over the last few years and made many submissions which resulted in a rich variety of datasets and a huge thanks to David Dodson and Brent Pritchard who handle most of the DevOps involved in it and keep it going. So here is a list of datasets that are available on QC archive that was generated by OpenFF. So for example, the OpenFF gen 1, 2, and 3 datasets curated by HESU and Simon they contain the HSCN soft mesh geometries and torsion scans of small molecules which we regularly use in our force wheel training and here is a recently curated industry benchmark set which we use to test our force wheel and we also have the OpenFF protein datasets for fitting our next generation of biopolymer and small molecule force wheels curated by Chapin and the OpenFF electrostatic potential datasets created by Simon to train and test graph neural net charge models, virtual sets and polarized force fields and David Dodson and I helped create a single point energy dataset along with bit management recently which contains more than a million conformers from the short datasets, webcam datasets and a lot of female assets and we call it the spice sets which really stress tested the capabilities of our data generation infrastructure and we can probably say that we can handle quantum chemistry data at scale and the main workhorse of our quantum calculations is Cy4 and the outmeasure we use is geometric from leaping one screw along with torsion drive for our torsion scans and QCFactile is the server package maintained by Ben Pritchard at MOLSI which makes these calculations possible and QC submit package developed by Josh Horton and OpenFF team makes the submission smooth and also smooth data retrieval within a few lines of code after the calculations are done. Here you can see like within six lines you can import the necessary packages access the remote to CR cable holes download the data and perform a few filtering tasks to pick the data you want. For example, we are filtering for molecules that match in nitrile from the optimization data set collection here. So when it comes to force field fitting one of the primary concerns is how good is our theory level for our quantum chemistry data and whether we are able to model diverse chemistries and whether the method is accurate in conformer energies accurately predicts torsion profiles and it is accurate for charged molecular prosthetics such as dipoles and accurately describes complex and diverse chemistry such as nitrogen, sulfur, et cetera and most importantly whether it is computationally efficient. You don't leaping one Daniel Smith and Victorio Lim did a benchmark of conformer energies back in 2018 and they proposed that B3DP3BJ within DCVPT basically is the optimal choice and that benchmark was on a smaller set of peptides medium sized micro cycles and canonical amino acids and some alkyl groups with different substituents and all of these were neutral molecules. So we extended the scope of this and did another benchmark recently on torsion profile energies including diverse chemistry such as hypervalence, sulfur, nitrogen in different hybridizations by alloys, halogens and a lot of charged molecules and we thought like we need to upgrade our level of PM theory but it turns out that what we had initially chosen provides a remarkably good balance of speed versus accuracy, despite its limitations and we plan on sticking with it for now. So our default B3DP3BJ within DCVP we produce the reference CCSDT in complete basis of limit energies within 0.52 kcal per mole and the top performing one is the drain separator hybrid with an error of like 0.42 kcal per mole and where was by only 0.1 kcal per mole and we don't think like there is a need to switch theory levels and dump all the data we generated over the last five years the overall gain in accuracy is not significant enough but for the subset of charged molecules the difference is a bit higher it is around 0.2 kcal per mole with respect to the best one so this can be considered depending on priorities but I have to mention here that it is super easy to add an additional compute specification and regenerate all the data sets and it is just a question of how much compute we can assign right now and what needs to be prioritized at the moment. Another major point is that CY4 does not have analytic gradients for rivet and dispersion calculations which makes these calculations really expensive as the molecule size increases. So we now have the quantum chemistry data to train the balance parameters in our force fields the next thing is how to efficiently use this data to bring parity between the MM and QM birds for this purpose we use force balance from leaping ones group to optimize our force fields and it is a really good program that solves the multidimensional problem and the reference data can be QM or fiscal property data and it supports a lot of MD engines and custom functional forms and as Josh showed this morning with the double exponential and it is also really simple to build custom objective functions as well. So when I say that I want to bring parity between the MM and QM birds by the I mean we want our MM geometries to match the QM reference and the same with conformal energies and torsion energy profiles. So in our force field fields we use optimized geometry target which minimizes the sum of residuals in internal coordinates the bond angle and torsion deviations within a molecule they are all scaled and summed up which goes into the objective function and the same with torsion profiles between QM and QM where the weighted energy differences go into the objective function. One thing to note here is that we prioritize the mesh to low energy regions for the torsion profiles. So using this infrastructure we made a few changes to say it's 2.0 and broadly they can be categorized changes pertaining to chemical typing and some related to fitting procedure that went to building stage 2.0 and this improved the accuracy of force fields and today I'm going to share a few details on the changes pertaining to hypervalent sulfur and some changes related to how we perform fitting and how choosing a good starting point really matters. So one thing we found is that careful refinement of typing is crucial for getting accurate geometries on sulfonamides otherwise we have terrible failures like the pathologies pictured here. The one on the left is a reproduction of a pathology with parsley 1.3 where the carbon sulfur nitrogen angle around sulfur is overly bent and thanks to the report by Chris Bailey and Guy Tannum we fixed that in 1.3.1 and stage 2.2 but again on the right we found few systemic issues which remain even in stage where a similar angle bending occurs in gas-based minimalist geometries when there is a non-carbon atom present beside a sulfonamide functional group and there are a lot of FTA-proof plugs that have hypervalent sulfur and it's a really prominent chemistry and we don't want to be doing bad in describing this and because of the unique electronic properties offered by sulfur which can lead to enhanced bending affinity and selectivity towards specific biological targets this is the chemistry of interest to drug hunters and almost a quarter of the hypervalent sulfur containing red sulfonamides and the initial stop of the reading fix was to avoid tangles A31 and A32 that were assigned for hypervalent sulfurs back to reasonable values manually as we observed created the parameters due to optimization and the values being shifted in one bad direction by the optimizer and this fix was in 1.3.1 and we didn't see the same pathology with sage but the other anomaly we observed with sage was that there is a subset of sulfonamides called sulfonamides and we were able to pinpoint the blame on the same angle parameters say 31 and A32 and also a few torsion terms T143 and T157 and I think the changes we made is all the issues permanently and I went to Pubchem to improve confidence in our force fields and tested molecules which have the same functional group and they are all doing fine with another iteration of force field. So another quick story I want to share with you today is on how changing the way we perform fitting can affect the force field. So we are doing a multi-objective optimization and if we are not starting from a good initial point we may end up with unphysical values for our parameters. For example, here is one of the earliest pathologies without force field where the optimist force constants were counterintuitive for carbon-carbon single bond and carbon-carbon double bond where the single bond has a higher force constant and much stiffer than the triple bond and propane molecules were exploding during MD simulations and a manual fix was made in 1.2, 1.2 to solve this issue. But in general, if we have a bunch of QM data it is really easy to derive a mean equilibrium value for the bond lengths and angles but it is difficult to derive the force constants easily and we were expecting the training to vibrational frequencies would set them right but matching the normal modes between QM and MM is not foolproof and as I mentioned, we were performing multi-objective optimization and we may end up with a paradox solution that gives a minimum of the potential landscape but may have highly unphysical values like this. So what's the solution for this? Fortunately, we can use the QM HACN data we already have and obtain the harmonic bond and angle force constants using a method called Modified Seminar View. This was shown by Danny Colant, George Wharton and we used it in training 2.1.0 and the graph on the right shows that parsley 1.2 that was trained explicitly to vibrational frequency training targets has a lower error in reproducing the vibrational frequencies and says 2.2 and 2.1 they don't use vibrational frequencies in their feeds and says use the previous force field as a starting point whereas 2.1 use the main values obtained from Modified Seminar View method on the training data and we do see significant improvement in reproducing vibrational frequencies closer to your force field that was explicitly trained to vibrational frequencies but without the pathologies that may arise from errors in frequency matching. And here in the left table, we can see that the parameters were propane pathology I showed before where the carbon-carbon triple bond is closer to its QM value in the range of 2,300 plus units and broadly the same changes have been propagated to all other chemistries and for example if you look at the force constants of single bonded aromatic double bonded and triple bonded nitrogen the parameters in sage all the optimal and produce the right gas phase geometries their auto fodder will compare to the physical intuitive values and now with 2.1 that use Modified Seminar View method as the starting point we see them all fall in line where the triple bonds are stiffer and the single bonds are softer in force constant values and we hope that this would avoid any new failure mode that we may not expect in market end mixtents and it is also an efficient use of QM data we already have and we don't need to generate any new data or build any new packages or invent any new methods. This is right out there for us to use with QP from Danny Cole's group and I think this is one of the beauties of open FF5 mind where people are contributing from their expertise and resolving systemic issues. So I talked about data the generation workflows and the quality of data we are producing and shared two quick stories about the changes we made in the recent force field and here comes the last and most crucial part of force field development, which is benchmarking. We have two sets of benchmarking one on the QM geometries in vacuum and the other being binary free energy calculations and most of the improvements we make in violence parameters show up in QM benchmarks and the global metrics we use are RMSG of MMOP mass geometries compared to their QM counterparts and then torsion fingerprint deviations and also delta, delta E like the conformer energies with respect to the QM and here we see improvements in all of these metrics globally with 2.1.0 and we can also slice and dice this data at a granular level to look into discrepancies and areas of further improvement. For example, we can look at the molecules with higher bond entangled deviations or we can also obtain parameter error metrics or better functional group strategies and many such things. As I mentioned before, our pre-release checklist also includes binding free energy benchmarks as well and luckily this time our force field is coincided with OpenFEs alchemist scale release and a huge thanks to Irfan Alibe for extending help to test the new iteration of force field. So is everything as good as it looks and what are some other measures that can help bring out any systemic discrepancies? There is one systemic discrepancy that is common to all the force fields in general discovered by Bilso from our industry partner Geninteg. Bilfoam that sometimes QM minima do not always match MM minima and sometimes higher energy QM conformers are stable on MM landscape. So for this molecule, if you look at the graph over here, on the left side in blue, I have the QM energies and in green, I have the corresponding QM energies starting from the same QM geometry and optimizing with Sage force field. And if we take the QM minima as the ground truth and expect the corresponding MM minima as geometry to be the MM minima, then we see that the ranking of energies does not hold and some of the MM geometries, they fall below the expected minima. And in purple on the right side, I have the RMSD with respect to the QM geometry for the MM openness geometry. And if we consider an RMSD cutoff of 0.6 as a threshold, then there are no matches between the QM and MM openness conformers. And we'll term them as orphans. And if we bump it to one angstrom, then we have three MM geometries that match to QM and three orphans that are unmatched MM conformers. And Bill has listed out a set of interesting metrics and for now I'm highlighting a few and see how we are doing in 2.1 when compared to Sage. I have to mention here that we didn't address this issue explicitly in 2.1, but we do see an improvement by improving some of the other areas in force field fitting. So on the left, I have a table of these metrics based on a TFD matching threshold of 0.1, where any confirm, any MM conformer that doesn't match with the QM conformer within the threshold of 0.1 is considered an orphan. And here out of total 72,000 place conformers the number of matches is pretty good at 67,000. And the number of orphans is like 5,000 with 2.1. And here we can see a drop from around 7,000 to 5,000. And the same can be observed with RMSD cutoff of one angstrom where the number of orphans drop from 10,000 to 8,000. So we are making good progress by addressing some of the other issues, but there is still a long way to go. So in conclusion, I have shown you that having benchmarked on diverse chemistries, our default B3-Lift D3-BJ within DZVP basis set is still fast and accurate enough for training and testing our force fields. And I also showed you that our QC dataset generation by plane can handle data at scale. And we resolved some chemistry specific issues related to cell phonomers, force weights and more through improved chemical typing and torsion periodicity handling in our force fields. And more for similar method provided a really good way to initialize our force field optimizations which led to greatly improved angle and bond force constants that are closer to physical intuitive values. And much more work can be done with the industry benchmark set and resolving the discrepancies. And what is the future for? Josh Wharton has shown you results with double exponential for intervals and we also have some really good preliminary results with virtual sites for addressing charges symmetry and polarizable force fields. And Rosemary, our self consistent small molecule plus biopolymer force field is on the horizon and Chapin will share more details and automated chemical perception, base marks from cover, it may improve parameter definitions significantly. And last but not the least, I want to thank David Mobley, my advisor for his support and guidance for the last three years and Chris Bailey as well. And I also want to thank other peers, Leaping Wang and Bill Swope and Peter Eastman, John Kudera and a lot of our industry partners. And I also want to thank the Opener 15, especially Simon, Jeff Wagner, David Dodson, Aisujang and Jessica Matt. And I also want to thank Mobley Lab members, especially Hannah Bowman and a special thanks to Trevor Goki with all his help. And I want to thank my funding agency and IH and all the awesome developers at Molesi and QC archive. So if you have any questions, I would take that. And I also want to mention that your feedback is highly appreciated and if you find anything wrong with our force, we'll please do raise an issue on GitHub or drop a message on Slack and we'll try to address it as much as soon as possible. Thank you all for your attention.