*Updated*COLVARS PMF extraction from separate simulations. Possible?

Dear LAMMPS forum,
I am setting up simulations to extract Potential of Mean Force using LAMMPS COLVARS. (Lammps version: 3Mar2020)

My setup involves calculating PMF between two alumina particles immersed in aqueous solutions, as a function of interparticle separation. For this purpose I’ll be using WHAM on the output.

My situation:

I need to cover a distance of 25angstroms - as far as I am aware there could be two ways to do this:
a) setup one long simulation with a moving restraint, pulling over the distance in stages, of ‘n’ timesteps each. With stages of ~10ns, this would involve restarting the simulation a number of times in my case.

b) setup individual simulations, using a fixed restraint, at each distance interval (e.g. one simulation at 30, one at 29, 28 angstroms etc). Then somehow combine the separate traj files in post-processing from the separate simulations.

My question:

can PMF be obtained using the second version? If yes, is outputAccumulatedWork or outputAppliedForce the appropriate option to use? This latter bit has been confusing me for a while since the COLVARS manual (page 86) talks about outputAccumulatedWork being used to extract PMF (for moving constraints) but no mention of outputAppliedForce in the context of PMF that I could find. I have been advised to use outputAppliedForce, but I’m not sure about it.

If anyone, please, could answer these questions clearly/explain it would really help.

Please find my colvars file attached - I’ve managed continuous pulling over the distance to generate starting configurations, but am running into problems with what options to implement for holding the interparticle separation at a fixed distance (presuming this version can be used to obtain PMF).

Kind regards,

PMF.colvars (611 Bytes)

Hi Olivera, it depends on what implementation of WHAM you are using. When you process only the trajectory of the variables, you need to provide the Hamiltonian of each window to unbias its computed probability distribution. This is also true if you just accumulate the histogram yourself and give it to WHAM.

Some WHAM implementations will require you to type in the restraint’s center and force constant. Others will used the matching trajectory of applied forces to infer the restraint potential. Either one is correct, so if you turn on outputAppliedForce you are definitely leaving yourself open to all.

Regarding (a) vs. (b), the first should give you better-equilibrated snapshots for the windows, but it would take longer because you would need to run them consecutively. Either way, WHAM won’t care if the windows are consecutive or are run in parallel.

Umbrella sampling is a fixed-restraint method (time-independent), so you don’t need targetCenters and targetNumSteps. You may have seen tutorials where they run WHAM over a steered MD (time-dependent) but that’s quite poor in my opinion, because in that case you would not only make the approximation that your finite histogram represents a probability, but also that that probability would be steady (when it’s not).


Hello Giacomo,
thank you, that has really helped :slightly_smiling_face:

So, am I right in thinking:

-> specifying outputAppliedForce, Centers and ForceConstant in the colvars file with version b (separate simulations for each distance) would be a sufficient combination, for at least one implementation of WHAM, to extract the PMF?
Kind regards,

That would work. Take care in setting up the initial structure for each window as accurately as possible to reduce non-equilibrium effects.

You may also want to take a look at the writeTIPMF flag, because distanceZ has TI gradients available:
It is an additional flag that has no effect on dynamics, but it will collect an PMF based on thermodynamic integration that you can compare the probability-derived PMF from WHAM. The PMF is for each window individually, but they may be combined together relatively easily. A PMF computed in this way requires lowerBoundary, upperBoundary and width specified in the config as well.


Hello Giacomo,
thank you for the suggestion - I should probably try that as a double check.

However, it’s raised another question I’ve had for a while so I may as well ask now:
what’s the physical meaning of ‘width’ and is it necessary in my case to change it from the default value?
Say, I want to restrain at 30angstrom separation in z, would ‘width’ of 1 be permitting a 1angstrom fluctuation either side of that? I don’t quite understand the ‘scaling’ aspect of it. I guess I also have the same question for lower/upperBoundary.

Any wisdom much appreciated :slightly_smiling_face:

Kind regards,

Hi Olivera, with respect to restraints “width” allows you to specify a consistent scale in case you apply multi-dimensional restraints that share a force constant. It also helps you compute easily how much force constant you need to set a standard deviation equal to the width itself.

For a 1D harmonic oscillator with mean = 0 and SD sigma = 1 at 300K and in real units: ~= 1/2 K <sigma^2> = 1/2 K, so to set = kBT you can simply use K = 2kBT. Thus, because the harmonic oscillator is implemented for the variable scaled by “width”, this relationship will always hold.