Reweighting from std. Metadynamics simulation

[Adip Jhaveri]

Hello Plumed users and developers,

I have performed a metadynamics simulation using two CVs called (d1,d2) and I am trying use the reweighting method to calculate average potential energies of the system and also the FES w.r.t an unbiased CV. (I following the protocol as given in Tiwary, Pratyush, and Michele Parrinello. The Journal of Phy. Chem. (2015)). I am facing a problem regarding this.

In my plumed file (used for the simulation), I did not calculate the reweights (i.e. c(t)). So I have to calculate this by using the driver tool on the trajectory. (plumed driver --mf_trr *file.trr --plumed plumed_post.dat). The problem with this is I have recorded the positions every 200 ps but the bias was deposited every 140 ps in the simulation. Therefore during post-processing if I recalculate the biases and c(t), it will not match with that of the simulation biases. Is there a way around this, perhaps using the COLVAR file (written every 20 ps) or it is best to calculate the reweighting factor by writing a script?

Regards,

Adip    

[Omar Valsson]

Hello Adip, 

You can use the READ action (https://www.plumed.org/doc-v2.6/user-doc/html/_r_e_a_d.html) to read values from the COLVAR file using the driver. You should then start the driver using the "--noatoms” flag and specify the time step used using “--timestep” (so that the stride used in the plumed.dat). Just write out all the values (e.g. CVs and bias values) to another colvar file and compare to the original one to make sure everything is consistent. 

Regards,

Omar  

--
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/e4418b90-3b54-4fad-bacf-0d3eb348b70d%40googlegroups.com.

[Adip Jhaveri]

Thank you Omar. It works this way.

On Wednesday, April 29, 2020 at 4:55:24 AM UTC-4, Omar Valsson wrote:

Hello Adip, 

You can use the READ action (https://www.plumed.org/doc-v2.6/user-doc/html/_r_e_a_d.html) to read values from the COLVAR file using the driver. You should then start the driver using the "--noatoms” flag and specify the time step used using “--timestep” (so that the stride used in the plumed.dat). Just write out all the values (e.g. CVs and bias values) to another colvar file and compare to the original one to make sure everything is consistent. 

Regards,

Omar  

On 29. April 2020 at 06:59:44, Adip Jhaveri (jhave...@gmail.com) wrote:

Hello Plumed users and developers,

I have performed a metadynamics simulation using two CVs called (d1,d2) and I am trying use the reweighting method to calculate average potential energies of the system and also the FES w.r.t an unbiased CV. (I following the protocol as given in Tiwary, Pratyush, and Michele Parrinello. The Journal of Phy. Chem. (2015)). I am facing a problem regarding this.

In my plumed file (used for the simulation), I did not calculate the reweights (i.e. c(t)). So I have to calculate this by using the driver tool on the trajectory. (plumed driver --mf_trr *file.trr --plumed plumed_post.dat). The problem with this is I have recorded the positions every 200 ps but the bias was deposited every 140 ps in the simulation. Therefore during post-processing if I recalculate the biases and c(t), it will not match with that of the simulation biases. Is there a way around this, perhaps using the COLVAR file (written every 20 ps) or it is best to calculate the reweighting factor by writing a script?

Regards,

Adip    

--
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...@googlegroups.com.

[Adip Jhaveri]

Hi Omar,

I have followed the procedure you described to calculate the reweights. However, I am facing a problem. After a certain number of steps, the reweighting factor c(t) is calculated as infinity. Therefore I am not able to get the reweights for the entire simulation. Can you help me out with this error?

Regards,

Adip

On Wednesday, April 29, 2020 at 4:55:24 AM UTC-4, Omar Valsson wrote:

Hello Adip, 

You can use the READ action (https://www.plumed.org/doc-v2.6/user-doc/html/_r_e_a_d.html) to read values from the COLVAR file using the driver. You should then start the driver using the "--noatoms” flag and specify the time step used using “--timestep” (so that the stride used in the plumed.dat). Just write out all the values (e.g. CVs and bias values) to another colvar file and compare to the original one to make sure everything is consistent. 

Regards,

Omar  

On 29. April 2020 at 06:59:44, Adip Jhaveri (jhave...@gmail.com) wrote:

Hello Plumed users and developers,

I have performed a metadynamics simulation using two CVs called (d1,d2) and I am trying use the reweighting method to calculate average potential energies of the system and also the FES w.r.t an unbiased CV. (I following the protocol as given in Tiwary, Pratyush, and Michele Parrinello. The Journal of Phy. Chem. (2015)). I am facing a problem regarding this.

In my plumed file (used for the simulation), I did not calculate the reweights (i.e. c(t)). So I have to calculate this by using the driver tool on the trajectory. (plumed driver --mf_trr *file.trr --plumed plumed_post.dat). The problem with this is I have recorded the positions every 200 ps but the bias was deposited every 140 ps in the simulation. Therefore during post-processing if I recalculate the biases and c(t), it will not match with that of the simulation biases. Is there a way around this, perhaps using the COLVAR file (written every 20 ps) or it is best to calculate the reweighting factor by writing a script?

Regards,

Adip    

--
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...@googlegroups.com.

[fede...@gmail.com]

Hello Adip,

I hope you solved yuor problem. But if it is not the case I have a sort of solution which worked for me, because I had the same problem about a month ago.

I corrected the routine which calculate the reweighting factor c(t). In the MetaD.cpp plumed file (in the version that I am using) there is this part:

**********

void MetaD::computeReweightingFactor()
{
  if(biasf_==1.0) { // in this case we have no bias, so reweight factor is 0.0
    getPntrToComponent("rct")->set(0.0);
    return;
  }

  double Z_0=0; //proportional to the integral of exp(-beta*F)
  double Z_V=0; //proportional to the integral of exp(-beta*(F+V))
  double minusBetaF=biasf_/(biasf_-1.)/kbt_;
  double minusBetaFplusV=1./(biasf_-1.)/kbt_;
  if (biasf_==-1.0) { //non well-tempered case
    minusBetaF=1;
    minusBetaFplusV=0;
  }
     const double big_number=minusBetaF*BiasGrid_->getMaxValue(); //to avoid exp overflow

  const unsigned rank=comm.Get_rank();
  const unsigned stride=comm.Get_size();
  for (Grid::index_t t=rank; t<BiasGrid_->getSize(); t+=stride) {
    const double val=BiasGrid_->getValue(t);
    Z_0+=std::exp(minusBetaF*val-big_number);
//    Z_V+=std::exp(minusBetaFplusV*val-big_number);
    Z_V+=std::exp(minusBetaFplusV*val);
  }
  if (stride>1) {
    comm.Sum(Z_0);
    comm.Sum(Z_V);
  }

//  reweight_factor_=kbt_*std::log(Z_0/Z_V);
  reweight_factor_=kbt_*(big_number+std::log(Z_0/Z_V));
  getPntrToComponent("rct")->set(reweight_factor_);
}
***********

I have corrected the red lines with blue lines. It is just a semplification but it works (at least for me).

I hope this is useful,

Best regards,

Federica

[Michele Invernizzi]

Hi Federica,

What was the value of BIASFACTOR when you had an infinite c(t)?

This implementation (my fault) is not robust in case of non-tempered metad, and maybe also in case of a high BIASFACTOR... 

In my defense I can say that the previous implementation did not support non-tempered metad at all :)

Anyway, it shouldn't be too complicated to fix this, I'll take a closer look and make a pull request.

Thanks for pointing it out.

Regards,

Michele

[Pratyush Tiwary]

Ciao Michele, Federica

For non-tempered metadynamics to get c(t) I would suggest doing (eq 18 of valsson et al https://www.annualreviews.org/doi/pdf/10.1146/annurev-physchem-040215-112229 )

where puting $F(s)=-V(s,t)$ as the running estimate one gets

I believe the denominator here can be ignored. This c(t) should increase linearly in t, unlike well-tempered where the increase is logarithmic in t.

Pratyush

--

Pratyush Tiwary

Assistant Professor, Department of Chemistry and Biochemistry & Institute for Physical Science and Technology 

University of Maryland, College Park, MD 20742

sites.google.com/site/pratyushtiwary/

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/2f8026b0-8358-4f5b-92d4-dfc634cb4601n%40googlegroups.com.

[Michele Invernizzi]

Ciao  Pratyush,

I totally agree with you! I was able with a small modification to avoid the rct overflow no matter if one is using well-tempered metad or not, and made a pull request  https://github.com/plumed/plumed2/pull/599 

Regards,

Michele 

[Federica Lodesani]

Hello Michele and Pratyush,

I used a bias factor of 100 (I also tried non tempered but the "infinity" problem was solved by changing the plumed version).

Thanks for your quick response and for the modification.

Regards,

Federica

You received this message because you are subscribed to a topic in the Google Groups "PLUMED users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/plumed-users/maD46ZOfss4/unsubscribe.
To unsubscribe from this group and all its topics, send an email to plumed-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/7938859c-7e4e-4eae-909b-4153dba184bcn%40googlegroups.com.

[Adip Jhaveri]

Hello everybody, thanks for the help!

Federica, thanks for the solution. This solves the infinity problem for me. Although I do have a little doubt regarding the code. For the non-tempered case (biasf_ = -1), the minusBetaF is set to 1 and Z_0 is then calculated. But the kbt_ factor is not multiplied before taking the exponential. Am I missing something?

Pratyush,

While calculating the c(t) you have suggested to use the following expression - 

I was wondering how c(t) would change if I use the expression given in https://pubs.acs.org/doi/full/10.1021/jp504920s - (for non-tempered case)

Regards,

Adip

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/CAL%3DUV%2B9GFpmAkFdu_XU5-WhLRbiGw-Z0f1dXgyavK_gT3LNxJg%40mail.gmail.com.

[Michele Invernizzi]

You are right, there was a missing kbt.

I added it back to the code.

Thanks!

Michele

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/CAMHp%2Bz3cYOkgbNhW0%2BTacPGJ8ya8PqruOujDCo35hLQ9V2Np8Q%40mail.gmail.com.

[Adip Jhaveri]

Hello, I just wanted to know the difference in the algorithms used to calculate the reweighting factor c(t). The expression which plumed uses is - 

There is another expression given in Tiwary and Parinello (2014), https://pubs.acs.org/doi/full/10.1021/jp504920s - 

 

How do these methods differ. And which one should I prefer to calculate c(t) for the non-tempered case?

Thanks,

Adip

On Saturday, July 11, 2020 at 10:58:29 AM UTC-4, Adip Jhaveri wrote:

Hello everybody, thanks for the help!

Federica, thanks for the solution. This solves the infinity problem for me. Although I do have a little doubt regarding the code. For the non-tempered case (biasf_ = -1), the minusBetaF is set to 1 and Z_0 is then calculated. But the kbt_ factor is not multiplied before taking the exponential. Am I missing something?

Pratyush,

While calculating the c(t) you have suggested to use the following expression - 

I was wondering how c(t) would change if I use the expression given in https://pubs.acs.org/doi/full/10.1021/jp504920s - (for non-tempered case)

Regards,

Adip

To unsubscribe from this group and all its topics, send an email to plumed-users+unsubscribe@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/7938859c-7e4e-4eae-909b-4153dba184bcn%40googlegroups.com.

--
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+unsubscribe@googlegroups.com.

[Pratyush Tiwary]

Hi Adip

Within noise they should all give same answer (if not, woohoo, you have a bug or simulation is really not converged! so, different expressions in principle are good for basic debugging)

Otherwise the formula to use is in this review what is marked as eq 22 

https://www.annualreviews.org/doi/abs/10.1146/annurev-physchem-040215-112229

Hope that helps

Pratyush

--

Pratyush Tiwary

Assistant Professor, Department of Chemistry and Biochemistry & Institute for Physical Science and Technology 

University of Maryland, College Park, MD 20742

sites.google.com/site/pratyushtiwary/

Adip

To unsubscribe from this group and all its topics, send an email to plumed-users...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/7938859c-7e4e-4eae-909b-4153dba184bcn%40googlegroups.com.

--
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/CAL%3DUV%2B9GFpmAkFdu_XU5-WhLRbiGw-Z0f1dXgyavK_gT3LNxJg%40mail.gmail.com.

--
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/72204eb6-1fee-4952-8cc0-8736915b855co%40googlegroups.com.

[Adip Jhaveri]

Hi Pratyush,

Thanks for the advice. Is there a way in plumed to calculate c(t) using these different expressions?

Regards,

Adip

On Wednesday, August 5, 2020 at 3:52:05 PM UTC-4, Pratyush Tiwary wrote:

Hi Adip

Within noise they should all give same answer (if not, woohoo, you have a bug or simulation is really not converged! so, different expressions in principle are good for basic debugging)

Otherwise the formula to use is in this review what is marked as eq 22 

https://www.annualreviews.org/doi/abs/10.1146/annurev-physchem-040215-112229

Hope that helps

Pratyush

--

Pratyush Tiwary

Assistant Professor, Department of Chemistry and Biochemistry & Institute for Physical Science and Technology 

University of Maryland, College Park, MD 20742

sites.google.com/site/pratyushtiwary/

Adip

To unsubscribe from this group and all its topics, send an email to plumed...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/7938859c-7e4e-4eae-909b-4153dba184bcn%40googlegroups.com.

--
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...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/CAL%3DUV%2B9GFpmAkFdu_XU5-WhLRbiGw-Z0f1dXgyavK_gT3LNxJg%40mail.gmail.com.

--
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...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/72204eb6-1fee-4952-8cc0-8736915b855co%40googlegroups.com.

[Pratyush Tiwary]

Dear Adip

I recommend going through the source code for PLUMED you are using to see which c(t) is used and connect to the published literature I cited. If you want to use other versions, I recommend writing a code that used HILLS and COLVAR as input (at least that is what I did and found to be a useful exercise). I will be happy to take a look at your code separately if you like.

Pratyush

--

Pratyush Tiwary

Assistant Professor, Department of Chemistry and Biochemistry & Institute for Physical Science and Technology 

University of Maryland, College Park, MD 20742

sites.google.com/site/pratyushtiwary/

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/c4c9a5c0-8a82-4954-9a50-4152bacc2989o%40googlegroups.com.

[jhaver...@gmail.com]

Hi Pratyush,

Yes, I have checked that and then wrote my own code to calculate c(t) using a different expression. But I have a question regarding the plumed reweighting. 

What kind of values are expected for c(t)? I mean, I have read that the values of c(t) is close to to the bias potential, but I am getting larger c(t) which reduces the reweighting factor drastically (exp(V(s,t) - c(t))). Therefore, I am not able to get a reasonable ensemble average after reweighting (kind of end up with zeros since the difference is so large). Is the large difference (Between V(s,t) and c(t)) indicative that the simulation has not converged?  

Regards,

Adip

[Michele Invernizzi]

Hi Adip,

I'm not Pratyush, but hopefully I can be helpful :)

You can find examples of typical c(t) behaviour on the plumed nest, for example here https://www.plumed-nest.org/eggs/19/068/

This is a typical alanine dipeptide metadynamics run, which is the ideal situation, but unfortunately not very representative: https://github.com/invemichele/opes/blob/master/ala2/Colvar-metad.0.data

Personally, in order to better understand the c(t), I found it useful to play around with a simple model where I could calculate the free energy and the c(t) factor exactly, without approximations. If you are interested, here are some of the results https://github.com/invemichele/opes/tree/master/model/metad-model  with both a Colvar.data file from the metad run, and a rct.data with the exact values. You'll find a description of the simulations here https://pubs.acs.org/doi/suppl/10.1021/acs.jpclett.0c00497/suppl_file/jz0c00497_si_001.pdf, in the paragraph "suboptimal double well".

Hope this helps.

Best,

Michele

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/d9bc063e-b192-4b6a-ab83-513dae752606n%40googlegroups.com.

[jhaver...@gmail.com]

Hi Michele,

Help is welcomed from anywhere! Yes, these are some good examples to try and I will do that.

Thanks!

Regards,

Adip