Project

General

Profile

Bug #73

distance restaints problem when using several molecules

Added by Jochen Hub over 13 years ago. Updated about 13 years ago.

Status:
Closed
Priority:
Normal
Assignee:
Category:
mdrun
Target version:
Affected version - extra info:
Affected version:
Difficulty:
uncategorized
Close

Description

When using distance restraits in several identical molecules with 'pair
restraits' that belong to the 'same restraint' (i.e. they have the same "index"
in the itp file) grompp makes a mistake.

What grompp does, is that he handels the restraints in the first molecule
correctly but for all following molecules he handels each 'pair restraints' as
an independent restrain instead of putting it into the same restraint. This can
be observed when using gmxdump -s topol.tpr | grep DISRE. You will see that the
pair restraints in the first molecule belong to the same functiontype while all
pair restrains the folloing molecules have their own functiontype.

The problem is caused in the bondcat routine in topcat.c. See:

for (j=0; (j<copies); j++) {
for (i=0; (i<src->nr); i++) {
memcpy((char *)&(dest->param[l]),
(char *)&(src->param[i]),(size_t)sizeof(src->param[i]));
for (m=0; (m<nral); m++)
dest->param[l].a[m] = n0;
if (bDisres && (j>0)) {
>>>>>>> max_index
+; <<<<<<<<<
^^^^^^^^^^
dest->param[l].c[0] = max_index;
}
l++;
}
n0 += dstart;
}

The param[l].c0 is later put into the "label" of the disres functiontype
variable in convparm.c

case F_DISRES:
new->disres.label = old[0]; (old[0] is c[0] in bondcat)

and if the labels of two pair restraints are different they are treated
independetly. So probalby one
should better write in the boundcat routine:

if (bDisres && (j>0)) {
/* max_index++; */
dest->param[l].c[0] = max_index + j;
}

Alternatively the max_index can be incremented inside the the j-Loop but NOT
inside the i-Loop.

In my case this way thing got correct. Please let me know if you would like to
have sample input files (pdb, itp, top and ndx) that reproduce the error.

Cheers from Goettingen,
Jochen

History

#1 Updated by Berk Hess over 13 years ago

Sorry, I only noticed this bug now, but better late then never.

I would suggest a different solution.
The current way renumbers everything, which makes things
confusing when one started with a non sequential index numbering.

So I suggest to add j*(max_index + 1) to each restraint
that should not be ensemble averaged.

#2 Updated by David Mobley about 13 years ago

OK, so if what you're saying is correct, this really comes back to the problem that, to apply restraints
between two molecules, they must be defined to be part of the same "molecule" in GROMACS, yes? Only
it is slightly worse than that, in that they only get treated as one molecule IF I apply restraints
between the two, but not otherwise, even though they are in the same molecule description.

I suppose I agree that this is not exactly a bug, but it seems that at least this stuff ought to be put on a
list of "improvements to be made" or something, as it seems that it's far from ideal.

I agree the topologies are not identical -- but it is EXTREMELY nonintuitive that the definition of, say,
an angle restraint in a topology should cause changes in computed dihedrals, don't you agree?

Maybe another work-around would be to run ALL of my simulations with a restraint applied between
the ligand and protein, using a spring constant of 0? Then the protein and ligand will be treated as the
same molecule consistently, right? That is, I can run "unrestrained" simulations using a spring constant
of 0, and restrained simulations by turning up the spring constant, and at least the treatment of the
dihedrals will be consistent across the two simulations?

I've been aware in the past that PBC can influence dihedrals, so usually what I do is try a couple choices
of reference atoms for choosing my dihedral restraints, and I pick the set that gives me a nice, stable
dihedral angle (indicating that PBC isn't causing flipping in the dihedral angle). I guess I just never
realized before that whether or not restraints are used influences how PBC is handled.

In summary:
1) I think I know how to proceed now, either by (a) always using a "nominal" restraint with a spring
constant of 0 to make sure the treatment of PBC is consistent, or (b) by choosing reference atoms
closer to the ligand.
2) I think this situation is far from ideal, as it is EXTREMELY nonintuitive that the use or non-use of
restraints should influence computed dihedrals. Maybe it isn't a bug, but it's something that should be
changed eventually.

#3 Updated by David Mobley about 13 years ago

Whoops, sorry, I just accidentally committed a comment that had nothing to do with this bug (pertained to
bug #100). Apologies.

#4 Updated by Berk Hess about 13 years ago

Fixed the bug as described in my previous comment.

Also available in: Atom PDF