Problem converging to bonded potential (IBI)

[Cecília Álvares]

Hey there,

I am currently trying to derive bonded potentials of a very simple CG system (containing only one bond type and one angle type) using IBI. However, I have been failing miserably at doing it: instead of reaching potentials that are better and better at reproducing the target distributions for the bond and for the angle, I end up having weider and weider distributions as I do the iterations. I am posting a plot of the bond and angle distributions to give a glimpse on the "weirdness". I have already tried:

(1) providing very refined (small bin size and a lot of bins) target distributions of excelent quality (meaning not noisy at all) for the bond and the angle. Similarly, I have also tried using less refined target distributions (larger bin sizes and less amount of bins).

(2) varied a lot the setup in the settings.xml concerning bin sizes for the distributions to be built at each iteration from the trajectory file. I have tried very small bin sizes as well as large bin sizes.

(3) increasing the size of my simulation box hoping that maybe it was all a problem of not having "enough statistics" to build good distributions at each iteration within the trajectory file I was collecting from my simulations.

None of these things has worked and I think I ran out of ideas of what could possibly be the cause of the problem. Does anyone have any insights?

I am also attaching my target distributions (this is the scenario in which I am feeding target distributions lot of points and smaller bin size) and the settings.xml file for what is worth it.

[Marvin Bernhardt]

Hi Cecília,

I once encountered similar problems with bonded and non-bonded interactions. See Fig. 9 of this paper. In short: The problem was that the potential update of the non-bonded has some influence on the bonded distribution, and vice versa. But the potential update is calculated as if they were independent.

The fix in my case was to update the two interactions alternately using `<do_potential>1 0</do_potential>` for bonded and `<do_potential>0 1</do_potential>` for non-bonded interactions. You could try something similar.

Otherwise, is your system liquid? Are there non-bonded interactions that you are optimizing at the same time?

Greetings,

Marvin

[Cecília Álvares]

Hey Marvin,

Thanks a lot for the reply!

I will have a look on the paper right now and do some thinking. In fact, I wanted to test the possibility of optimizing the bonded potentials first and, after its optimization is done, optimize the non-bonded. So basically there is no optimization of non-bonded whatsover being done in my simulation. To build the target distributions, I sampled an atomistic system in which the non-bonded forces were artificially removed. After having a trajectory file of this AA system, I built the corresponding target distributions to be used in VOTCA with csg_stat. For what is worth it, the target distributions of angle and bond don't seem at all weird relative to the "real ones", of when non-bonded forces exist. And then, after having the target distributions, I set up the CG MD simulations within the IBI to have only bonded potential also. So, besides there being no non-bonded potential optimization, there is also no non-bonded forces at all in my CG system. But I dont think this should be a problem, right? It makes sense to entrust the CG bonded potentials to reproduce the target distributions of the AA bonded potentials.

What I did try also, and that is in allignment with your idea, was to set up two IBI runs: (1) one run to optimize only the potential for the bonds and keep the angle potential active (in this case the latter comes from a simple BI) and (2) one run to optimize only the potential for the angles and keep the bond potential active (in this case the latter comes from a simple BI). In the case (1) it seems that I converge to a potential for bonds that is quite able to reproduce the corresponding distribution, while in the case (2) I converge more and more to potentials that give super weird distributions (like with three weird peaks, as I showed in the figure above)

Concerning the phase of the system: it is a solid system. More specifically, it is a coarsened grained version of ZIF8 in which the whole repeating unit was assumed to be one bead. I know that IBI has not at all been developed for solids and even further not for MOFs - the goal is actually to derive potentials in the CG level using many different strategies (IBI, FM, relative entropy) and evaluate the results. In any case, I dont think that the fact that my system is a xtalline solid could be the reason why my results are super weird (right?). It seems like such a simple system when in the CG level.

For what is worth it, I am also assessing different mappings. Following the same strategy of optimizing first bonded-potential for a less coarsened mapping (2 beads), I am able to reach less weird results. For example, you can find below the evolution of the corresponding distributions as I perform more iterations for this system (it has one bond type and two angle types). I think there is still a problem since we can see some tendency of the distributions becoming non-smooth as I do more iterations, but the results are definitely less weird.

[Marvin Bernhardt]

Hey Cecília,

Oh ok, then it is probably not the interaction with the non-bonded terms, that causes issues. But I believe something similar is going on, that indeed has something to do with your system being a solid/crystal:

IBI is a very good potential update scheme, when the degrees of freedom are well separated. For molecules in liquids, angles and bonds are usually well separated, i.e. changing the potential of one, does not affect the dist of the other much. But multiple occurrences of equivalent DoFs also need to be well separated for IBI to work well. In your case, consider a single angle potential between three beads in the crystal is changed, but all the others are kept constant. It will change the distribution of that angle, but also have  effect on different angles. In that case IBI is not providing a good potential update at each iteration.

What is happening in detail, I believe, is that the angle potential of all angles is updated by IBI, but this leads to an “overshoot”. The next iteration, IBI tries to compensate, but overshoots again in the other direction. You can easily test if this is what is happening, plotting even and uneven iterations separately, i.e. compare a plot at iterations 10, 12, 14 with 11, 13, 15.

This has happened to me before with ring molecules, where the situation is similar. A possible solution is to scale the update, by some factor between 0 and 1 (I'd try 0.25).

Also test this for the bond potential, I guess this is happening there too, otherwise it should converge within ~20 iterations.

Greetings,

Marvin

[Cecília Álvares]

Indeed, this could be the reason why I have this weird non-smoothness in the plots I sent in my 2nd message (the ones concerning a less coarsened mapping), because indeed in this case I was optimizing all the three bonded potentials at once. I will try not doing them at the same time and see if the smoothness-issue improves.

But then this would not explain the issues I had in the original post I made, which concerned another mapping (a highly coarsened one). If the problem was a matter of optimizing more than one bonded potential at once, I should have had good results when I tried to do IBI only for one angle type and kept the potential for bonds constant (at a BI guess) throughout the procedure. But unfortunately my angle distribution still converges to something ultra weird with 3 peaks.

PS: maybe my last message was too big and maybe it was confusing, but the figures I sent in my 1st message and in my 2nd message are for different mappings. In the first one (let's call it mapping A), I have only 1 bond type and 1 angle type. For this one I did try optimizing separately to see if it would fix the problem and yet I reached weird results. The second message had figures of a less coarsened mapping (let's call it mapping B) in which I somewhat successfully converge to potentials that yield more or less rightful distributions (apart from the smoothness issue). I only brought up the results of the second mapping to show that the same strategy "worked" for deriving bonded potentials via IBI for another mapping. Sorry if I made it more confusing!

[Cecília Álvares]

(In any case let me try your factor idea, some other stuff that came to mind + finish reading your paper so that maybe I have more useful info on the problem)

[Cecília Álvares]

I think at this point I may be ready to just say that indeed IBI cannot be used to converge to a potential that is able to reproduce the structure of xtalline materials (or at least the material I am studying).

I've tried
(1) diminishing the factor used to update the potential (as you mentioned) and it did not work.
(2) updating literally only one potential at a time in the IBI and keeping the others literally constant either in the BI potential or in analytical forms that are able to reproduce perfectly the probability distributions. This would discard the possibility of dependence on the degrees of freedom in that sense that the update of one potential is affecting the distributions related to other potentials.
(3) Although the result is not meant to be bin-size-dependent, I tried playing with the bin size of both, the references I am feeding to VOTCA, and of the distributions it is meant to built as the iterative process runs for the different potentials. I thought maybe I was not setting up "proper" bin sizes for the algorithm.
(4) I tried dividing the angles lying within each of the two peaks in the initial figure I showed into two different angle types and it also did not work.
(5) I read your paper and tried to be more careful with issues that you raised in section 2.9 related to the smoothness of the distributions in the onset region (although VOTCA is supposed to take care of this internally apparently via the extrapolation methodology). Although section 2.10 bring up issues related to IMC, I also tried some more ideas that came to mind from reading that section and it didnt work.
(6) I've tried keeping analytical forms for the bonded potentials (I happen to have analytical forms that perfectly reproduce the distributions) and optimize the non-bonded and it also doesnt work.

Naturaly, in all cases, together with weird distributions, my potentials are also going to hell as the iterative procedure goes on (which explains why the corresponding distributions are weird).

For sure the problem doesnt have to do with the "sharpness" of the probability distribution curves (due to the xtalline material being highly ordered) cause I tried to feed "artificial" target distributions that are wide and thus less step and I dont converge to anything reasonable either.

Maybe the shape of the distributions for xtalline materials is not friendly to be used within IBI to converge to a potential, idk...
Well..

[Cecília Álvares]

Let me just ask one more question if I may:

In the section 2.9 of your paper, you talk about how the algorithm is set to create an "alternative RDF" which cherishes an interpolation in the onset region, where the values of the original RDF tend to be very small and the region tend to be poorly sampled (which is a quite good idea btw :) ). In the paper it specifically says range of values that you guys have had good experience with applying this interpolation procedure. In the abstract of the paper it says that the methods are implemented in VOTCA. Do you mean only the specific numerical methods you are using to do the iterative process or do you include also other specific things such as the interpolation protocol you described in section 2.9?

I am asking because in my case, sometimes, the distribution coming from the CG simulation ends up having small values that sometimes oscillates a bit back and forward in the onset region but the g(r) has values a bit larger than the value you mentioned in the paper for which the itnerpolation is done (1E-4). This causes weird potentials to happen which could be the reason why everything is going to hell. I am attaching a figure to illustrate the point. Is there a way in which I can change myself the value of the threshold for which I want to apply the interpolation? Maybe in my case I would need to use values higher than 1E-4. It could totally save the day and also make sense: since I am simulating a xtalline material whose superatoms are allowed less movement compared to a liquid, the setup of my interpolation needs to be more strict for the IBI to work.

[Marvin Bernhardt]

Regarding optimizing non-bonded potentials in crystals, just a list of things I would check:

Are the distributions at the iterations oscillating around the target distribution? Or is it rather a slow approach that never gets there? Or is it chaotic?

Are you working at constant pressure? If so, I would try at constant volume.

Is your thermostat working and your time step small enough such that the temperature is always as expected in each iteration?

Well possible, that it just does not work for your system, however, I am really surprised, that separating out a single potential in the whole system did not work.

Regarding the implementation in Votca:

It is still in the branch csg/mulit-iie2 at GitHub, you can build it from there. It has all the methods from the paper.

Regarding the bonded potentials:

For this situation it helps to restrict the range such that the problematic regions are not included. Votca should then extrapolate bonded potentials linearly.

[Cecília Álvares]

(1) indeed I spotted that in some cases they oscilate back and forth around the target distribution (I am attaching a pic as an example). However, this is not something that putting a factor < 1 was able to solve.

(2) no, I am working in the NVT ensemble.

(3) my thermostat is working: the temperature is quite well equilibrated (no weird spikes). The timestep us also small (I am using 5fs atm).

(4) me too :"D

R: Regarding the implementation in Votca: I saw that link in the paper. So indeed the interpolation scheme at the onset region that is mentioned in the paper is not implemented in the basic VOTCA installation and we need to use those codes in the branch you mentioned, right?

R: Regarding the bonded potentials: Good idea. That is actually something I did not try. I test it.

Photo below: evolution of the angle distribution in a scenario in which I am optimizing only one potential (i.e., the angle potential) + using a factor of 0.25

[Cecília Álvares]

PS: sorry, the y axis says g(r) but it is the angle probability distribution

[Marvin Bernhardt]

Regarding your last picture: I observe that all the even iterations (10, 30, 100) have spikes at different positions compared to the odd iterations (15, 19). I really would like to see a plot with consecutive iterations, i.e. 30-36 to see if it goes forth and back. However, this should be solved by a scaling factor. Did you try a small scaling factor like 0.1 or smaller?

I can offer you to run IBI on my computer and have a closer look. If you don't want to share the files here, you can also send me a direct E-Mail.

Cheers

[Cecília Álvares]

Hey Marvin,

In fact, maybe I am not setting the scaling factor correctly. I had seen page 56 of VOTCA's manual and understood the factor to be an option <scale> [value you want] </scale> " that should be input inside the <post_update_options> (which is inside <inverse>). I mean, it does say in page 56 of the manual,  "post_update_options.scale : scale factor for the update (default 1.0)". But now that I took a look comparing the results with and without the factor, instead of simply looking at them two separately, I realize that the curves are exactly the same. So the factor that I am setting is not doing anything....

Sure, I have no problem sharing the files here. If you want to run, I suggest that you try optimizing the bonded potentials (and also only the angle, leaving the bond potential constant throughout the IBI) because it is simpler than doing it for the g(r). I will put it all here in a zip file. Thanks a lot for the offer.
PS: In fact, after I re-read your previous email, I realized I had misread the first sentence of your phrase: in my case, I was much more surprised that the optimization of the bonded didnt work. The g(r)s have very complicated shapes, so in fact for me it is more shocking not to be able to reproduce the angle distribution than the g(r) - at least not without the interpolation in the onset regions you mentioned in the paper.

Based on your previous suggestion: yesterday I tried narrowing the min and the max to do no accomodate the onset regions and splitting the angles into two types, but upon looking at the results quickly, the iterations are not getting better either, but then I need to analyse this more carefully still.

[Cecília Álvares]

Also, here comes potential curves and probability distribution curves that are of consecutive steps as you asked. You can see that indeed it is going back and forth (at least in these three steps that I prepared). And despite understanding how the "pre-factor" idea can help with the cause, I dont think that it will really save the day without implementing your codes to take care of the onset interpolation. This is because in my case I have a lot of onset regions with very tiny values which will result in huge values when boltzmann inverted (specially in the g(r), if I go for optimizing non-bonded later on, since my material is xrystalline).

You can see that I have these weird peaks that are artificially created without the interpolation. I think this is what causing everything to go to hell. But I havent managed to use the codes you mentioned in that VOTCA branch to give you a feedback.

Anyways, thanks a lot for helping me with all of this ! :)

[Cecília Álvares]

(more pictures lying of distributions in farther steps.

PS: In the figure that I showed you before the last message you sent, step 30 is a weird distribution whilst in the image below it isnt. This is because I ran the simulation to generate the results shown above with a smaller amount of atoms in order to get it done faster. That's why they dont match. In any case, the problem seems to be there regardless of how big the box is and I do gather enough microstates to do the statistics.

These curves here and the one from my message above with steps 5,6 and 7 are from a same simulation containing the same amount of superatoms as in the file I sent you (1500), and with optimization occuring only in the angle, so that maybe it is more comparable with your case.

[Marvin Bernhardt]

Hi Cecília,

unfortunately I did not have time to run it myself. But I had a quick look at the files and I think I see what is missing to activate post-update scaling. It should have "scale" in post_update:

```     

<post_update>scale</post_update>
<post_update_options>
  <scale>0.25</scale>
</post_update_options>
```

Maybe you can try that out, I keep my fingers crossed :)

[Cecília Álvares]

Hey Marvin,

It seems to have worked! I guess it was not very intelligent of me to have taken for granted that my previous setup for the scaling factor was correct: if I had compared the results with vs without my factor instead of analyzing it isolate to see if the problem had been solved, I would have noted long ago that there the factor setup was not on and saved some time and energy :"D

It seems though that I will indeed need to use your codes to do the interpolation/extrapolation mentioned within the section 2.9 of the paper you sent since my potential curves are still awkward in the angle range corresponding to the onset regions of the rdf. So that will definitely save the day also - thanks for sending me that paper.

Other than that, hopefully IBI will keep working later on, when I move to the non-bonded in other mappings cherishing more than one bead type, which will be more challanging since the effect one pair potential has in the spatial arrangement of other pairs is meant to be more intense in xtalline solids compared to liquids.

In any case, thanks a lot once again! It is no scientific prize, but your name will totally be in the acknowledgements of my thesis!! :D

[Christoph Junghans]

On Mon, May 8, 2023 at 12:44 PM Cecília Álvares <cecilia...@gmail.com> wrote:

Hey Marvin,

It seems to have worked! I guess it was not very intelligent of me to have taken for granted that my previous setup for the scaling factor was correct: if I had compared the results with vs without my factor instead of analyzing it isolate to see if the problem had been solved, I would have noted long ago that there the factor setup was not on and saved some time and energy :"D

It seems though that I will indeed need to use your codes to do the interpolation/extrapolation mentioned within the section 2.9 of the paper you sent since my potential curves are still awkward in the angle range corresponding to the onset regions of the rdf. So that will definitely save the day also - thanks for sending me that paper.

Other than that, hopefully IBI will keep working later on, when I move to the non-bonded in other mappings cherishing more than one bead type, which will be more challanging since the effect one pair potential has in the spatial arrangement of other pairs is meant to be more intense in xtalline solids compared to liquids.

In any case, thanks a lot once again! It is no scientific prize, but your name will totally be in the acknowledgements of my thesis!! :D

Celillia, sorry for being late to the party! I am glad you found the problem! And thanks Marvin for helping out here.

There is one example in the tutorials about scaling,see:

https://github.com/votca/votca/blob/master/csg-tutorials/urea-water/ibi/settings.xml#L52C1-L58

And the VOTCA pdf manual is obsolete, please use our online documentation at https://www.votca.org instead.

Christoph

--
Join us on Slack: https://join.slack.com/t/votca/signup
---
You received this message because you are subscribed to the Google Groups "votca" group.
To unsubscribe from this group and stop receiving emails from it, send an email to votca+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/votca/220aeac6-1ae0-4ad7-ae1e-6124aa435469n%40googlegroups.com.

--

Christoph Junghans
Web: http://www.compphys.de

[Sanjeet Singh]

Hello Cecilia,

I am trying to follow your procedure of optimizing the bonded potential first. I am also using LAMMPS for running all my atomistic simulations.

I am bit confused with removing the non-bonded forces artificially. Are you initially running your simulations with all the interactions (bonded + nonbonded) present in the system, and then using rerun command to artificially remove the non-bonded forces from your system and create a trajectory with only the bonded forces.

Greetings,

Sanjeet

[Cecília Álvares]

Hello Sanjeet,

So, in the end of the day, I, myself, ended up using IBI only for the non-bonded. I am using a simple BI to derive the bonded potentials.

I never explored this rerun command to remove non-bonded forces. I am also skepitcal on how this would work (although I never read much about it). When I was with that plan of optimizing bonded first in head, I ran a simulation in the AA level without the non-bonded force-fields (so only with bonded), saved a bunch of microstates of the equilibrated system, and derived the reference bond and angle distributions for the CG representation from these saved AA microstates. Note that this may not work for all the systems universely: sometimes the resulting reference CG structure you get from simulating the AA only with bonded forces is very "irrealistic" when compared with the actual, rightful structure your system should have in the corresponding CG level (i.e., the one corresponding to when you simulate the AA with both bonded and non-bonded). For example, angle distributions in the CG level that should have 2 peaks when simulating the AA with the complete set of force-fields may turn out to have only one peak (as if the two peaks merged) in the AA simulated only with bonded force fields. So if you see that the structural results you get from running the AA simulation only with bonded force-fields are very unrealistic compared to the rightful structure your system should have, I would not recommend doing the optimization of the CG bonded potentials using those structural results of AA sampled with bonded contributions only. It is better to use as reference the structural results of the AA system simulated with bonded and non-bonded.

[Sanjeet Singh]

Hello Cecilia,

Thank you for the info. It is very helpful indeed.