[gmx-users] pull code, constraint force

Chinmay Das chinmaydas71 at googlemail.com
Thu Nov 23 11:18:53 CET 2006

I am trying to calculate free energy and diffusivity of a guest
molecule at a given height from lipid bilayer. I am using pull code
with runtype=constraint, reftype=com, pulldim=N N Y
in version 3.3.1 compiled in double precision. The ensemble is NPT
(Nose-Hoover/ Parrinello-Rahman).

To save on storage, I included extra lines in src/kernel/md.c to
calculate displacements of the guest molecule and force felt by the
guest molecule. This is included just before updating
    xx = (do_per_step(step,inputrec->nstxout) ...

To get the force components, I am summing the force for all atoms in
the molecule - which should cancel out intra-molecule contributions.

The z component of the force is about the same as the output in
pull.pdo - but not quite so. The differences are in the third decimal
place - which is too large for double precision.

Am I missing something and the constraint force is not same as the
force felt by the molecule? Or the difference is because I am
calculating the force at the beginning of the time step, while the
pull code, as far as I can make out, uses iterative scheme to fix the
constraint as part of the updating scheme? Or errors are accumulating
in the iterative scheme?

Finally this is probably not at all important for calculating
diffusion constant (from force auto-correlation) or free-energy
derivative (average z force) because the difference is small - but a
personal curiosity about what is going on and to find out if I am
missing some subtle point.

With best regards

More information about the gromacs.org_gmx-users mailing list