Post-processing process of WHAM Umbrella Sampling

[高奋娥]

Dear developers and users of Plumed,

I'm a new user of Plumed. I would like to do an ab initio simulation using PLUMED2.5 interfaced to the CP2K MD engine refer to the Exercise 4 in the tutorial: https://www.plumed.org/doc-v2.5/user-doc/html/lugano-2.html.

It can be successfully calculated using this bash script run_us.sh:

for AT in -3.00  -2.75 -2.50 -2.25 -2.00 -1.75 -1.50 -1.25 -1.00 -0.75 -0.50 -0.25 +0.00 +0.25 +0.50 +0.75 +1.00 +1.25 +1.50 +1.75 +2.00 +2.25 +2.50 +2.75 +3.00

do

cat >plumed$AT.dat << EOF

# vim:ft=plumed

phi: TORSION ATOMS=5,7,9,15

psi: TORSION ATOMS=7,9,15,17

# Impose an umbrella potential on CV 1 and CV 2 with a spring constant 

# of 500 kjoule/mol at fixed points on the Ramachandran plot

restraint-phi: RESTRAINT ARG=phi KAPPA=250.0 AT=$AT

# monitor the two variables and the bias potential from the two restraints

PRINT STRIDE=10 ARG=phi,psi,restraint-phi.bias FILE=COLVAR$AT

FLUSH STRIDE=10

EOF

sed -e '2s/diala/diala'$AT'/g' -e '132s/plumed/plumed'$AT'/g' cp2k.inp > cp2k$AT.inp

sed '24s/'cp2k'/cp2k'$AT'/g' cp2k.pbs > cp2k$AT.pbs

qsub cp2k$AT.pbs

done

However, I have some problems with the post-processing process about the number 2.49. In the tutorial to run the iterative WHAM optimization and get a weight per frame we need to implement: python wham.py biases.dat 25 2.49. where 2.49 is the temperature in energy unit. The free energy profile along phi can be obtained after the following two step. 

1. paste allphi.dat weights.dat | grep -v \# > allphi-w.dat.  

2. python do_fes.py allphi-w.dat 1 -3.1415 3.1415 50 2.49 fes.dat.

If the system is simulated at 310 K, 400K or 600K(set TEMPERATURE 400 or 600.0 in the MD module in the input file of cp2k.),  Is it even okay to use 2.49?

I hope someone can point me in the right direction here as I'm a bit confused.

Thanks in advance, 

Fen‘e Gao 

[Shaunak Badani]

Hi,

I believe that the value of kB in kJ / (mol K) is 0.0083144621 = (~1.38 * 10^(-23) * 6.023 * 10^(23) * 10 ^ (-3)). So when you multiply this value by a temperature of 300, you get ~2.4 as kBT, which is the input, I'm assuming. 
So if you were to use a temperature of 600 K, the value would just change to 4.8 .

I believe this should be the answer.

Kind Regards,

Shaunak Badani,

Undergraduate Researcher in CCNSB, IIIT Hyderabad.