sum_hills multiple HILLS from independent simulation

[wu bohua]

Dear supporters，

I am running 40 independent well-tempered Metad simulations with the same plumed input file. The difference of all 40 independent is only the initial random velocity in the simulation setting input in .mdin file in AMBER software. After simulation, I got multiple HILLS file. 

I process these 40 HILLS like this way:

plumed sum_hills --hills 1/HILLS,2/HILLS ,3/HILLS,....,40/HILLS --idw cv1,cv2 --kt 2.53 --mintozero --outfile cv1-cv2.dat

When I checked the cv1-cv2.dat, I find the values in the third column are many times larger than those calculated with a separate file. I am wondering if I want to get the correct calculation, what should I do?

The results from the single HILLS:

   -2.873840000   -2.816270000   35.874444239
   -2.768320370   -2.816270000   35.874444057
   -2.662800741   -2.816270000   35.874443455
   -2.557281111   -2.816270000   35.874441625
   -2.451761481   -2.816270000   35.874436651
   -2.346241852   -2.816270000   35.874424798

   ......

The results from the 40 HILLS:

   -2.925840000   -2.856170000 1349.727101088
   -2.820087455   -2.856170000 1349.727100444
   -2.714334909   -2.856170000 1349.727097968
   -2.608582364   -2.856170000 1349.727090117
   -2.502829818   -2.856170000 1349.727067416
   -2.397077273   -2.856170000 1349.727009534
    ......

Best regards,

Bohua

[alx.th...@gmail.com]

Hi Bohua,

What you are doing here is you are simply summing all individual potentials. So it is no surprising that energy in this summed profile is ~40 times the energy in the single profile.

I suppose that what you really want to do is to produce an average profile, and you made 40 replicas to report errorbars. It is tricky though because you first need to understand how to align these profiles. If you just make the global minimum zero in each separate profile and then take average, you will have apparent 0 error in this region and high errors in remote regions. Are these errors actually meaningful? Well, no. So, your aligning scheme influences what information your errors confer.

What is an easy way to overcome it is to perform Tiwary reweighting since profiles it produces come naturally aligned. Then, you can simply take average and std or sem of all your reweighted profiles.

sum_hills on many HILLS files actually works when these HILLS come from dependent replicas. For example, when using file-based multiple walkers scheme.

Hope it is helpful,

Alex

[wu bohua]

Dear Alex,

Thanks for your quick and kind reply. Based on your suggestions, I check the plumed tutorial to find the Tiwary reweighting method, but I can only find the umbrella-sampling like method in this one. And I checked his paper (A Time-Independent Free Energy Estimator for Metadynamics.J. Phys. Chem. B 2015, 119, 3, 736–742). This paper illustrated theories about how to derive a time-independent and locally convergent free energy estimator for metadynamics. It is still hard for me to deal with my case. Could you  give me some suggestions how to process these HILLs files with Tiwary reweighting method.

Best regards,

Bohua

[alx.th...@gmail.com]

Dear Bohua,

Reweighting is performed by histograming the information from the trajectory itself, not HILLS file. In order to do it, you need to have a weight assigned to each time point. In order to assign a weight, you need to calculate rct factors. Then you subtract them from your metadynamics bias, so you have reduced bias. You can then perform histograming by reading the column with reduced bias values altogether with the column for the data that will serve as the space in which your profile will be built. 

However, I do not know how to calculate rct factors aposteriori (though there may be some way, I just never needed it). What I do is that I calculate rct and reduced biases (rbias) in place, while running metadynamics simulations. It is activated by CALC_RCT keyword. Then, if your label for MetaD data is, say, M, you can output PRINT FILE=CV ARG=CV1,CV2,M.rbias STRIDE=1000 for example.

And then you write your reweight.dat file:

CV1: READ FILE=CV IGNORE_TIME IGNORE_FORCES VALUES=CV1

CV2: READ FILE=CV IGNORE_TIME IGNORE_FORCES VALUES=CV2

RW: READ FILE=CV IGNORE_TIME IGNORE_FORCES VALUES=M.rbias

weights: REWEIGHT_BIAS ARG=RW.rbias TEMP=300

HISTOGRAM ... 

  ARG=CV1,CV2

  GRID_MIN=-0.2,-0.2 (adjust for your space)

  GRID_MAX=0.2,0.2

  GRID_SPACING=0.01,0.01 (adjust)

  BANDWIDTH=0.01,0.01 (adjust)

  LOGWEIGHTS=weights

  LABEL=h1

... HISTOGRAM

ff1: CONVERT_TO_FES GRID=h1 TEMP=300                                                                                    

DUMPGRID GRID=ff1 FILE=fes.dat                 

And you get your reweighted free energy profile. Please also note that you need to throw away some initial part of your CV file ("initial transient"). The question of how do you know which part to throw away is tricky. One way that was proposed by Michele Invernizzi I suppose was to look at M.rct plotted against time and note the time point after which the graph becomes roughly linear (not flat, but just rising steadily). Than you throw away everything before this time point. You can see what errors you can get by retaining this data by simulating a toy symmetrical system (say ethane eclipsed/staggered switch, with torsion as CV). You know a priori that the truth is 3 equal minima, but if you perform reweighting on the full trajectory, the free energy minimum for the state you start your simulation from will be lower. It happens because these calculations (rct, rbias) are approximate, and a degree of this approximation is higher at the start of the trajectory, so you get weights calculated with much more uncertainty.

But, even more easy way to solve your troubles is to switch to OPES (OPES_METAD or OPES_METAD_EXPLORE). It works similarly to metadynamics but is easier to setup, and it is much easier to compare profiles and derive errorbars from it. For it, you need to compile Plumed 2.8.x with OPES module enabled at the time of configuration. These methods output "already reweighted" biases (this is technically incorrect, but works as an analogy to what was descibed for metad), and histograming as I wrote higher is a default way to build profiles in these methods (you just need to replace metad.rbias for opes.bias). This way you can run 40 simulations, get 40 profiles and compute average profile and your uncertainty in calling this average profile the one true profile with minimal effort.

Best. 

Alex

[wu bohua]

Dear Alex,

Thanks for your quick and kind reply. And I will attempt to use the reweighting  procedure as  you said in my case. And could you please tell me which reference this method  is from？Because maybe I need to cite this literature in my paper. 

Best regards,
Bohua