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