# Colvar output force and work

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:

1. 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.
2. 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.

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