GROMACS version: 2020.6
GROMACS modification: Yes
Dear,
I am doing SMD using collar, and I want to obtain the force (in pN) and the accumulated work and know how the work is calculated (to quote the equation in my article),
I also assume that the work will have units of kJ / mol, right?
best,
Geo.
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.
Giacomo
Professor,
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
indexFile segid_whole_index.ndx
#! Restrain in cylinder during pulling.
#! objetive is not permit that atoms out of box.
colvar {
name restraint_xy
width 1.0
lowerboundary 0.0
upperboundary 1.0
distanceXY {
main { indexGroup PROB_CA }
ref { dummyAtom (0.0, 0.0, 0.0) }
}
}
harmonicWalls {
name hw0
colvars restraint_xy
lowerWalls 4.4
upperWalls 8.4 #! to restrain into an cilynder of 4.0 nm of diameter.
lowerWallConstant 1000 #! kJ/mol x nm^2
upperWallConstant 1000
}
#! Pulling Colvar & Protocol
#! pulling from z+ to z-
#! velocy (v): v = (targetCenters - centers)/(targetNumSteps * dt)
#! where dt = 0.002 fs
colvar {
name move_to_z
width 1.0
outputAppliedForce yes
distanceZ {
main { indexGroup PROB_CA }
ref { dummyAtom (0.0, 0.0, 0.0) }
}
}
harmonic {
name harm0
colvars move_to_z
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
targetNumSteps 22500000
outputCenters on
outputAccumulatedWork yes
}
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:
Giacomo
Prof. Giacomo,
Thanks so much for it, now it clear all.
Best,
Geo.
Professor,
I want to know how I can use some enhanced simple method by Gromacs+colvars, such as well-tempered metadynamics, well-tempered eabf and so on.
Well-tempered metadynamics is supported (there is a specific keyword for that, check the doc).
I’m not aware of a well-tempered eABF method, but if you mean the concurrent use of eABF with WT-metadynamics, that is of course also supported. See the doc for details.
Giacomo
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.
What was the problem with index groups? Have you mentioned it to me before? Since I don’t know your name, I don’t know if we have talked.
Jérôme
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.
Giacomo
Indeed, as Prof. Giacomo mentioned, I was the one who solved that question, which thanks to Prof. Giacomo’s suggestion I was able to solve.
Here I attach a python script to make the index list in a correct way to work with colvar, maybe it will serve for other users as a template for their own proposals.
best,
Geo.
(Attachment make_index.py is missing)
Just paste the relevant code snippet for what you need to do when generating the index file in MDAnalysis.
You’ll need that same syntax for generating a correct index file for use in GROMACS itself, which sets the standard by not allowing spaces in groups’ names.
Giacomo
Hi Prof Giacomo,
I’ve view strange thing into my colvars taj file:
22500000 1.86979699134827e+00 -8.93909740447996e+01 1.84000000000000e+00 7.31263625161090e+02
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?
Geo.
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.
Giacomo
I did multiple reps of 200 to be exact. and in all of them I have final accumulated values > 700 kJ , I am using k = 1000 kJ/mol nm and pulling v = 1 nm/ns.
How can I fix it?
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.
Giacomo
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 ?
Geo.
The unit conversion looks okay as far as I can tell. Whether the protocol itself is appropriate for the problem should be evaluated by somebody working together with you on this project.
Giacomo
Professor Giacomo,
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.
Giacomo