 Yes, looks like we're good. All right. Okay, great. Thank you. Here we go. We can get started. So welcome everyone. Thanks for joining the fifth open force field workshop. Today we will focus on benchmarking open force field by Aether protein ligand binding free energy calculation and geometry optimization. Yeah, so in the second part of the workshop we will also have interactive session I will guide through that will be focused on using custom analysis to identify FF torsion parameter russia coming. So I'm Lorenzo Damore. I joined the open FF team on March this year as a postdoctoral scientist employed by Jansen. And yeah, the idea is essentially that I will be kicking things off. Then David will tell us more about benchmarking open FF with the protein ligand binding free energies. Then I'll present the current status of small molecule conformer benchmarking. We will have then time for discussion, five minute break and then the interactive session of approximately 45 minutes and again time for question and answer. So I'd like to start showing you are a kind of map representing overall the open force field framework. So we have an initial force field to start with which can be translated by the open FF toolkit to get into the infrastructure where it can be optimized. And there are essentially two ways to optimize our force field. So the main way is pulling data from quantum mechanics for bond angles and torsions. And this is done at the QC archive. Whereas the other way is using open databases containing essentially condensate phase physical properties which together with the QC archive constitute the open FF data set. Now all this data is fed into the open FF evaluator for automated selection and correction through the force balance developed by the living one group. And then the optimization essentially is done. We get force field to assess. So essentially like in every map where there is a pin to indicate your position today we will be talking about things happening at this point of the infrastructure. Benchmarking is an essential part of the whole process as it allows us to detect any problem related to force field and we can give feedback so that the entire process can be repeated, iterated and the force field can improve in performance. As you know, season one benchmarking has been remarkable and also also collaborative effort from pharma partners in geometry optimization benchmarks. Protein ligand binding free energies were not in the scope of season one but obviously many, many people are interested in using open force field in FIE calculation and would like to see open FF assessed against also this metric. So David Dunn did this benchmarking and will now today present us the result. So thanks David for joining us and sharing the result with us. Please the screen is yours. Thank you very much Lorenzo. Here as Lorenzo said I will present the protein ligand benchmarking results where we basically want to answer the question how does open force field versions and parameters perform in protein ligand binding free energy calculations and all this work was part of my postdoc which I did during this postdoc work with the open force field and also at Janssen and this will be taken over then in the coming month also by Lorenzo. Can you go on please Lorenzo? I created a benchmark set consisting out of 22 targets, 599 ligands and more than 1100 alchemical perturbations between pairs of ligands. These 22 targets are listed here and they are available on GitHub and this is now a versioned data set. I didn't create it myself. It's coming from various sources which are listed also here in the bottom right. So the Schrodinger Jax data set which everyone knows quite well then Christina Schindler's FEP benchmark set coming from work from Merck and a couple of targets are also from Janssen and from other public data sources. All this work led to a best practices publication which is also listed here where we tried to summarize the most important learnings from this project and how to build good benchmark data sets where we can benchmark protein free ligand free energy calculations. The content and also the structures might change in the future of this benchmark data set because we always want to improve the quality of the data and maybe also extend the data set by different chemistries for example or also for other purposes benchmarking for example absolute binding free energy calculations because so far it only focuses on relative binding free energy calculations and that's something where I still have questions which we can discuss maybe later. So for example what do we want to see in a future benchmark set and how can we ensure that still it's still possible to compare new calculations with older one because if you change input structures it will also change the results. Next please one more that's yeah here these are now the results of all the more than 1,100 perturbations done with the open force field 1.0.0 so parsley and here in a general overview where you see as in the radial axis with the experimental free energy difference between pair of ligands and the polar coordinates shows the difference between the experimental and calculated values using the free energy calculations with the PMX non-equilibrium workflow and the parameters. You can see here that the targets have different performances they are targets performing very well for example the CDK2 in the top then there are targets with more outliers like the f2 alpha in the bottom and also you have different experimental spreads dynamic ranges a number of data points for targets making some targets better for benchmarking and some others less less good. Next please we can also illustrate the difference between experimental and calculated value in such a histogram where the different colors you know the values results for the different targets. This is more or less a normal distribution with a standard deviation of 2.1 kcal per mole meaning yeah that's more than half of the perturbations deviate less than 1 kcal per mole from experiment. But we are interested to learn from these results where the where might be problematic parameters or what are the problems leading to these deviations and basically we have four possible origins of errors which are the setup meaning the starting structures for example or the charges or charged states and protein and ligands then the sampling and the model accuracy so the force field and finally the experimental data could also be wrong or not suitable for for comparison to calculations and to distinguish between the different origins of error we have a couple of things to look at for example we can look at the convergence of the calculations and this can be done without even knowing the experimental values so you have convergence criteria next please which makes it possible to filter these original 1100 perturbations and only look at converge ones where in this case I look at the overlap of the work distribution from the non-equilibrium workflow as well as the standard deviations over three repeats which I always run with these calculations and if it's larger than 1.5 kcal per mole I consider it not converged or sampling issues in the calculation. If we do this we can reduce the standard deviation of the remaining perturbations to 1.5 kcal per mole and increase the number the ratio of successful or low deviations in the results and it's more likely to find issues with the force field in this filtered perturbations than all the perturbations because there we might have more sampling issues which lead to the errors. Here I show the same data more or less in a different plot it's something like a rock curve where in the x-axis we have the deviation between experimental and calculated errors and the y-axis shows the ratio of perturbations which have a deviation lower than the threshold. The thick old pink line is the total dataset consisting out of all targets and we see that if you look at this target separately that there are quite some differences so a couple of targets work really well or have lower errors and other targets are more problematic and this could also have could be either preparation issues or also maybe force field issues in the specific target so and we can start also looking at the less good targets or also conclude from this that some data sets or parts of the data sets are not that suitable for this benchmarking exercise. Next please. On the left this plot I also showed just before that's the it's again all perturbations and then we have the filtered perturbations on the right where you see that in general the curves are steeper but also sometimes the order of the perturbations the order of the targets change. Next please. For example if you look at the green and blue curve so for SHP2 and HIF2α it changes between all perturbations and the filtered perturbations so for SHP2 a lot of calculations or values are filtered out leading to a steeper curve and this might be that there's just convergence issues so difficult perturbations which would need longer sampling times for example and whereas in the HIF2α the simulations in general converge but still have large deviations which might be force field issues or wrong starting poses. Now let's compare the open force field parameters to other force fields and we can as shown here is compared to three other force fields so the Zgen FF and GAF2 force fields both calculated with PMX by Vitas Gapsis and then the FEP plus calculations with the OPLS3E force fields. In general all four force fields agree with each other and provide similar performance with non-significant differences which makes it very difficult to draw general conclusions about the force fields and how for example to improve the open force field parameters so the differences are more or less in the noise of the sampling and also issues with the starting conformations potentially. One thing to note here the CMED open force field results look quite bad but that's actually a starting structure issue which I have resolved but not yet updated in this plot and as a side note it's still preliminary data here it's not the final which hopefully will be published soon. Next please we can also look at the different open force field versions again here we see only differences within the error bars between the force fields but it looks like we generally improve with the force field versions so the release candidate of the new recently released SAGE force field shows slightly better results in general for the Merck dataset in this case. Next I have what I should say to the plots before that was really the comparison the Wutmein square deviation based on the free energy relatively energy differences meaning the delta-delta Gs which is more the raw data which where we could rather see force field differences if we go to free energy differences so delta Gs there's some post-processing step which like cycle error corrections cycle closure corrections involved which might again shadow a bit force field differences. Here for example if you look at the SHP2 case you see quite some differences between the three open force field versions but now if we go to the next slide and look at the delta G correlations you actually don't see that big differences between the delta Gs for the three open force field versions at the top and Opel S3e also not that much better in the raw delta delta Gs has a very good performance when looking at the delta Gs so that's something also to consider that this post-processing step generally improves the the data or the comparison to experiment and also can improve therefore or shadows the parameter differences. Here we compare the outliers and also successes outliers at the left between the three force fields so parsley Opel S3e and GAF2 and we see that for the outliers we have many which are specific to the to the different force fields whereas for the successes on the right we have many perturbations which are similar or yeah the common successes are quite larger in general and we can start looking at these different subsets and maybe find issues then for the parameters. In the next slide Lorenzo please so here we can see for example one edge on the top left which is only an outlier in the parsley force field is in a JNK1 dataset and that might be I'm just yeah guessing it could be that there's something with the description of the bond which still needs to be improved or was also improved in the following open ff versions. Then here in the bottom left we have an edge which is only an outlier for the open force field and the GAF force field and that's something which is yeah which is a mutation in heterocycles and here we could also it could be also that just FEP plus does a better job in sampling this and therefore if it's not an outlier with Opel S3e or that's the the force field parameter description and finally if it's if they are outliers for all three force fields it could also mean that there's some issue with the experimental value or that they are set up errors which which lead to wrong results for example in this MCL one case we have a meta substituted fennel ring at bottom right so this methyl could be on both sides of the fennel ring and depending on this starting structure you will get different results and if you pick the the wrong one you might yeah you have an outlier basically in all force fields and with this uh Lorenzo please next I'm coming to the conclusions we have got encouraging open ff results for probably in ligand free energy calculations across all the versions and even sharing slight improvements with the newer versions but the force field differences are often in the noise of the calculations and there are still the other origins of errors and deviations from experiment which can also be found in this which were also found and in during this analysis and especially like the structure preparation issues further work is needed to automatize this benchmarking to to increase the reproducibility and do it more systematically with every new release and finally yeah as you have seen there are other issues with the free energy calculations which lead to to errors and so the force field development is only one part of many aspects which have to be developed in parallel and that's also illustrated on the right that we have next to the phosphate parameters we have to find good datasets with good uh high quality experimental values structures we can still do a lot of work with the simulation planning meaning which uh perturbations to simulate at all then the sampling is always an issue or there's a base to improve the sampling and finally you have the post processing and another step um where you also can still improve the raw delta delta g values and uh yeah i'm open to questions or even if you have if you have time we can have a discussion what kind of further analysis you would like to see or calculations and but what would be a good way to continue the protein ligand benchmarking in the future and how to improve the benchmarking dataset thank you so let's maybe stop for a second and see if anyone has questions off the top of their head for david before we move on yeah then let's keep going and we'll have another q and a at the end of lorenzo section yeah so well thanks david it was a really nice nice talk i guess that we will have time to for for discussion also at the end of the talk session so yeah let's now focus on small molecule geometry optimization benchmark that have been done in uh in season one so uh season one has been officially kicked off on january 22 and now it's complete i already mentioned that has been an awesome collaborative effort to which the pharma partners of the consortium also took part we had bsf buyer voringer engelheim bristol meyers genentech jansen merc rush vertex and excel pi taking part to this effort and contributing to constitute a proprietary dataset of drug molecules actually used in their projects where over 65 000 optimization have been run and then the pharma partners that are highlighted in blue also contributed constituting a public dataset including also a contribution of from william sloop sloop burning set where overall more than 70 000 optimization have been run the effort was also assisted by four different developers david dotson david dan jeff wagner joe shorton and by also for different subject other experts again joe shorton bill whooper leaping one and myself so season one included benchmarks with the gaff smirnov and parsley forcefield smirnov was first released back in 2016 and is essentially an initial port of the part 99 forcefield used but using the smirnov format parsley was first released in 2019 and it's the first incremental generation of open forcefield containing a balanced parameter retrained against a fully redesigned quantum chemical set that essentially enable a better coverage of the chemical space and better performance obviously in regions that are of interest of such chemical space with essentially pharmaceutical relevance not included in season one but released only one month ago in july i guess yeah of this year we now have sage that is the next and current incremental generation of open ff which is focused on a refitted set of undervalued parameter trained against physical properties data alongside also with many other improvements in in the in the balance parameters we have run geometry optimization on the public open ff industry data set with the smirnov parsley and obviously also sage assessing them against quantum mechanical data yeah molecular energy were assessed with the delta delta e consisting in relative mm energies minus qm relative energies and then geometries were also compared through rmsd and torsion fingerprint deviation dfd which uses gaussian we did the torsion angles between the ff optimized geometry and the reference geometry rather than the yeah the cartesian difference that we see already in the rmsd and should be less dependent than rmsd on molecular sites for structural comparison purposes we have been able to see that sage shows excellent performance if you look the histogram of these these three metrics in particular in rmsd and ntfd the mm minimized structure we have sage which is in in red and then open ff 1.2 parsley in violet so there has been a significant improvement in sage reproducing the reference qm structure whereas in terms of delta delta e we do see relatively well relatively good performance of sage because as we go from from parsley 1.2 to 1.3 there is a slight regression in performance but again from going to from 1.3 back to sage it seems that this regression has been solved bad and things were finally also improved also yeah for the open ff public sorry proprietary industry dataset we haven't collected so far enough results with sage to show them but we can see from the same type of histogram that parsley was significantly was a significant improvement over smirnov with open ff 1.2 showing a good performance in the assessment of both yeah conformational conformational energies and geometries against quantum quantum mechanics yeah meanwhile season one was running new features were also added to the workflow for instance the workflow now supports optimization with opls thanks to today we've done was been working to to roll out this this ability and to execute optimization using the shredding and minaries so essentially the open ff benchmark shredding commentary allows the optimization using either both opls for for shredding releases from 2021 on and with opls 3 with shredding releases since 2024 and and earlier also the optimization can be executed with both custom and running the ff builder and default parameters and essentially it's very easy to use it use the same input output behavior as the previous open ff benchmark optimized execute commentary so it's kind of familiar obviously it requires shredding binaries in particular the ff builder for custom parameters and macro model to run the optimization and obviously also an active license obviously yeah result obtained with the with the opls force field cannot be shared unless approval has been obtained and granted from from Schrodinger so we want to to recall that all the partners that want to share their results should yeah independently get approval before before doing doing so and yeah here to recapitulate a bit so season one is complete as next step you know from already from the last pharma partner call we aim to publish the result of our benchmark and Janssen will lead this manuscript effort but we would require from all partners to benchmark open ff 2.0 sage on their proprietary dataset to include also this result to the already available aggregated result of gaff smirnoff and parsley and yeah we recognize that benchmarking with opls for especially yeah perhaps running the ff builder to build the custom parameter might not be possible for all partners anyway we we strongly encourage doing so to collect as much opls data as possible and we would like setting up eventually our proprietary subset data subset of the main dataset to uh well i have to admit that it will be definitely great to have opls result from from all all partners also other features that were made available during season one i'm here referring to additional analysis that were proposed by bill swap and chaviers chavier lucas so bill swap suggested looking at relative ff energies between each ff optimized conformer and the ff minimum on the other end chavier lucas suggested looking at relative ff energies but this time between the ff conformer which is closer to the qi minimum by rmsd and the the ff minimum and these two proposed analysis have been added to the already available compared force field where we look to uh delta delta energies between each conformer so between yeah the relative ff energy and the qm relative energies and and match minima where we look to the ff relative energy but this time of the jate ff conformer which is closest to the qm i conformer by rmsd minus the ff reference conformer which is close closer to the qm the qm minimum so we we have this yeah additional analysis available and also while adding these two analysis in general code refactoring of the the whole analysis allows us now to to run them on different methods we can run them as a separate tasks which with with a remarkable improvement in time consumption for instance uh we have uh here we our uh folder uh containing the the result so we can ls over this this folder store the result of the ls uh essentially the all the sub the director is containing the the result with all the ff methods storing the result of ls in in this um uh variable and we are able uh essentially to um to run uh all the analysis independently as a separate task on on each method uh and for instance for the open ff industry public data set uh this analysis run in approximately 20 22 hours considering five different force field rather than uh 22 hours times five because times five because for uh it would have run uh let's say yeah uh first on the first force field the method and then on the other so as a serial let's say execution uh yeah that's that's not all from the community because we also see partners independently developing their own custom analysis of the result of the geometry optimization and i'm here referring to thomas fox that um was able to identify structural shortcoming uh bonds angles and torsion for specific uh chemical groups that were optimized with the open ff 1.3 and yeah i'm showing uh just some of the uh torsional issues that identified and also uh afterwards power um uh deeply uh look at uh to those chemistries uh and essentially show that uh in uh some cases uh open ff to uh sage it's is able to uh yeah to to solve to solve these issues uh here for instance we can see for the uh vinyl uh cio uh moieties we see here in uh in in xiang the open ff uh 1.3 uh structure that is aligned with the qm optimized geometry which is in green and we see that uh this uh torsion angle between one and zero uh as uh essentially uh deviation of uh almost 50 degrees uh with respect to the qm structure but uh the same optimized with uh with sage uh showed almost uh no no no deviation again we need uh cio uh torsion issues for this particular uh uh torsion between nine and and 10 again we see that in uh open ff 1.3 we have a deviation of uh yeah about 35 degrees with respect to the qm data uh but in in sage um this is uh reduced so the same torsion deviation reduces to uh yeah 18 15 degrees so um let's say that the deviation degree is at uh substantially and instead for this other particular bond uh again for ff 1.3 we have uh almost 35 degrees uh sorry yeah 35 and 40 uh degrees deviation whereas in in sage this values uh lowers to 20 and 24 degrees for this aryl metoxy uh moiety instead uh which is the the bonds are alighted between uh atom 8 and 9 uh yeah the deviation was about 24 25 degrees for ff 1 1.3 uh but this again has been almost completely uh solved but then we can see that the in in sage we can see that uh sage optimized structure almost overlapped perfectly uh with the qm optimized uh structure and indeed the values are uh yeah about just five degrees and another aryl metoxy uh moiety that exhibited some some issue uh this time uh we can look here either to this bond from three to two four or this other one from four to two five uh yeah identifying a torsion shortcoming of uh yeah from 30 up to 42 degrees for off f 1.3 whereas for for sage again these this values um lowers substantially to uh two 10 degrees or or even less uh final aryl metoxy moieties that were identified uh as offending for ff 1.3 uh this one highlighted here between uh atoms one and and two uh exhibited for ff 1.3 uh deviation of yeah almost 50 degrees with respect to the qm optimized structure but again lowers to just five four degrees with uh with sage uh yeah this this um uh sulfonyl uh vinyl moiety seems still to be yeah kind of problematical in in sage but the aryl metoxy moiety uh yeah almost overlaps perfectly uh in in sage and uh this last one yeah uh so the the bond highlighted here from atom zero and and two uh with the deviation of uh almost 40 degrees for ff 1.3 reduces to just 30 about 30 degrees with with sage so again we do see that um many of this uh concerning chemistry for ff 1.3 where um yeah fixed or improved substantially in in sage and yeah it was my last slide yeah i draw uh just very general conclusion so uh we have seen that benchmarking in us in is an essential part or within the open ff infrastructure uh it may allow uh the early detection of of problems and shortcoming that can be fixed early by before the ff release but can lead also to force field improvements then thanks to the continuous feedback from the community also we are pretty happy because uh in general geometry optimization uh benchmarking of small molecules has been uh fully automated we are also happy because uh as i shown sage is showing uh really promising results solving some of the issues we have seen for uh partially 1.3 regarding uh ocean of specific chemistry and yeah it probably constitute nowadays the best free open source alternative uh and we also have seen um that important inside uh can be shared by the open ff community as uh for the case of the proposed analysis or also the uh the code to identify geometrical issue that was shared by tomas fox and yeah so we finally have uh time for uh discussion now also including uh the part uh and this part if there are uh any question i see that alberto had a question for david haun okay yeah i i just saw the question in the chat by alberto so uh when running fpp with op less did you use the ff builder to build customer force field custom force field for all compounds so the answer to this is yes and these all the fvp plus calculations were not actually run by me but uh i took them from the published uh results so they were even run with uh different um showing our versions and different uh yeah uh on different uh environments so by by kristina schindler by people here at jansen and oh that's it thank you thank you but uh that makes of course um also the comparison to open s uh quite difficult because first you have there you have the custom parameters then you have a different sampling and a different workflow which might improve the results over the pmx non-equilibrium workflow or even the other way around depending on the perturbations and um also yeah it's difficult to ensure with these benchmarking that every calculation was really done with the same input structures still very impressive work both lorenza and you david thanks very much can you go to hide this engine can you go to slide number 27 about comparison the yeah yeah um have you thought of a way of comparing the shape of the curve uh no no we no we don't have any way to sorry can i ask you what what do you mean by the shape of the of the cure in the sense of that maybe there's a difference in the height of the potential yeah or maybe it's a different location so you'll be good if like if they are correlated so in fact that we do we have a way to try to compare the shape i look at at each of the angle whether there's a correlation between qm and ff potential are you doing this with torsion scans or how are you doing yeah yeah for the torsion scan yeah so yeah indeed let's say that here the conformer are separated by a torsional barrier we are essentially developing and trying to include the automated torsion scanner within the the open ff benchmarking yeah common tree so this is something that is now well David Dotson is working on on this it's not yet fully automated but we look forward to to it just to just to expand upon what lorenzo said um we've started to do some crude analysis at least for torsion energies in terms of both shape and magnitude um and so josh horton did a lot of the work looking at this but our current approach is just to kind of normalize two torsion profiles by the maximum barrier heights and then compute the rmsc and we find that actually gives a pretty reasonable metric of how well do the shapes agree so so definitely we are starting to think about how can we look at kind of correlations and especially the shapes of the potential energies but it's not directly incorporated into the benchmarking workflow just yet okay yes so far we only have looked really at single points on the potential energy surface and all the torsion benchmarking is upcoming okay i have another question for this slide do you have any plots for these different energy comparisons uh you mean uh like showing different methods uh based on these comparisons that yeah you mean uh showing the the plots also for biswup and charviers luga's analysis yes and the other one too as well right to compare force fields and match minima so all the plots that all the histogram that i showed today are done with the compare force field uh i well there was not enough time to to run all the analysis over the uh what this extended the data set uh yeah the the public and the and the probatory uh but for the sense i i i do have this kind of plots for the Janssen probatory data set uh and i i can see that yeah well it's art without showing any any any plot but uh they might answer to slightly different question regarding the the energy comparison metrics and i have one more question and that is so on the on the following slide you showed the several structures that improved between 1.2 1.3 and 2 yeah did this analysis help you improve 1.3 or did the results not go into the improvements what i'm trying to ask is is this is the improvement just purely based on the methodology of open ff 2.0 and that would be great or did this data actually go into the improvement and then maybe we are overfitting uh no i maybe someone can can answer better this question but as far as i know uh sage was not uh uh refitted or yeah was not uh don't uh based on these short coming of 1.3 yeah that's correct this was uh this was something we looked at after the release um yeah the timing happened to work out well that it happened around this around the same time i think it happened we in fact probably looked at it while we were preparing the release but so we were pleased to see that it fixed many of these issues i think there's maybe one or two that um some of this data brought up that are things we should still work on improving but these all just were fixed by the release but then improve fitting using for 2.0 and another question is uh is the code to identify these issues is that available should we all run it on our internal data sets to to try to identify some more issues uh yes so thomas fox uh very kindly shared uh his code uh with with us uh it's a jupyter notebook it's um available on on uh a slack channel which now i don't remember exactly which one but uh moreover we we will uh so during the interactive uh session we will use let's say a slightly modified version of his code to just focus on torsional uh uh deviation okay that would be great yeah uh well for the sake of um um yeah of the knowledge uh his code includes also uh part to identify uh bond angle so rather than only torsional issues uh with this code you can also identify bond angles uh and proper torsional issues and bond length any more question maybe i have a question for the whole group so you have seen now uh benchmarking of protein ligand free energy calculations and small molecule benchmarking these are no two different aspects of benchmarking force fields but one could imagine many other ways of benchmarking force fields and maybe has someone suggestions what uh you would like to see what what kind of benchmarking would be still important to be able to say in the end yeah we have a good and general force field which can be used in many applications one thing we were talking already and i think this is basically coming uh are the the torsional profiles um one could also think about conformational ensembles because these are i think the main main things i at least would use the or use uh force fields for to get and get an idea how how what are the conformations this this molecule can adopt and uh or as a as a starting point before before actually doing a quantum calculations if the force field is good enough that i don't need to do the the lengthy quantum things um this is certainly something that would be interesting for me so it comes down to a relative getting do we get the relative uh a minima uh correa or the ranking of the minima somewhat correct yeah correct sorry that data we should have what the this data set right it's just a question of doing the analysis to answer your question Thomas yeah but i think right now what what the analysis does it it throws out the relative uh energies of the of the the confirmers this is a confirmation this is an information that's not captured maybe not even is is not even in the in the jason's that are produced i'm not quite sure but i had the impression this this this information of the relative energies of the conformations is somewhat somewhat lost at least in the the way it's analyzed right now of course you're right the the data is there because you have the quantum quantum quantum energies and the and the force field energies but maybe this is uh i'm not sure i'm not sure if i totally understand what you would what you would look at in that um analysis much other than what we're already doing so maybe i need you explain it i need you to explain it a little better because i mean i like the i if if there's something there we're not doing i like the idea of it and the thing that we are doing in in the deep the relative conformer energy analysis we currently have is that we look at the energy differences between conformers in m m and we look at the error of those relative to the qm energy differences for the same conformers which is a measure an attempted measuring how accurate the conformer energetics are but there may be other ways of doing that that we haven't yet hit on you were honest i haven't i've never had the time to look into these new these new energy measures that um you that were discussed and you just presented this you know stuff by um um so my my feeling was that i would i would i would love to have is kind of an energy rank or kind of a ranking is the ranking of the conformations in qm the same as it is in the in m m so are the is the lowest energy minimum the lowest energy in qm is the second lowest the second lowest in qm but probably this is something that's captured there you know i'm not sure that we've looked at ranking specifically so that might be something to check i don't know um what are the guys involved in that benchmarking thing yeah i think that'd be a good idea like a a row um or a towel like a ranking television that's actually that's a very good idea yeah it's it's a bit captured in the dde because if the ranking is wrong then the dds will be off quite well but yeah not completely but what i uh wanted to add is what i understood also the first from thomas is that it sounds a bit like ensembles sounds like md so it's it would be a more applied benchmarking that you for example run molecular dynamic simulations and then see uh whether you end up with different ensembles with different clusters of conformance which also would tell you whether the barriers are about the same size yeah almost the same height um which we don't capture right now in the benchmark so but of course it's there then a larger uncertainty where again it's a bit more like the protein ligand benchmarking because there's much more more work to be done kind of getting getting five five or ten confirmations for one molecule that's that's relatively easy but doing setting it up with md simulation do thing is the md simulation analyzing the md simulation i'm not quite sure how feasible that is for for thousands of compounds yeah so if there are no more questions i guess that we can move farther um yeah so uh yeah as you see we will have a five minute break uh before the interactive session uh what i was i want just to show you is to uh yeah this page to download the the material i hope that you already uh downloaded and did the the installation otherwise you can use this quick break to to do so uh yeah and essentially we will be looking to this particular chemistry that has been identified with this analysis it's um essentially out of sync torsion profile between the qm and of nff 1.3 torsion profile of this three-floor material pyridine and we will look into other desert containing this this chemistry to see whether we are able to to identify the the issue we can see we can be back at uh at 10 minute yeah 11 yeah let's uh restart this at 11 minutes past the hour okay yeah uh is there anyone is there anyone that needs uh output uh installation or running the the notebook essentially there are two ways to to running it uh you can run either on your local laptop or computer or by uh yeah install all the environment on a remote machine and then use a remote uptr notebook that will be the the the approach that i will be using so i can show you how to how to set up the the uptr notebook remotely because yeah i do need to do it myself and i'll also add for anyone running locally who's having trouble you can go ahead and um if you mention it in the chat i can pull you into a breakout room and help you get everything set up okay so uh first things i think i'm doing it's to create an ssh tunnel to the remote machine where the notebook will be running so to do so i do so essentially i'm here specifying the local port and the remote port local host is in here in the middle and then my credential and then you type just your username at the remote machine you are going to connect and might take a while okay now i cd into the directory that contains the the material and i conductivate the the environment uh my environment name is called a little bit differently but you should type konda activate yeah 2021 benchmark workshop and then once activated just type uptr then i'd like to specify no browser to avoid opening a browser window on the remote machine and also we need to specify the the remote port that should be the same that we specified here so everyone got the the same so when when you start running the uptr notebook you should get something like this essentially i'm going to copy paste these entirely on my browser i'm ready to offer the with uh with the uptr notebook uh yeah any question help needed to reach this point well if not i can go ahead so yeah okay welcome uh to this benchmark workshop live session uh as i mentioned uh we will be looking uh to uh this torsion shortcoming regarding this three floor medilla peridine group and yeah so the learning objective of this uh live session are to obviously identify the torsion violation that are contained in our data frame and then to identify which force field parameter are associated with these specific torsion violations uh and finally we want to visualize uh the the violation frequency so how many times this violation is uh committed as a function of the specific torsion parameter uh yeah well the notebook it's already compiled i will go briefly through it explaining it um you can stop me anytime for asking me question this is my first time ever uh doing a live coding session so uh please uh stop me at any time uh to have more uh yeah explanation if something it's uh not completely uh uh clear and yeah so we start loading this bunch of yeah could you zoom in a little bit so the text is larger oh yeah thank you great that's perfect okay yeah that's better thanks jeffrey uh yeah so we can start so uh essentially what you need to do it's to shift enter in each cell to execute it if you just type enter you see that you just get a new line so to execute each cell you need to to hold down the shift key and to present yeah by doing so in this first cell we load a bunch of python packages which are required to to run the the analysis so uh first part identify torsion violation uh so this first part of the code uh section uh we load uh our sdf files of the qm and the dmm methods then it will go through all the torsions of these qm and mm molecules to compare them to compare the qm and dmm values and will output the the conformers that exhibit torsion deviation which is larger uh to add event threshold that we will we will provide uh now since uh a given torsion can be defined well it's defined by for auto indices but actually the bond which is uh rotating is always uh defined by the two middle indices among this for torsional indices uh then a single rotatable bond might be identified as an outlier multiple times uh since it is central to multiple torsional to multiple torsions uh so the last part of this uh first code section uh essentially will will be used to remove this uh duplicates this uh this outliers so uh here we uh we define essentially uh a list of uh sdf files uh using this um yeah sorry no here we we we define the the location where the our sdf files are are stored and then uh we create uh a list to then go uh over uh with glob over this uh this folder containing all our sdf file uh we will use panda tools load sdf to load each single sdf as a separate data frame and then we will uh recursively uh append this uh this data frame to our qm data list so essentially this will be a list of uh data frames and then we uh concatenate this uh this data uh and we assign the the following uh columns to to this uh data frame and then we well this is just to modify the column name to call uh essentially conformer id this column conformer id rather than conformer index this column mall id rather than molecule index and and so on and then we do the same for the dmm data it is exactly the the same is done for the dmm data we concatenate also the dmm data we slightly modify the column of uh of this data frame and then we merge these two data frame based on the uh on the id column so at the end we we we get uh a data frame containing both our qm and our mm mm data uh yeah well here i'm just uh uh renaming this uh this column um this way so energy qm kigal per per mall and then i'm uh also assigning a float float number for for this for this particular value uh and then uh we we are sorting the the values of these resulting data frame that merged the qm and dmm data we sort the the values uh based on these two columns mall id and the energy of the qm method and then we group by uh by mall id mall id first so we can execute it you will see that nothing is happening but essentially we got this data frame which is our initial data frame we can have a look to it yeah it's not very informative because it contains uh yeah uh the the mall is uh yeah well we can forget about this it's not very informative uh but now we know that we have loaded all our data in this uh in this data frame and we can proceed with uh with the analysis so essentially this cell will uh it's the analysis the actual analysis part that will go uh through all the all the torsion to identify to essentially yeah looked to the to the qm structure and the mm structure and output all the conformer that uh exhibit torsion deviation of more than 30 degrees this is the torsion threshold that we are defining obviously you can play around and yeah lower and or increasing this this threshold but for now we will use 30 degrees as uh yeah torsion threshold and yeah essentially afterwards you can also consider uh uh uncommenting this this line that will filter out the the torsion that contains um this this atom uh and in in this way the code will not look for violation uh if the torsion ends with uh yeah hydrogen fluorine chlorine or the the atoms that are defined here you can play afterwards but uh for now we will uh go set for what with the with the analysis so what this uh cell and this code does so it goes uh over each rows of our data frame that stores the the qm and the mm data uh then we we define our uh mold one essentially will be uh our qm molecule mold two the uh open ff uh molecule we get the the properties of uh the qm molecule and then we get the the conformers of the of the two molecules and then we can uh search for all the dihedrals of ijkl uh indices uh by iterating over all the as i said before the middle indices so over all the jk uh bonds so uh we get we first iterate uh here over here getting all the bonds of the qm molecule and we get the atom j by getting uh the the begin of of the uh the yeah sorry the first atom involved in in this bond and the k atom so the one involved in this indices by getting the the end the ending of this of this jk bond and then uh well so we have essentially defined the bond now we we want to extend this bond uh to neighbors also to include the the the other bonds that will form the the the torsion so uh we go we iterate over the bonded neighbors of j uh that first that are not k so essentially we iterate over over here and then we do the other way around so we iterate over all the the the neighbors of k that are not j so essentially we we iterate over uh we iterate here with this part of the of the code and then well once we have uh essentially defined at the our our torsion we want to compute the the difference between the the qm torsion value and the dmm torsion value and we do uh this here essentially so we define t1 if you remember more one it's the qm molecule so t1 will be the qm torsion we get the the dihedral of the qm conformer and then we get all the the atom indices i j k and l we do the same for the dmm molecule and then uh essentially we want to test whether this the difference in in angle between the qm and dmm uh is uh exceed yeah it's more than than the given three shoulder uh and if this happened we we record the evaluation so we compute the the difference between the qm torsion and the dmm torsion uh yeah here essentially we are just rounding the uh by two decimal numbers by two decimal values and uh yeah this is a bit called artifact the actual difference might be more than 180 degrees and well if this is the the case we want to renormalize this this difference so we take the absolute value uh removing from this difference uh 36 uh 360 degrees um and yeah if uh the the difference it's uh exceed the our specified the three shoulder value then we we append uh to this um to this list we append the uh the name or the name of each i j k and l atom so the symbol of the atom and also the the indices we are here uh adding one value to the indices because essentially in the uh in python uh yeah number are zero basis uh whereas in yeah many of the program that you might use to visualize a molecule essentially when we check and we look to the labels we see one basis uh indices so i uh i think that it's um well it's uh added value to get uh yeah at the end of the different that contains uh the indices in one base at the number so that we can directly visually compare uh with uh a program like uh yeah uh apogadro or i don't know any any program pymol uh but then you will see that we will also use uh zero basis uh we will come back to also to zero basis uh indices when we will compare uh these indices with the the open ff uh parameter because the open ff parameter uh are phytonic and also uses uses um zero basis uh indices uh yeah so at the end we can uh convert this uh this list into a data frame we do this here and we specified uh the the the columns of this data frame uh id the name of the four atoms the indices of the four atoms and yeah the qm torsion value the mm or the open ff torsion value and then the difference um yeah well essentially i think we've got a pretty expert audience on and so uh it's probably okay if you want to move a little bit faster yeah yeah i'm see that timing is running out uh so we can have a look to to these uh data frame how it looks like so essentially we we have all the conformer the taxi be that uh torsion um so yeah torsion difference uh between qm and mm of more than 30 degrees but as we already mentioned uh there are torsion duplicates here for instance we see that the central bond that is rotating is is the same but just uh yeah this last atom indices is is changing so we want to get rid of these um uh farther in the wood we want to get rid of the of these duplicates we can uh actually preview uh which are the uh yeah the the molecules that are exhibiting a torsion violation for instance i'm here previewing uh uh the index index of the data frame 137 which would be this one so this this molecule we can have a look and and indeed you see that it it exhibits these uh difference between uh see in this cf3 group between the qm and open ff 1.3 optimized structures so uh yeah essentially what i do here is just to uh yeah rename some column and drop the the values of the qm torsion the ff torsion and the the difference since we are not interested anymore in in this and yeah well i here i i've put some challenge but maybe for the sake of time it is better if we skip this this part also the solution uh is already provided and was needed to further uh advance in in the code uh so we can essentially have a tuple containing now the the zero-basic uh indices so i'm removing one to to get back a tuple with the the indices of the torsion in zero base at the number because then we can further compare uh this tuple uh with the the open ff parameter indices and this uh this cell so this code block will actually remove the the duplicates the torsion duplicates that we we see we see here it will remove all this this kind of duplicates and i can finally have a data frame that contains the unique violation so each rotating bond now it's taken only only once and you see we we get this unique data unique violation data frame containing our conformer the uh the atom types the atom indices in one base at the indices to be comparable with any visualization visualization program and zero base at the indices that will be further compared with the parameter well also we can skip this this challenge you can play afterwards with the the notebook and yeah so in this uh second part of the of the code we won't so once we got our data frame containing all the unique rotating bond violation we want to essentially identify which force field parameter are associated with the specific torsion violation so in order to do so we build a dictionary containing all the parameters that are exercised by the conformer that produce a violation as we saw in the briefs data frame and afterwards again we have to we need to search the middle indices of the violating torsion that will match the middle indices of the torsion parameter to associate uh yeah the the the scene violation with the with the torsion parameter so to do so we need to uh to load the open ff uh molecule and topology and also to import our open ff uh force field so uh yeah we want to build this this dictionary and so we first load one conformer for each molecule that contains a violation here and then we load the this molecules in in upon the data frame and then essentially we yeah we we load our our force field and we create this uh this dictionary uh containing the dff parameter for all the unique conformer that produce a violation this will take uh a while uh this is obviously a tweaked and reduced the data set but if you execute this uh analysis and this you better not book over a real data frame containing uh a thousand maybe dozens of thousands of uh different conformer you will see that uh essentially this loading uh all the ff parameter will uh into a dictionary will take will take a while so now uh we want to essentially um associate to each uh torsion violation uh that was identified by the first part of the of the code to a force field parameter which is uh present in this in this dictionary that we just created and we do we do so in uh in this in this cell uh and we can have a look to what we we get so we have uh essentially yeah associated the the uh the violation uh identified by by the analysis exhibiting uh certain violent violation indices uh with the parameter uh so the force field parameter IDs uh not with the same indices but uh having so matching the two middle indices since we are interested uh obviously in the rotating bond so it might happen that for instance for a specific rotating bond there are uh more um ff parameter uh associated because the the terminal bonds are are different so uh here we match um any ff parameter that as the the the two middle indices and once we get these uh this data frame we can essentially visualize the the the result so we can visualize um the frequency the violation frequency as a function of torsion parameter we import the packages that is needed for uh for plots and we can uh essentially visualize this first plot uh which just takes into account the the number so the number of uh violation that are uh counted for a specific uh um open ff uh torsion parameter and here we see that uh t16 which is uh indeed the parameter that is associated uh with these um yeah this particular uh chemistry uh produced over uh 250 violation and essentially this is because the data frame is uh is triggered that and contains uh many of these uh of molecules containing these moieties on on purpose to to show you clearly the the result but obviously uh now we can argue that um yeah this plot is counting uh the the number of times uh a parameter as produced at the violation but at the same time this parameter uh has been uh used many times so it's uh this parameter has been over represented in in our data frame because uh our data frame contains many molecules with this uh with the moiety that is associated with this parameter so we are sure that t17 looks bad uh but we are not still sure it's because just t17 it's a parameter that has been assigned to too many torsion in in our data frame so we can normalize this violation count by the number of times that uh each parameter has been assigned and to do so uh we we can use the coverage uh report uh that we normally um we normally run uh in uh during our uh benchmarking uh and we can uh yeah we can load the the um um yeah the coverage report for the for the torsions and essentially we can finally normalize uh the absolute violation count dividing by the number of times the parameter has been uh exercised and multiplied to by 100 to have a kind of uh yeah relative normalized violation and if we now plot the result taking into account these sorry this normalized uh violation count uh we see so essentially uh here we see that only t17 is uh producing uh it's uh well concerning but uh once we normalize the uh the the violation count also by the by the number of times uh that the each parameter has been exercised in in our data frame well we can still see that the t17 uh is uh concerning parameter but there are also other parameters which are uh not very used very much in in our data frame uh but are producing a certain number of violations so at the end the the relative violation count will be um uh in some way will be uh also also i uh when uh so when we do this this plot uh more than having a normal bar plot uh we also want to uh to color the the bars uh with this palette so take into account uh again the the violation so the absolute violation count as we did uh essentially here so uh the colors will be assigned uh uh by the number of times that the the violation has been produced so that at the end uh let's say that we can we could focus on parameters that exhibit either um well both uh i relative violation but also produces so our color red red because they produce a um an i violation violation count number and yeah so we are really looking forward the working with this analysis and uh well one thing that i'm i'm doing uh right now that is not in the scope of obviously of this uber notebook is uh essentially to run a torsion scan uh we can run a torsion scan for uh over all the atom that are contained in parameters that we uh we think that might be concerning or we can simply uh run a torsion scan for uh all the torsion that are identified in the in by the the the analysis so uh yeah we really look forward to to it and yeah maybe in in the next near future uh we will show we do the the results of torsion scan uh benchmark done uh using these these analysis uh yeah i know i i was rushing a bit so i don't know if you have uh any any question regarding this uh notebook i will be happy to do yes lorenz i actually have one uh thanks a lot it was really nice and uh especially all these these the visualization and the the kind of pulling out the the torsion is something i never got around to do so now it's basically presented to me on a plate um so um it's the it's the interpretation just to be sure i understood what you do here that for t 17 you have uh 260 violations and 35 percent or 30 percent of the of the of the torsion 35 percent of the when this this this torsion shows up somewhere it's it's it's it's it's more than 30 degrees off from the qm yeah so uh regarding the absolute violation count so 260 this is exactly the number of times that uh this parameter has been uh seen in in in a torsion violation regarding the this actually percentage so the relative violation uh it's a bit controversial to to give uh uh let's say a straightforward the interpretation saying for instance that 30 percent of the time that the parameter so it's not actually exactly that 30 percent of that or 35 percent of the time this parameter um as produced a violation with respect to the times that has been exercised uh but it's uh kind of weighted so it's a kind of weighted matrix uh to normalize uh somehow by the number of times that uh the the the the parameter has been exercised so it has uh i'm afraid that i and i cannot provide uh what to say a physical straightforward uh um yeah but yeah Lorenzo the value would be one if it occurred no more often in molecules with violations than it occurs in the rest of the data set right sorry can you say again the value will be one i think the value on the vertical axis wouldn't that be one if it occurs at this in the same frequency in molecules with violations as in the rest of the data set yeah in principle uh yeah should be should be one i think it's like an over representation factor so it's 30 times more likely to occur in molecules with violations than in the data set as a whole so yeah i might add that there's two kind of confounding factor here so number one is uh this data set was selected to exhibit problems with t17 uh so in a in an unbiased data set this frequency wouldn't look so bad but this is kind of made for for ease of the workshop number two is that the number is a little more complicated than that because we're counting by a number of appearances of the torsion parameter so like the same torsion parameter could be assigned multiple times to the same central bond and so i'm i think it's a little more complicated than than kind of assuming that one is a baseline um i i would advise against trying to directly interpret this number but rather to focus on using it as a relative measure yeah i agree it certainly points you to the two areas where you you have to go for a deeper look yeah exactly and also yeah it's a kind of way to uh by the taking into account also by the number of times that the parameter has been exercised because uh otherwise we may end up with uh yeah some uh really high um number of violations but obviously that are due to the fact that this parameter has been used uh many many times so many times used many chances also to produce a violation uh whereas uh with this kind of normalization we want to um essentially yeah have a bit more clear representation of the of the picture anyway i i yeah no i i i would just want to say that uh this weight can be defined so i mean if anyone uh anybody finds uh any better way to define this weight it can be definitely uh weighted by well taking into account some different uh mathematics or um yeah so the other question i have is this t 157 this is kind of the i checked um above that's the generic atom one bonded to atom two bonded to atom three bonded to atom four so basic this is where the the force field doesn't produce anything reasonable mm-hmm so this is uh something where has anybody looked into what what this what weird bonds these are uh i guess we can have a look uh so i um yeah i'm i'm storing the this mix here and i guess this is a generic fallback parameter that if nothing else is found then it's still can uh i think it's further down oh yeah yeah here for instance yeah that's true yeah we should see what that gets used for but well it it should be i guess it should be seen uh which for what it's used among uh this um this data set so yeah i might not be able to to answer your question right now i think this was great um i think it would be great to have a something like a command line tool that runs over a full data set um that does just this because we probably have a larger a larger set of molecules right then we want to yeah the queuing system and maybe create st files with the qm and the mm examples so that we can overlay them and take a look definitely but the indeed the the idea would be to to use these analysis as a starting point to then uh uh launch uh uh at at ocean uh at ocean scan uh benchmarking with all the violation identified in the same manner as we did for the um yeah organ ff optimized execute so definitely this is something that we uh we would like to to have at some point yeah for for us because we have the qm data already but we probably have to pick out specific molecules and see if if the ip situation allows us to share it with you with you probably some tool that lets us run this on our whole data set and then visualize the outliers yes yeah i think so right now the benchmarking team's engineering effort is focused on um basically refactoring everything you may have noticed that some steps are like painfully slow and it's because for for scientific consistency in the data set we have to stick with the architecture we started out with uh and that ended up maybe not being suitable for some of the later analyses that we added so for the next few months we're going to focus on re-engineering things and then we'll be contacting you about a season two in 2022 and I think the violation analysis would be a great tool to include in season two so we're we're compiling requests for that now and I think the violation analysis command line tool will be a great one great any other discussion or rather that's the the end of your presentation right yeah yeah the end of the live session and yeah also the of the presentation so if there are that was really useful thanks Renzo yeah thanks well I I just want to thank you all for joining this workshop session yeah there will be a follow-up workshop in about two weeks if I remember well on bespoke feeding uh one week on one week sorry yeah same time slot next week yeah you are all invited uh and yeah so thanks again for for for joining uh and yeah well thanks also for to the open ff team for the support provided for yeah setting up this this workshop uh it has been awesome so yeah thanks again thanks for preparing things that's really helpful it was very helpful okay thanks all bye bye bye okay thank you bye