 So today I want to talk about the question of how many Leonard Jones atom types of force field ought to have. And to introduce this topic I'm going to go back to a slide that David Mobley used in his talk for this workshop, which is shown here, which is about automatic parameter fitting in open force field. So the basic idea here is that one starts with a force field with a set of initial parameters. And then uses the open force field infrastructure and in our case, focusing on force balance parameter optimizer to compute a set of quantities that one can compare with reference data. So these are experimentally available data like heats of vaporization and densities of pure organic liquids. As well as quantum mechanical, the available quantities like the torsional energy profiles of the rotatable bonds was so one computes these quantities using the force field. Does a comparison, and then does it an adjustment of the parameters that should yield an improved force field. So here goes back runs the calculations again compares with the reference data again and iterates to convergence at which point one has a set of presumably better force field parameters one can then benchmark the force field and if it works well, have a new force field release. But what this process doesn't address is exactly what parameters one is optimizing. You know, what are the appropriate Leonard Jones types for which one has numerical parameters to adjust what are the bond types, and so forth. This has been the topic of some work at open force field really led by Caitlin Bannon when she was a graduate student in David Mobley's lab. The process the idea of coming up with a method of automatically assigning types, but it's still, I think an unsolved problem. So we'll look at Leonard Jones Adam types in a recent force open force field release it turns out there are 20 for the elements hydrogen carbon oxygen and nitrogen. It turns out there are actually a number of types for hydrogen and fewer types for carbon, as shown in green oxygen and red and not and only one for nitrogen, and to get a sense for what's going on here. So like at hydrogen, there's a default string so all hydrogens will match this string, but if the hydrogen you're typing in your molecule happens to be bound to an SP3 carbon as shown here, then it'll be assigned to specialized parameters for that. If that carbon happens to be bound to one, two or even three electronegative elements, as shown in the next three lines then it'll be assigned further specialized parameters and so forth and you go down the list doing these assignments. Similarly for carbon there's a default parameter which could get modified if it turns out to be an SP3 type carbon. Putting this set in context so as I said there are 20 HCO and Leonard Jones types and open force field currently Gaff 1.8 has 48 as best I can tell and CGNFF which is a charm general force field recently has 130 but it's important to note that it's not really clear certainly for Gaff and I imagine for CGNFF how just really different all these parameters are in Gaff certainly many of the nitrogen parameters are identical and probably should be not be viewed as independent types even though they're listed as independent. There may not really be so many independent Leonard Jones types in these other force fields but they're listed that way and when it comes to doing an optimization like this it's not entirely clear which one should be considered as identical and kept together and which one should be allowed to vary independently. With that as background the question is how many Leonard Jones types should a force field have and I would argue that we want the minimal number of types required to reach maximum accuracy. And of course that depends on the data set you're using to assess accuracy and on what the other force field terms are but the basic ideas diagram in this little sort of conceptual graph with the y-axis shows accuracy and the x-axis shows a number of Leonard Jones types one can imagine that if you treat all elements as identical so you have one type for all elements you won't get very high accuracy. If you allow each element one type and you'll see that we're going to try that in a moment you should get better accuracy if you allow perhaps two types per element it should get better and so forth. But there's probably also some sort of asymptotic limit because at a certain point a carbon is Leonard Jones parameters probably are not very dependent on whether there's an oxygen or a sulfur five bonds away. So why do we want to minimize the number of Leonard Jones types you know why we want the minimum number to achieve maximal accuracy I think there are two reasons. Aside from just plain simplicity, one of them is to reduce the dimensionality of the optimization problem so this is a diagram of some objective function some measure measure of error. So we want to find a global optimum the epsilon and sigma as shown here with lowest error that may or may not be a simple problem but as you add more and more dimensions so more epsilons and more sigmas, it's likely to get more problematic. In addition, having fewer parameters reduces the risk of overfitting where you get good looking results on your training set but they. It doesn't really force field doesn't work well on data outside your training set. So I'm going to talk about here is a study that was led by Michael shell Pearl postdoc in my lab, where we define a set of at Leonard Jones typing models, really keeping them quite simple. For example, one of them would be each to see when so to hydrogen types, one carbon type, one oxygen type and one nitrogen type. So I think Wang's force balance gradient based optimizer to optimize the parameters associated with this model against experimental training set data, and then we evaluate the accuracy of the resulting force field against test set data. And then we compare with other models so was thrown off 99 frost a fairly recent version and gaff 1.8. The other parameters and these calculations are drawn from Smirnoff 99 frost. And then we look at a number of the Leonard Jones typing schemes we looked at, and they're, they're in several groups, the top group starts with HCO and so that's one type per element so all hydrogens are the same, all carbons are the same etc. And then we looked at splitting the hydrogen types in a chemically motivated manner, but keeping the others the same splitting the carbons splitting the oxygens or splitting the nitrogens and then taking all the splits putting them together to make a 12 type model 303 and two. The next set starts with a split between polar and a polar hydrogens, keeping the other elements the same so one type for carbon one type of oxygen, and then did similar splits for carbon, oxygen and nitrogen. And finally here we have Smirnoff 99 frost, which has 15 types represented in the training and test sets that I'll share in a moment, and gaff 1.8 which has 28 types in our content set. So I'm not going to go through this table in detail, but this is really just to give you a sense for how we did the splits. So for example, for carbon when we split carbon into three types, they were simply sp2 sp3 and sp, when we split oxygen it was simply carbonyl alcohol and ether, and so forth. And then at the bottom we have the a polar polar hydrogen splits based on whether they're bonded to a nitrogen or an oxygen or not. But we use a set of compounds with 78 molecules in it that capture 15 of the 20 Smirnoff HCO and Leonard Jones types and we did the calculations with two splits training test splits of the set of 78 compounds I'll be showing the results from only the first split today because time is short, but the results are virtually identical for the other split. So the first split is shown here the training set on the left test set on the right and the other training set. We just drew the compounds colored in orange out and use those for training and the rest for testing. For training we use heats of vaporization and densities with these pure organic liquids and for testing we use the same types of data plus dielectric constants. So let's look at the results. One simple result is that when we use the ridiculously simple model of one Leonard Jones type of element it's not that bad it's not great, but it's not as bad as one might have imagined. These are test set results. In the table we see the pictorial representation of the types here's the name of the model so we'd see when Smirnoff and Gaff. This is mean percent error on densities mean percent error on heats of vaporization mean percent error on dielectric constants and the objective function is a nice summary. It's the objective function that force balance minimizes to reduce error so a low number is better and the cells are colored where orange is worse and blue is better. So here what you see is that the objective function for HCON is 243 Smirnoff and Gaff are better. Perhaps not dramatically better if you look at the actual numbers here heats of vaporization are similar in accuracy dielectric constants are similar in accuracy maybe it's a little better but the densities for this for this model or a little bit worse. A very clear result was that simply distinguishing between polar and non-polar hydrogens gives you a huge boost in accuracy and it gives you almost the almost the entire boost in accuracy that we've seen in any other model. So you can see that in a few comparisons. Let's look at the top two lines. We go from HCON to H2CON so that polar a polar split and the objective function falls by not quite a factor of two. And these other results tend to get better so heats of vaporization get better densities get better dielectric constants that's going to be a little worse. But here are some more comparisons let's go from HC3ON and split the hydrogens we go to H2C3ON. And we get a drop in the objective function from HC3N to H2C3N we get a small drop in the objective function. And from HCON2 so split nitrogen is and now we split the hydrogens H2CON2 we get another nice reduction. And so in all cases, splitting the hydrogens into polar and a polar improves the results. And in fact gives us results that are at least on this you know test set better than Smirnoff and Gaff. I'm not trying to claim that this is a better model necessarily which is a very limited test set but it's certainly something to think about. We then can look at more generally we looked a little bit at just a very simple split for the hydrogens but if we split the hydrogens into four types we get improvement so that's 243 to 153 this first green arrow. If we split oxygen types from HCON to HC3N we also get improvement actually about the same as splitting hydrogens. Interestingly for this both for both test sets, splitting carbons makes things worse on the test set and splitting nitrogens makes things worse on the test set. Not a lot but not but it's not better. So that's losing to split seem to be of secondary value, if any value at all in this context. And finally, we find that some runs and don't find the global optimum. So here we look at a training set because that's where we're going to see whether the optimum was found or not. And I want to give you two comparisons. Let's start off with HCO3N and add parameters to make this 12 parameter model H4C303N2 objective function is worse after training. And or if we just split the hydrogens in the last line the objective function also gets worse than for the simpler model. And we know that because the more complicated models the ones with more parameters have more freedom in principle. The achievable objective function with a more complicated model has to be at least as good as that with a simpler model and probably a little bit better because there's more freedom for optimization. And we're seeing the opposite. So what happened here and actually Michael did a really nice analysis of this. What we noticed is that so HCO3N is considerably better than HCON which I don't show here. And what's happening of course is that the three different types of oxygens are being assigned distinct parameters. When he looked at H2CO3N from this optimization he noticed that the oxygen parameters were actually pretty similar to each other. And basically what had happened is we got in a polar-ray polar-hydrogen split. So that raised the possibility that somehow this optimization wandered off into a false optimum or a local optimum where the hydrogens were split but the oxygens weren't. So to check that he did a fresh optimization of H2CO3N starting from HCO3N. So we start with the split oxygens and re-optimize to get this re-optimize model and now the objective function is a little bit lower maybe the same in effect as HCO3N. So it's fallen relative to the prior optimization where we had 84. So indeed when we looked at the parameters for H2CO3N re-optimize the oxygens were again split and the hydrogens had very similar parameters. So there seem to be two pretty good ways of optimizing H2CO3N parameters. One is to split the hydrogens the other one is to split the oxygens and the optimizer risks getting trapped in the false minimal. So let me wrap up by summarizing and suggesting some implications. It looks you know from this preliminary study as though using surprisingly few atom types may be advantageous. We get competitive test set accuracy. It's a way of avoiding overfitting and it limits the search dimensionality. And indeed, even for quite few parameters so H2CO3N is what is it, you know, take the number of types and multiply by two that's so many parameters it has. There is danger of getting caught in local optima. So this suggests that in future optimization work it may be worth at least using multiple starting points in our optimizations and maybe considering other ways to improve the global search. So going back to our graph from before we can maybe put sort of an impressionistic set of data on here. This is not meant to be quantitative. But if we look at accuracy versus number of Leonard Jones types. It's not so great for HCON we get a nice boost and accuracy by adding by splitting hydrogens and we get modest boosts potentially by splitting nitrogen and oxygens. So I'll wrap up by acknowledging primarily Michael Shopperl who's really done all this work and done a beautiful job of it. Unfortunately, he'll be for us he'll be moving on to greener pastures or different shade of green pastures soon. Sophie Canton and a grad student with us who developed some of the initial ideas and leaping long who's helped with force balance and interpretation of the data. Thanks everyone in open force field especially for earlier discussions of this work and to NIH for funding. Thank you. That's great thanks. So, yeah, I guess if we were continuing. Well, I want to make sure Simon is thinking about this if we were continuing on with, you know, maybe trying to connect some of this with the new Leonard Jones refitting going into release force fields it might make sense to do a couple of slides that start with a pretty minimal model, like one of the ones you guys have tried and see how that fair is which one were ones, which would be your first picks for for trying. I would. I totally depend. Okay, so if we were only going to do one. I would do either a, I would probably do H2CO3N. And if we were going to do two or three, I would probably do the three at the top here. So H2CO when, even though I say it's lower accuracy than the other two. It's not impressively lower I mean so if we look at where is that each to see when is 139 on the objective function. That's amazing. If we add oxygen split, we get to 142 which is not as good on the test set and end to is not as good but I'm tempted to split the oxygens and nitrogen anyhow because we've seen that on their own they can improve things. What do you Michael Schapro is on the line I think he may have some thoughts on this as well. Michael Schapro you're muted if you're trying to say yes. Yeah, I'm more or less would have the same suggestions so H2CO3N would be my first choice if you only test one. And if we would test two or three. I would probably take this one and then probably the ones Mike suggested those were the ones which are standing out. You could even go further and set the bare minimum as a HCON so the just the element wise atom type or type model which would be interesting. Do you are the smart patterns you guys use for those in the paper. So, um, they wait they're here to, they're here, you see the right. Yeah, yeah, yeah. Okay. We've, I think what we did is essentially take sort of top levels, sort of parts. Yeah, smirk, whatever they are. Cool. Yeah, so I'll make sure you know Simon is thinking about this. I think it'd be neat, you know, to some extent this just depends on compute because automating it is trivial so be need to repeat whatever fitting he's doing after the 1.2 release. You know, as he finishes, as he finalizes how he's going to be fitting and testing his LJ refits to also just do the same thing with your, you know, one or two of these types of typing and see what works best. Yeah, this is a larger scale version of your same experiment. Yeah, exactly. And I think that would be fascinating.