Questions regarding PBMETAD Parallel Bias Metadynamics

[KPC]

Dear all,

I am running a simulation using PBMETAD with Gromacs 4.6.7. Following is my PBMETAD options:

_________

wallcom: UPPER_WALLS ARG=cv4 AT=1.7 KAPPA=150.0 

PBMETAD ... 

ARG=cv1,cv2,cv3,cv5,cv6,cv7,cv8 SIGMA=0.1,0.1,0.1,0.1,1,1,1 HEIGHT=1.2 PACE=2500 LABEL=pb 

FILE=HILLS_1,HILLS_2,HILLS_3,HILLS_5,HILLS_6,HILLS_7,HILLS_8

BIASFACTOR=8 MULTIPLE_WALKERS 

GRID_MIN=0,0,0,0,0,0,0 GRID_MAX=3.142,8,1.1,3.142,40,20,20

INTERVAL_MIN=0,0,0,0,0.1,0.1,0.1 INTERVAL_MAX=3.1,8,1,3.1,38,18,19

... PBMETAD

PRINT ARG=cv1,cv2,cv3,cv4,cv5,cv6,cv7,cv8,wallcom.bias STRIDE=1000 FILE=COLVAR

_________

As I am using multiple walkers, I noticed that Plumed prints the value of all the replicas at any timestep into the same HILLS file, so there will be multiple lines for time 0 for example.

My question is, does sum_hills work correctly for these HILLS files? Or do I have to split my HILLS files into multiple files before running sum_hills?

Secondly, as I have not asked plumed to print the pbmetad bias into my COLVAR files, I will have to write a script to calculate the PBMETAD bias from my HILLS files. If my understanding is correct, the PBMETAD bias is just the sum of the exponential of the bias from all the HILLS files (with the Boltzmann factor, of course), right? And to reweight my trajectory I can simply multiple the conformations using the exponential of the PBMETAD bias according to the paper? 

Thank you for viewing and I appreciate your help!

[Massimiliano Bonomi]

Hi!

> My question is, does sum_hills work correctly for these HILLS files? Or do I have to split my HILLS files into multiple files before running sum_hills?

You don’t need to split the HILLS file into multiple files, you can use it as it is with sum_hills

>
> Secondly, as I have not asked plumed to print the pbmetad bias into my COLVAR files, I will have to write a script to calculate the PBMETAD bias from my HILLS files. If my understanding is correct, the PBMETAD bias is just the sum of the exponential of the bias from all the HILLS files (with the Boltzmann factor, of course), right? And to reweight my trajectory I can simply multiple the conformations using the exponential of the PBMETAD bias according to the paper?

One way to reweight PBMetaD runs is to use the final bias potential, as it were static along
the entire course of your simulation. This is what has been done in the first PBMetaD paper.
In order to use this approach, you don’t need to print out the instantaneous value of the bias potential.
So, no worries ;-)

Max

>
> Thank you for viewing and I appreciate your help!
>
>

> --
> You received this message because you are subscribed to the Google Groups "PLUMED users" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to plumed-users...@googlegroups.com.
> To post to this group, send email to plumed...@googlegroups.com.
> Visit this group at https://groups.google.com/group/plumed-users.
> To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/87f8bdc1-acac-465b-9928-fcc274ba9448%40googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[specialkender]

Hello everyone, 

Sorry to reopen an old post, but i have a similar question to the opener. 

I am running a simulation with one Umbrella sampling variable (a COM-distance) with 5 PBMetaD CVs (Well-tempered) and I want to unbias my ensemble from the PBMetaD CVs to then obtain a PMF over the Z CV using good old WHAM. Now i have a moment of confusion in the sense that i know how to do it with a classical Ndimensional gaussian potential, but I'm not exactly sure what form would the potential have in the PBMetaD.

In PBMetaD i have 1 potential for each of the CVs, so 5 VGs for 5 CVs all printed in each respective HILLS files. All the CVs are well tempered, so in the HILL files i know there is the gaussian height already scaled by the bias factor.
Now if I want to unbias the ensemble what should i do with the values in the hills files?   
In equation 8 from the Pfaendtner and Bonomi paper (https://www.plumed.org/doc-v2.6/user-doc/html/_p_b_m_e_t_a_d.html) it seems that they express the Parallel Bias potential as :

VPB(S1,S2,t)= -1/B log[ exp(-BVG(S1,t)) +exp(-BVG(S2,t)) ] 

Which would suggest just to get the exponential of the height column in the HILLS1, HILLS2...etc files, sum them and take the logarithm. Except that in the HILLS files i think there is the VG(S) multiplied by the bias factor. 

So the algorithm that seems to me right would be the following:

1 pick a time tk in which all my monodimensional potential has converged.

2 obtain the --negbias of every CV from the chosen time point tk so that the bias factor is taken into account.

3 do the exponential for each of the CVs, sum them up and take the negative logarithm divided by Beta. This is my VPB(S1,S2,t)
4 Reweight my ensemble from tk on with this VPB(S1,S2,t)

5 (case specific) Do the additional WHAM reweighting

What did I get wrong :P ?

 

[Massimiliano Bonomi]

Hi,

if you use plumed driver to re-read your trajectory and your HILLS file, you can directly print out the PBMETAD bias,
which is as you said equal to

> VPB(S1,S2,t)= -1/B log[ exp(-BVG(S1,t)) +exp(-BVG(S2,t)) ]

by using the PRINT instruction and mypbmetad.bias, where mypbmetad is the name you gave to the
PBMETAD action in the PLUMED input file. No need to calculate the individual biases and do the exponential summation by yourself.

This is quite similar to what we do with METAD in "Exercise 3: reweighting” in this tutorial:

https://www.plumed.org/doc-v2.6/user-doc/html/master-_i_s_d_d-2.html

Once you get VPB(S1,S2) you can calculate the PBMETAD “weights” as exp(+VPB(S1,S2)/kbT) and proceed with the WHAM

Max

> To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/040d003f-d650-4b58-9ef3-24f886279f20n%40googlegroups.com.

[specialkender]

Ok thank you for your clarification Max, that definitely saved me some time. I have one more doubt, previously I used c(t) to assess the convergence of the metadynamic and also to obtain the metad.rbias potential to reweight my histograms, so i would have something like: 

metad: METAD ARG=CV1,CV2 PACE=500 HEIGHT=2.0 SIGMA=0.03,0.03 BIASFACTOR=10.0 GRID_MIN=-1,-1 GRID_MAX=1,1 REWEIGHTING_NGRID=180,180 TEMP=298 FILE=HILLS  RESTART=YES 

w: REWEIGHT_METAD TEMP=298 ARG=metad.rbias

h: HISTOGRAM ARG=ALFA,BETA GRID_MIN=-1,-1 GRID_MAX=1,1 GRID_BIN=180,180 BANDWIDTH=0.05,0.05  LOGWEIGHTS=w

Now i know that in the new version of plumed REWEIGHTING_NGRID has been changed to RCT_USTRIDE, but running PBmetaD i cannot find any of the options to calculate the c(t). Is there any way to obtain this quantity? I know i could just use something like:

pbmetad: PBMETAD ARG=..... etc etc

w: REWEIGHT_METAD TEMP=298 ARG=pbmetad.bias

h: HISTOGRAM yadda yadda

and i might be overly cautious here, but I noticed that the final PMF that i obtain running WHAM on the  reweighted histograms is very sensitive to small changes in the parameters used to obtain the histograms (like for example change from the default normalization of the histograms to NORMALIZATION=TRUE), so I'm worried that using a bias rather than a rescaled bias ( bias-c(t) ) might affect my results.  

[Massimiliano Bonomi]

Hi,

as of today, the reweighing algorithm developed by Pratyush Tiwary et al. (the c(t) you are referring to) has not be extended to PBMETAD.
We discussed this possibility some time ago with Pratyush, but we did not push this forward.
Personally, with PBMETAD I use the umbrella-sampling static-bias exponential correction, using the PBMETAD bias obtained
at the end of the simulation, as described in my previous email. I simply obtain the weights from PLUMED and then build histogram/free energy
with external python scripts (not the HISTOGRAM functionalities).

Regarding the WHAM part, are you using the scripts provided on the PLUMED website?
If I remember correctly, these scripts are not using any HISTOGRAM functionalities either.

Max

> To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/f2cb280a-0046-4822-b9c6-630f5d80e39bn%40googlegroups.com.

[specialkender]

I see, i was fearing it. Well i guess I'll rerun the old calculations with the normal bias and see if there are any big differences in the resulting PMFs, hopefully not.

Regarding the WHAM part, are you using the scripts provided on the PLUMED website?
If I remember correctly, these scripts are not using any HISTOGRAM functionalities either.

->Nope, I'm using an in house code of WHAM that uses the histograms as input rather than the time series, as I was not really sure of how to reconstruct a proper time serie off the reweighted ensemble.   

Cheers,

Carlo

[Massimiliano Bonomi]

Well, I believe you can use the WHAM scripts in which the bias is assumed to be an harmonic potential centered at different positions
by adding also the V_PB contribution (calculated with the final HILLS file as the bias were static).
I think it is formally correct, apart from neglecting the time-dependence of the PBMETAD bias.
Maybe Giovanni Bussi has tried something similar with other flavors of METAD + umbrella sampling.

Max

> To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/1d6da4b2-ed90-4ef9-bed1-b255d3abc665n%40googlegroups.com.

[Giovanni Bussi]

Hi,

As far as I understand, Carlo wants to combine multiple simulations, each done with a restraint on a CV located in a different position, plus PBMETAD acting on another set of CVs. I am afraid this cannot be analyzed with a "standard" WHAM script on the first CV alone, since the replica will also have difference in the final PBMETAD bias for the other CVs. I would use the WHAM included in one of the PLUMED tutorials, which is binless and thus has maximum flexibility.

Giovanni

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/EEB3DD8A-5968-472D-B2B6-22539B903FEF%40gmail.com.

[Massimiliano Bonomi]

Well, I meant you can use our (PLUMED) WHAM scripts on the concatenated trajectory from all pbmetad/restraint simulations ;-)

Sent from my iPhone

On 6 Oct 2020, at 17:49, Giovanni Bussi <giovann...@gmail.com> wrote:



To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/CAPLm8ZK0YC6Y7%3DPvwhQZcgkMeBVOGgoE9CXGE9%2BpiEc2r6Te-w%40mail.gmail.com.

[Giovanni Bussi]

Ah sorry, yes that can be certainly done. Basically you just have to reprocess the whole traj using your separate plumed files, this will implicitly include also the (different) harmonic restraints.

In general, one should be careful when combining with WHAM trajectories that were obtained separately. If you are sure they are sampling well, the result would be correct. But if one of them is stuck, then there will be difficult-to-detect problems. In general, the same identical setting with replica exchange would be much more robust.

Giovanni

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/7B2FD951-0F5C-4D56-B223-8C62582B8726%40gmail.com.

[Massimiliano Bonomi]

I think you can add a COMBINE line to each PLUMED input file to sum up the PBMETAD bias and RESTRAINT bias and print it out.
This should be the bias needed in our WHAM scripts.

Max

[specialkender]

Thank you both very much for the constructive feedback, I'm still a little bit left wondering on some of the points so I'll try to clarify (even for future people to come :D). 

To clarify about my system: My system is made up of 15 umbrella window of a peptide each at restrained at a different z distance from the membrane with a RESTRAINT. Each of these window has also 5 coordination CVs to enhance orientational sampling. For each umbrella window i was going to unbias the trajectory of the window using REWEIGHT_METAD to obtain the unbiased histogram on z and run our wham script using all the Hi(z) (where i is the window) obtained this way. So let's say that nUS_CVs=1 (the z RESTRAIN) and nPB_CVs=5 (coordination cvs)

Now having said this i have some doubts:

 I would use the WHAM included in one of the PLUMED tutorials..

Any of them? wham.py that i can find in https://www.plumed.org/doc-v2.5/user-doc/html/lugano-2.html , or wham.cpp that i can find in https://www.plumed.org/doc-v2.5/user-doc/html/belfast-4.html are them all fine?

I think you can add a COMBINE line to each PLUMED input file to sum up the PBMETAD bias and RESTRAINT bias and print it out.

So does this mean that for each window (15) i should print a file with 1(time)+6(nUS_CVs+nPB_CVs)  + 1 (combined VPB+VRESTRAINT) =8 columns and feed these 15 files to the wham script, then to obtain th? 

Well, I meant you can use our (PLUMED) WHAM scripts on the concatenated trajectory from all pbmetad/restraint simulations ;-)

If the previous statement is true where does the concatenated trajectory come into play, don't i need only the COLVARs + (VPB+VRESTRAINT) ? (I am thinking at "concatenating" as the action of action of pasting the trajectories frame one after the other, which might be not what you meant) 

In general, the same identical setting with replica exchange would be much more robust.

As in that for every window the 5 coordination CVs would be instead run with replica exchange and the z still with a RESTRAIN? Why would this be much more robust?

Cheers,

Carlo

[Massimiliano Bonomi]

Hi Carlo, here a summary of the procedure that I would apply to your situation:

1) concatenate all your 15 trajectories (xtc, dcd, or any format supported by plumed driver) into one global trajectory (example: concatenated.xtc).
The order is not important

2) modify each of your 15 PLUMED input file used for production in order to:
- read the HILLS files produced in the corresponding run -> add RESTART keyword at beginning of file
- do not add additional HILLS -> set PACE of PBMETAD to a number larger than numbers of frames in traj_concatenated.xtc
- add a COMBINE line to sum up the PBMETAD and RESTRAINT bias
- add a PRINT line to print to a file (let’s call it COLVAR_WHAM.i, with i=1,15): time, nUS_CV, combined VPB+VRESTRAINT bias

3) run “plumed driver" on concatenated.xtc using each of your 15 PLUMED input files created at 2)

4) create a file (biases.dat) with 16 columns: time, 15* combined VPB+VRESTRAINT bias extracted from COLVAR_WHAM.i (3rd column)

5) run wham.py (https://www.plumed.org/doc-v2.5/user-doc/html/lugano-2.html):

python wham.py biases.dat 15 2.49

where 2.49 is the temperature in energy unit (here simulation is at 300K)

6) the output of wham.py is a file (weight.dat) with one weight per frame that can be used to calculate any possible property of the system.
For example, now you can:
- extract the time serie of nUS_CV from one of COLVAR_WHAM.i - they are all the same as far as nUS_CV is concerned!
- paste nUS_CV and weight.dat into a 2 column file
- use the script do_fes.py (https://www.plumed.org/doc-v2.5/user-doc/html/lugano-2.html) to finally calculate the free energy as a function of nUS_CV


I hope this clarify the procedure a bit more!

Max

> --
> You received this message because you are subscribed to the Google Groups "PLUMED users" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to plumed-users...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/17c96b8a-1848-4f8d-82ef-d7790a76c89cn%40googlegroups.com.

[Giovanni Bussi]

Regarding this point:

In general, the same identical setting with replica exchange would be much more robust.

The reason is that for "multiple restraint umbrella sampling" to work correctly you need that different simulations overlap in the conformational space. The standard check is that the restrained variable overlaps in neighboring windows. However, this is not sufficient. Basically, conformations from window i should be similar to conformations from window i+1 in all coordinates. This is not easy to check.

Therefore, a more reliable alternative is to run the windows simultaneously and allow exchanges of coordinates with replica exchange. You will be then able to follow the "demuxed" trajectories, that are continuous in the conformational space but move across the ladder of windows. If these trajectories behave in a diffusive manner, then by construction the windows will be overlapping (as they will originate from the same set of trajectories). The demuxed trajectories can be analyzed looking at any suspicious slow conformational degree of freedom to check if it's correctly sampled.

So, the advantage of replica exchange is that it is easier to detect problems, plus it usually converges faster. The disadvantage is that it requires you to run the simulations simultaneously, with some overhead.

Regards

Giovanni

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/AA2F66B1-2F4C-452B-B55F-B8CF475CBE39%40gmail.com.

[specialkender]

I hope this clarify the procedure a bit more!

That definitely does it, thanks a lot! I'll try it out and hopefully get it to work.

Basically, conformations from window i should be similar to conformations from window i+1 in all coordinates. 

Do you know of a paper that formally explain why this is true by any chance? I didn't come across a paper talking about this requirement and this might turn out to be a major concern for past simulations.