GROMACS version:
GROMACS modification: No
Hi everyone,
Is it possible to logically use the minimum distance between two groups as the pulling coordinate for umbrella sampling? Also, if yes, is it possible to use gmx wham for these simulations?
I ran stationary umbrella sampling simulations with the minimum distance as pulling coordinate between two groups with PLUMED, but now, when I want to calculate the PMF, I get this:
Command line: gmx_mpi wham -ix metadata.dat -o pmf_output.xvg -it tpr-files.dat Fatal error: This is not a tpr with COM pulling
I am not sure now if using minimum distance as the coordinate was right in the first place. I appreciate your help and insight.
Hi!
I haven’t tried running umbrella sampling simulations with PLUMED, but it seems like you’re mixing a couple things up here.
It would be helpful to see your .mdp and plumed.dat files to see exactly where your misunderstanding lies, but in any case, a minimum distance is not a valid pull coordinate. At least I can’t think of a case where it would be. A distance is, though! Think about what you’re trying to do here: you want to calculate the PMF/free energy profile along a given reaction coordinate. How would that work for just one (local) minimum along a distance measure? But even calculating the minimum of your pulling coordinate in your simulation is not meaningful, because that depends entirely on where you positioned your umbrella potentials. A more meaningful quantity would be the distance at which the PMF assumes its minimum, for example.
Regarding the actual error you’re getting, I would think this is because you performed the umbrella sampling with PLUMED, so GROMACS doesn’t know anything about how it’s set up. You don’t need PLUMED to perform umbrella sampling in GROMACS, but if you want to I would suggest looking at one of their tutorials. Here’s a link to a GROMACS tutorial on umbrella sampling.
Thanks a lot. I appreciate your help. My umbrella sampling simulations are stationary and independent. I then want to combine them with WHAM to calculate their PMF.
The reason we considered using the minimum distance was that it would provide more informative insights into the binding properties in my system and yield more accurate windows. However, as our system was very slow, we ultimately decided on the stationary US.
Ah sorry, I misunderstood what you meant by minimum distance, I see how this could be useful now.
If the ZDISTANCES function is differentiable (which I think it should be, looking at the documentation), then this could indeed be meaningful.
However, it seems like gmx wham will only work on tpr files where pulling is active (i.e. pull=yes is specified in the map), as would be the case for a purely GROMACS umbrella sampling simulation. You’ll have to combine your simulation using an implementation of WHAM e.g. in Python.
I think you are right as probably the .tpr is used to get the information of the systems, length of simulations, and position of starting windows?
One thing you can try, @kimiasnn , is to compile those tpr files and pass them to gmx wham. So, if window i of a PLUMED simulation starts at position z_i, use the starting frame to compile the corresponding tpr file with a mdp that has the pull code active with the same biasing potential as PLUMED (umbrella with the same k for the spring) so you pass the information to gmx grompp about the restraint in place and then try to use these tpr files to run gmx wham.
So if I understood right, you mean I can once generate the .tpr files with the umbrella.mdp file using the conf0. But then, what would happen to my minimum distance for my reaction coordinate? Would it generate the right windows?
I am obtaining WHAM from Alan Grossfield’s site. I am testing it now, and I will get back if it works.
In theory the grompp command writes the distances of the groups when you pass an mdp file with the pull code on, so you should be able to verify that the distances are correct. But there are too many ifs, I would say go ahead with the script, as it might be safer in terms of analysis.
Thanks a lot. In what sense do you think using minimum distance would cause degeneracy as a good reaction coordinate?
We were thinking for understanding binding/unbinding to a surface, the minimum distance would be more informative.
It seems to be working. I created a metadata file to calculate the PMF where it has the address to the COLVAR files (time vs min distance), restraint position and spring constant.
Taking the minimum distance to the surface would make sense. I would need to think more about the implications of taking the minimum distance over atoms in the molecule.
But I would think that using the minimum distance as a reaction coordinate in biased simulations would be problematic, as the force only works between the closest atoms and would instantaneously jump from one atom to another when the distances change slightly. This will lead to extreme jump in forces and bad integration.