PMF profile from PLUMED vs gmx wham (are not same)

[srdjan125]

Hello, I did my US simulations only with GROMACS and now I am learning PLUMED.

I have concatenated trajectory of my US simulations and used following command: (trajectory has 9000 snapshots)

plumed driver --noatoms --plumed plumed-wham.dat --mc mcfile --timestep 20 --trajectory-stride 3000

plumed-wham.dat file has following format:

Dim1: GROUP NDX_FILE=index.ndx NDX_GROUP=DimerA
Dim2: GROUP NDX_FILE=index.ndx NDX_GROUP=DimerB
com1: COM ATOMS=Dim1
com2: COM ATOMS=Dim2
r: DISTANCE ATOMS=com2,com1
RESTRAINT ARG=r KAPPA=1500.0 AT=3.933
RESTRAINT ARG=r KAPPA=1500.0 AT=3.995
.....
PRINT ARG=*.bias FILE=biases.dat STRIDE=1
PRINT ARG=r FILE=All_COMs.dat STRIDE=1

Upon that I have used wham program (python wham.py biases.dat 54 2.49) and first I thought that I have got 0 for each snaphots, but later I realized that python wrote them like 0, so I have wrote them like weights (return np.log(w) ).

Then I did FES by using following script:

**_com1: READ FILE=bias_wham2.dat VALUES=com1 IGNORE_TIME
lw: READ FILE=bias_wham2.dat VALUES=logweights IGNORE_TIME

we use a smooth kernel to produce a nicer graph here

hhphir: HISTOGRAM ARG=com1 GRID_MIN=3.0 GRID_MAX=8.5 GRID_BIN=100 BANDWIDTH=0.1 LOGWEIGHTS=lw
ffphir: CONVERT_TO_FES GRID=hhphir TEMP=300
DUMPGRID GRID=ffphir FILE=FES_unbiased.dat_**

and I got some free energy profile which kinda resembles that free energy from gmx wham, but minimum is twice lower compared to gromacs. I don't know did I make some mistake or WHAM was evaluated only on 9000 snapshots while GROMACS evaluated it on 54*100000 snapshots.

My bias_wham2.dat files look like this:
! FIELDS time com1 b0 logweights
0.000000000000000000e+00 3.874699999999999811e+00 -9.505879999999999864e-05 2.516569416381209976e+02
0.000000000000000000e+00 3.874699999999999811e+00 -9.505879999999999864e-05 2.508529784958790003e+02
....

If I supply logweights as negative numbers, FES is totally weird (like some triangle) and totally shifted from real minimum.

I would appreciate if someone can help me with advice if I am doing something wrong here.

I also tried to compute bias by another approach, just to make plumed driver read COM distances from GROMACS generated pullx files (I used python to make COLVAR.dat file which plumed can read). Now I am trying to repeat the same by taking more snapshots from US (other than 9000 ones which I have saved in XTC file), and see if free energy profile will resemble the one from gmx wham. I guess I can use this too because plumed driver will either compute COM from trajectory, or it will read COM from pullx files? Or this is a wrong way to do it?

**_#SETTINGS MOLFILE=user-doc/tutorials/lugano-2/diala.pdb

vim:ft=plumed

r: READ FILE=COLVAR.dat VALUES=com1 IGNORE_TIME IGNORE_FORCES
RESTRAINT ARG=r KAPPA=1500.0 AT=3.933
...._**

plumed driver --noatoms --plumed plumed-wham.dat