Hello Geo, yes, the units are GROMACS units so kJ/mol for energy/work and kJ/mol/nm for the forces. Note also that you can use outputAppliedForce to report that force during that simulation:
and you can verify that it matches what you expect.
Assuming that you’re using a moving harmonic restraint, the “SMD formula” would just have the center of the restraint linearly interpolated between the initial and target value.
With respect to the GROMACS’ pull code, the restraint would not keep moving indefinitely with a constant velocity, but only reach the value assigned by targetCenters, and then stay there. It’s a good idea to use outputCenters for the moving restraint, so that you can see it going (but it’s almost trivial to calculate anyhow).
If you use the outputAccumulatedWork the integral of the restraint force times the CV’s net displacement is computed at every MD step, and printed periodically.
this is my code (bellow), at this point I have an additional question, I want to apply a restriction on a cylinder of radius 20 and diameter 40, I want to know if I correctly defined the collective variable and it values of lower and upper boundary.
colvarsTrajFrequency 100 #colvarTrajAppend off #analysis off
#! Restrain in cylinder during pulling.
#! objetive is not permit that atoms out of box.
centers 6.4 #! in nm com
targetCenters 1.9 #! in nm
forceConstant 2100 #! in kj/mol nm^2 → ~2.5 kcal/mol x A^2
Hi Geo, this is a bit confusing because you are giving bits of information, but asking to validate a whole input. It is you and your direct collaborators who know the purpose of the computation, and it’s your job to divide them up into specific questions: for most of these you’ll find the answer in the doc. For a few of these specifically, a question on the forum may help.
The distanceXY function computes the distance from the axis that passes through the center of mass of the ref group. The direction of that axis is by default parallel to the Z axis. For cylindrical symmetry, the lowest possible value of distanceXY is zero, when main and ref both lie along the axis. If you want to confine main to be within a cylinder of given radius, use that radius as upperWalls for a walls restraint.
Most importantly, there are two things you need to be careful about with GROMACS:
Please use either Å or nm consistently, be explicit about your choice when you are asking for help and please only use nm in GROMACS, as that’s the unit that it uses.
I already mentioned that the unit cell for GROMACS has a corner on the origin, but its center is (Lx/2, Ly/2, Lz/2). I don’t know what’s in your unit cell but you may have to adjust that setting.
Okay, thank you very much, Professor. I just mean the varianst of eABF method(meta-eABF and well-tempered-meta-eABF). Previously, I have used the Gromacs-2018.8+colvars. However, I think there are some bugs during chosing the index atoms in that version. And I I gave the question back to Dr. Jérôme Hénin.
For the record, in case if you refer to this issue with index group names:
that issue specifically was due to having used spaces inside a group’s name.
The GROMACS index group format (as produced by e.g. make_ndx) does not allow for spaces. The main difference between GROMACS and Colvars when a space is there is that GROMACS will just read everything up to the first space to use as the group’s name and discard the rest, whereas Colvars will raise a parsing error.
this is output from running 4.5 ns de md ( 2 fs) y 1 nm/ns velocity pulling and colvarsTrajFrequency 100. but the outputAccumulatedWork yes , givenme work on ~730 kJ/mol, and when I try calculate PMF using jarzynski using 2do order cummulants (PMF = - Beta/2 [<w^2> - ^2] + …) but ussing data from outputAccumulatedWork yes get high values ~3000 kJ/mol. how can I get rigth?
If this is a biomolecular system, your simulation is probably out of equilibrium (because of the >700 kJ/mol of accumulated work).
Also, please correct me if I’m wrong, but shouldn’t the cumulant expansion require more than one trajectory? Jarzynski’s equality is valid for a single trajectory (trivially so), but to compute the variance of the work you should need multiple samples.
Ok, if you have multiple replicas with the same result, then getting the variance makes sense. But at the same time if all replicas give very similar result, the variance will not matter much.
Going back to the first part of my answer, there are probably some unphysical configurations being visited consistently. Take any two replicas and analyze the atomistic trajectories, and cross-reference what you see on the screen with when does the collective variable begin lagging behind the restraint. That may inform a better choice of variable.
Dear Prof, change a bit the pulling atoms etc. Now I have a maximum force of 316 kJ/mol nm , so I would like to calculate the work myself, from what I see the work is basically: w = force * velocity * dt
here my speed was 1 A/ns which is equivalent to 0.1 nm/nm, the frequency of the print data in colvars was 100 steps, this is equivalent to 0.0002 ns (2 fs timestep), these conversions are correct, to calculate the work ?
I am reviewing literature, and I am trying to change the protocol for Umbrella sampling, because I am not getting improvements with jarzynski, I have corrected the pulling, changed forces, atoms to pull, etc, but the work results are always similar, very big.
My question is, I saw that with COLVARS, umbrella sampling could be done directly in namd/gromacs, could you tell me how I could add these parameters to my colvars file?
Generally speaking, umbrella sampling may involve any kind of steady-state restraint that does not change over time. The most popular example is a static harmonic restraint, which you can obtain from your existing input by removing targetCenters and instead setting multiple values of centers to cover the landscape. The theory behind how to choose the spacing with US is the same as any other implementation (i.e. you need overlap between windows).
To obtain a PMF, the probability distribution of the resulting energy landscape (including the restraint) is sampled over time, and is then unbiased using a post-processing method. WHAM and M-BAR are two very widely known examples of this. One of my colleagues has also developed recently another method that uses applied forces.