Regarding selection on min,max for getting the bonded potential using Boltzmann Inversion

[Sanjeet Singh]

Hello Christoph,

I have run a 1000 ns atomistic vacuum simulation to get bonded distribution with exclusions obtained using the following command:

csg_boltzmann --top topology.xml --cg mapping.xml --excl exclusions.txt

Then, I am getting my  bonded distribution using the following command:

csg_boltzmann --top topology.xml --trj Trajectory2.dump --cg PEO-10.xml --first-frame 500000 < boltzmann_cmds_dist

After getting my bonded distribution I want to use Boltzmann Inversion to get my bonded potential using the following command:

csg_boltzmann --top topology.xml --trj Trajectory2.dump --cg PEO-10.xml --first-frame 500000 < boltzmann_cmds

Now, suppose I want to do the Boltzmann Inversion for a Bond-AB for which I am using the boltzmann_cmds file with the following options:

tab set T 300
tab set scale bond    #SCALING IS NEEDED TO GET THE VOLUME NORMALIZED DISTRIBUTION FUNCTION i.e. dividing the distribution by 4pir^2
tab set auto 0
tab set min ______
tab set max _____
tab set n 100
tab bond-AB.potential-file-100.ib *bond-AB*

Then in the above lines, what should be the value of min and max? 

I am attaching the bonded distribution file for BOND-AB here, and also the plot of BOND-AB distribution.

If I go by the BOND-AB distribution file, then maybe I can use min = 0.348425, and max = 0.369061. As before this min and after the max value, I have 0 probability for some bond length.

Or, should I select min and max from the plot of the distribution? And select min = 0.351, and max = 0.367, i.e. the range where I have the normal distribution in the plot.

 Thank you. I would appreciate any advice you have for me.

 

[Christoph Junghans]

Short answer is for csg_botlzmann it doesn't really matter.
However, for csg_inverse the min and max matter, so the boltzmann
inversion from the target distribution (<name>.dist.tgt) is done
correctly, if you don't provide an initial guess for the potential
(<name>.pot.in)

Christoph

>
> Thank you. I would appreciate any advice you have for me.
>
>

> --
> 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/6187e55a-2925-4a5c-992e-8f3c44112b1bn%40googlegroups.com.



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

[Sanjeet Singh]

Hello Christoph,

Thank you for your prompt response. But, I am not able to follow you when you say  "for csg_botlzmann it doesn't really matter".

For example in my dist.tgt file here, I have initial data points as follows:

0.347465    0.00833478
0.347705    0.00833478
0.347945    0.0166696
0.348185    0

0.348425    0.00833478

 Therefore, when I am using this dist.tgt file to get bonded potential using Boltzmann Inversion I am setting the min at 0.348425 as before this I am having a point with zero probability.

Or if I set my min at the starting point 0.347465, then will this be OK? 

[Christoph Junghans]

On Tue, Oct 18, 2022 at 6:30 PM Sanjeet Singh <sanje...@gmail.com> wrote:
>
> Hello Christoph,
>

> Thank you for your prompt response. But, I am not able to follow you when you say "for csg_botlzmann it doesn't really matter".

What I meant there was that csg_boltzmann will invert things either
way, but what I would do is to pick min and max inside the peak of the
distribution, so 0.351 and 0.367, as you suggested, should work.

Christoph

> To view this discussion on the web visit https://groups.google.com/d/msgid/votca/60ada3ad-6d4f-473c-91ca-1a703f12891en%40googlegroups.com.

[Sanjeet Singh]

Dear Christoph,

I would like to ask a follow-up question to this. Now, suppose I set my min = 0.351 and max = 0.367, then I will be getting my potential in this range. What if I want to extend my potential outside this range, Say for example the range of my dist.tgt file. The range of my dist.tgt file is 0.347465 to 0.371461. It's just that in the boundary regions of dist.tgt file I don't have a proper sampling.   

[Christoph Junghans]

On Wed, Oct 19, 2022 at 9:24 AM Sanjeet Singh <sanje...@gmail.com> wrote:
>
> Dear Christoph,
>
> I would like to ask a follow-up question to this. Now, suppose I set my min = 0.351 and max = 0.367, then I will be getting my potential in this range. What if I want to extend my potential outside this range, Say for example the range of my dist.tgt file. The range of my dist.tgt file is 0.347465 to 0.371461. It's just that in the boundary regions of dist.tgt file I don't have a proper sampling.

csg_inverse would to the extension automatically, for manually modification see:
https://www.votca.org/csg/preparing.html?highlight=extrapolate#post-processing-of-the-potential

Christoph

> To view this discussion on the web visit https://groups.google.com/d/msgid/votca/86dd6512-6372-4061-a1d8-130f429ebc6fn%40googlegroups.com.

[Sanjeet Singh]

Hello Christoph,

Sorry for the incessant questioning.

I got the point that csg_inverse will do the extension automatically. But when I am using Boltzmann Inversion (not IBI) to get the bonded potential, then do I need to do the modifications (extending the range of the potential) manually by using the post-processing tools in VOTCA? Or can I use csg_inverse for simple Boltzmann Inversion also? 

Thank you.

Sanjeet

 

[Christoph Junghans]

On Wed, Oct 19, 2022 at 13:00 Sanjeet Singh <sanje...@gmail.com> wrote:

Hello Christoph,

Sorry for the incessant questioning.

I got the point that csg_inverse will do the extension automatically. But when I am using Boltzmann Inversion (not IBI) to get the bonded potential, then do I need to do the modifications (extending the range of the potential) manually by using the post-processing tools in VOTCA? Or can I use csg_inverse for simple Boltzmann Inversion also? 

yes, when you do BI then you need to do it manually! Of course you can do IBI for just one step and steal the potential from there ;-)

Christoph

To view this discussion on the web visit https://groups.google.com/d/msgid/votca/96264f49-d275-440a-9ee3-76dece6c69ecn%40googlegroups.com.

[Sanjeet Singh]

Dear Christoph,

Thank you for the clarifications.

I have another question. I was just verifying by hand the bonded potential which we get by using the Simple Boltzmann Inversion method. Say for example for a Dihedral distribution. The first few data points from dist.tgt file is shown below:

-3.14159     0.456519
-3.07876     0.897376
-3.01593     0.866437
-2.95309     0.825311
-2.89026     0.760535

Now, for the first data point if I try to get the potential by hand by using the formula U(q) = -kBT ln P(q), I am getting the value 1.95587. But the value which VOTCA gives me is 1.72411.

Can you please explain the reason for this discrepancy?

Here are the first few data points from the potential table:

-3.14159        1.72411                     40.8857
-3.09942        0                                  20.0311
-3.05725        0.0347262                 -1.32404
-3.01508        0.111667                   -1.92226
-2.97291        0.196846                    -2.01054

Thank you.

Sanjeet

[Christoph Junghans]

On Thu, Oct 20, 2022 at 6:25 AM Sanjeet Singh <sanje...@gmail.com> wrote:
>
> Dear Christoph,
>

> Thank you for the clarifications.
>
> I have another question. I was just verifying by hand the bonded potential which we get by using the Simple Boltzmann Inversion method. Say for example for a Dihedral distribution. The first few data points from dist.tgt file is shown below:
>
> -3.14159 0.456519
> -3.07876 0.897376
> -3.01593 0.866437
> -2.95309 0.825311
> -2.89026 0.760535
>
> Now, for the first data point if I try to get the potential by hand by using the formula U(q) = -kBT ln P(q), I am getting the value 1.95587. But the value which VOTCA gives me is 1.72411.
> Can you please explain the reason for this discrepancy?

It seems the grid spacing is different and hence some kind of
interpolation and /or spline fitting might have happened during the
procedure. (For dihedrals there is even some consideration about the
periodicity)
Usually VOTCA tells you all the steps it is doing along the way, but
without these details it is hard to say for me right now what exactly
triggered the difference.
If you post the commands you ran and output you got here we can have a
look at the details.

On a related note: dihedral interactions are a bit tricky due to
periodicity and grid spacing and when you search this mailing list you
will find some other cases where VOTCA wasn't smart enough to handle
all cases and users had to fail back to manual inversion.

Christoph

> To view this discussion on the web visit https://groups.google.com/d/msgid/votca/f954a415-b916-49e7-aace-d7ab2ab41aa6n%40googlegroups.com.

[Sanjeet Singh]

Hello Christoph,

For calculating the dihedral potential I have just use the following command:

csg_boltzmann --top topology.xml --trj Trajectory2.dump --cg PEO-10.xml --first-frame 500000 < boltzmann_cmds

Where the boltzmann_cmds file I am using is following:

tab set T 300
tab set scale no     #DIHEDRAL DISTRIBUTION DOESNOT NEED TO BE NORMALIZED
tab set auto 0
tab set min -3.07876
tab set max 3.07876
tab set n 150
tab dihedral-ABBB.potential.ib *dihedral-ABBB*

To be on the same page, when I am using the potential of mean force formula U(q) = -kBT ln P(q), I am using kB value as 0.0083144  kJ/mol K.

But, I also think that the deviation will be mainly because of the different grid spacing and some fitting which might have taken place.

Also, can you please point out What do you mean by manual inversion of the dihedral interactions?

I mean after performing the Simple Boltzmann inversion to get the dihedral potential, postprocessing of the potential should be done manually by running each command one by one.

Is that what you mean?

Thank you.

Sanjeet 

[Christoph Junghans]

On Fri, Oct 21, 2022 at 7:49 AM Sanjeet Singh <sanje...@gmail.com> wrote:
>
> Hello Christoph,
>
> For calculating the dihedral potential I have just use the following command:
>
> csg_boltzmann --top topology.xml --trj Trajectory2.dump --cg PEO-10.xml --first-frame 500000 < boltzmann_cmds
>
> Where the boltzmann_cmds file I am using is following:
>
> tab set T 300
> tab set scale no #DIHEDRAL DISTRIBUTION DOESNOT NEED TO BE NORMALIZED
> tab set auto 0
> tab set min -3.07876
> tab set max 3.07876
> tab set n 150
> tab dihedral-ABBB.potential.ib *dihedral-ABBB*
>
> To be on the same page, when I am using the potential of mean force formula U(q) = -kBT ln P(q), I am using kB value as 0.0083144 kJ/mol K.
>
> But, I also think that the deviation will be mainly because of the different grid spacing and some fitting which might have taken place.
>
> Also, can you please point out What do you mean by manual inversion of the dihedral interactions?
>
> I mean after performing the Simple Boltzmann inversion to get the dihedral potential, postprocessing of the potential should be done manually by running each command one by one.
> Is that what you mean?

Yes, sorry I meant the post processing of the potentials, that is
described here:
https://www.votca.org/csg/preparing.html?highlight=extrapolate#post-processing-of-the-potential

> To view this discussion on the web visit https://groups.google.com/d/msgid/votca/9049da08-df98-465c-92cf-dab49c10b434n%40googlegroups.com.