Colvars module: details on the pmf output computation

LAMMPS LAMMPS General Discussion
,
1

Hi,
I am using the colvars module with lammps to perform steered molecular dynamic.
Thank you very much the developper to provide us such tools.

The method works well in my case, but I would like to be sure of the formula used calculated the pmf. After defining my CV, I perform the steered MD using a harmonic restraint at position \xi_0. In order to check if my calculation is correct, I do a cycle on this harmonic restrain \xi_0 varying first from -0.75 to 0.98 (file atom-up.colvars) and then from 0.98 to -0.75 (file atom-down.colvars).
Both increase and decrease of \xi_0 are done during the same lammps calculation. Below, you will find the main instructions I used.

My question: What is the formula used to calculate the pmf at the end of my simulation ( down.colvars.pmf)? I ask this question since I have seen that the down.colvars.pmf involves some statistic of the whole cycle.
In addition, the up.colvars.pmf seems to have been smoothed compared to the direct plot of the accumulated work as a fonction of \xi_0.

Doing some bibliography, I have seen that the pmf could be computed using

  • the Jarzinsky equality (average over several trajetories)
  • the evaluation of the reversible work (and thus pmf) from irreversible works of cyclic paths (J. Chem. Phys. 122, 104106 (2005) )

Could please tell him how the pmf is evaluated in colvars (I did not find the precise answer in the manual file)?

Thank you very much in advance,
Nicolas

Lammps input file

  fix 14 gcell colvars atom-up.colvars output up  tstat 16
  run  7000000
  unfix 14
  fix 14 gcell colvars atom-down.colvars input up.colvars.state output first tstat 16 
  run  7000000

atom-up.colvars file:
harmonic {
name har
colvars cv1
forceConstant 20
centers -0.75
targetCenters 0.98
targetNumSteps 7000000
outputCenters yes
outputAccumulatedWork yes
writeTIPMF yes
}

atom-down.colvars file:
harmonic {
name har
colvars cv1
forceConstant 20
centers 2.71
targetCenters -0.75
targetNumSteps 14000000 (-> since the timestep is not reset, atom-down.colvars is called at timestep 7000000…which will correspond to \xi_0=0.98)
outputCenters yes
outputAccumulatedWork yes
writeTIPMF yes
}

2

You’re welcome (also on behalf of everybody else who wrote it).

With SMD, you typically get non-equilibrium work, and how well that approximates the PMF along that variable really depends on the system. Typically, people would use other methods, like umbrella sampling, metadynamics, ABF… The latter two will produce a .pmf file directly, using what is documented in the respective sections. Umbrella sampling requires post-processing using WHAM, MBAR, or similar tools.

Additionally (and this is more unique to Colvars), the thermodynamic integration free-energy estimator that ABF uses can be used in combination with some other methods, and produce a .pmf. Here is a link for other readers’ convenience, but it can be easily found it by looking up the writeTIPMF keyword that you have in your config file :slight_smile:

Giacomo