Project

General

Profile

Bug #100

interference between restraint definition and internal dihedral angles

Added by David Mobley about 13 years ago. Updated about 13 years ago.

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

Description

I'm attaching a tarball to reproduce a very crazy bug I've been running into.
Basically, it seems that the definition of certain restraints in a topology can
influence computed dihedral angles. First I'll describe how to reproduce the
problem, then I'll explain more of what I know about it.

In the attached tarball, find "wkg.tpr" and "notwkg.tpr", which were both
generated by running "grompp -f min_lbfgs0.mdp -c complex_restr.gro -o (name) -p
complex_restr.top", the only difference was for "notwkg.tpr" I uncommented one
of the angle restraint lines in the topology, and for "wkg.tpr" I left all of
the restraint lines commented out. Uncommenting any of the restraint lines
seems to cause the same problem.

The basic problem, then, is that the computed dihedral angle from the .gro file
depends on which tpr file I use, which it shouldn't, since the only difference
for generating the tpr files was whether restraints were defined in the topology
file or not. That is, if I run "g_angle -type dihedral -n anglegps.ndx -f
complex_restr.gro -s (name)" and enter "4" (or select any of the groups named
"phi", actually), I get two different answers depending on which tpr file I use.

This problem, however, extends beyond just g_angle. That is, for this particular
system, whenever I run with restraints defined in my topology file, I get
huge dihedral restraint energies output to g_energy. Also, g_angle gives me
wildly different dihedral angles than if I start a run from the same starting
structure without restraints defined, even beginning from the very first frame.

My best guess is that somehow the definition of restraints is messing up the
atom indexing or something internal that affects how dihedral angles are
computed for at least the atoms I'm interested in, which (a) leads to big
restraint energies, and (b) shows up in differences in output dihedral angles.

Please take a look ASAP, as (as usual) I'm stalled until this is fixed.

bug.tar.gz (3.57 MB) bug.tar.gz Tarball to reproduce the problem David Mobley, 08/22/2006 09:33 PM
tprs.tar.gz (2.59 MB) tprs.tar.gz tpr file with pbc full with and without restraints in topology David Mobley, 08/23/2006 05:46 PM

History

#1 Updated by David Mobley about 13 years ago

Created an attachment (id=64)
Tarball to reproduce the problem

As described in the bugzilla

#2 Updated by David van der Spoel about 13 years ago

I have reprodcued the problem. Do you know which angle is correct? I presume the
one without angle restraint?

#3 Updated by David Mobley about 13 years ago

Yes, I believe the one without the angle restraint.

This ends up causing a bunch of my simulations to "blow up": I run a simulation with no restraints to
calculate the preferred values of a set of angles and dihedral angles, then I turn on restraints to keep
those angles and dihedrals at their preferred values. But then I get huge dihedral restraint energies, since
the computed dihedrals no longer match with what I get without the restraints turned on (even though the
structure is exactly the same)...

Thanks,
David

#4 Updated by David van der Spoel about 13 years ago

This is a periodicity problem. The coordinates are shifted differently when the
angle restraints are present. While we try to fix it, please try to use pbc=full
as a workaround.

#5 Updated by David van der Spoel about 13 years ago

OK, after some mailing with Berk and Erik, we conclude that your problem is
ill-defined, because the peptide is half a box away from the protein. That means
that is somehow arbitrary which periodic image is being used to compute the
interaction (pbc = full won't help either!). A very small displacement of one of
the atoms (leading to using another image) can henceforth lead to very different
angles because there are intermolecular vectors involved. If you would define
both vectors inside one molecule each it would work.

Try to reduce the intermolecular distance with respect to the box size.

#6 Updated by David Mobley about 13 years ago

The peptide isn't actually half a box away from the protein -- at least, if I use trjconv with "-pbc inbox"
to make a pdb of it and look at this in a viewer, the peptide is still bound in the binding site on the
protein (which is what I want).

However, if I take you to mean "the distance from the peptide [ligand] to the atoms in the protein is too
large", you might be right. The method I'm using requires a choice of three arbitrary atoms in the
protein to use for reference for restraints, so I will try again using atoms nearer the ligand.

This does, however, still seem like a bug to me: If I run g_angle to compute the dihedral angles I'm
interested on a trajectory run without restraints, I CONSISTENTLY (for every snapshot!) get dihedral
angles of a certain value (near 80 for phi_A, if I remember correctly), and if I run g_angle on a trajectory
run with restraints (or presumably, on the SAME trajectory but using a tpr file that has restraints
defined in it, although I haven't tried this yet, but I can if you want) I CONSISTENTLY get an entirely
different set of dihedral angles, even though the ligand is still in basically the same place. I can buy it is
a periodicity problem if you tell me that it means the dihedral angles are well defined, but I can't buy
it's entirely a periodicity problem if you tell me that it means that the computed dihedral angles depend
on whether I use restraints in my topology or not.

What if I ran a trajectory of this and wanted to use mdrun -rerun to evaluate the energies with and
without dihedral restraints turned on? It seems like you're telling me I'll get garbage out, but this
doesn't seem OK: Even if there are periodicity problems, you should be able to process the same
coordinates with or without restraints and get the same answers, yes?

That's what's going on here, after all: We're evaluating the angle in the same gro file and we get two
different answers depending on the topology file. If I wrote my own code to compute the dihedral
angle, this would not be the case, and it doesn't seem like a desirable situation.

#7 Updated by David van der Spoel about 13 years ago

You have defined the angle between vectors from (atom numbers from 0)
1954 - 982
2236 - 982
x[ 982]={ 6.61300e+00, -2.34900e+00, 3.42100e+00}
x[ 1954]={ 5.50100e+00, -2.41400e+00, 2.95000e+00}
x[ 2236]={ 6.29100e+00, 4.37000e+00, 2.11100e+00}
we hence have
dx1 = 1.102, 0.065, 0.471
dx2 = 0.322,-6.719, 1.310

Now the problem is in the y coordinate of atom 2236 (part of your ligand). If
you define the angle restraint, gromacs will treat protein + ligand as one
molecule and move the ligand in the right place (i.e. it will subtract a box
length of 8.3). If you don't define the angle restraint these molecules will be
treated separately.

To avoid this problem you will have to provide starting coordinates in which the
complex is one piece even without periodic boundary condition. If you analyze
a trajectory in which the molecules were treated separately they may dissociate
in PBC (but not in reality) and you can have the same problem.

You are somehow implying that your topologies should be treated as identical,
even though they are not, as explained above.

I think that your results are correct if you add the angle restraint to the
topology and not otherwise. Please try computing some angles in an independent
manner to convince yourself.

#8 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.

#9 Updated by David van der Spoel about 13 years ago

I absolutely agree it is non-intuitive. It's not clear to me how to improve it
though. Of course one can take the four atoms and compute the angle as
indicated. However, you also expect the programs to respect the periodic
boundary conditions, that is to put the molecules in the box.

In formulating these thoughts, it occurs to me that maybe there is a bug
(missing feature) anyway, in the case that you do not specify the extra
interaction between protein and ligand, you may still expect the programs to
compute the dihedral correctly. Basically, intermolecular interactions are
treated different from intramolecular ones, which should not be the case. (here
molecule is defined in the chemical sense, not the topology sense).

Allow me to come back to you later.

#10 Updated by Berk Hess about 13 years ago

This problem with pbc for two non-interaction molecules CAN NOT be solved.
If the molecules do not interact via an interaction that keeps their
distance within some bounds, there is no way to tell how they should
be "made whole".
If their distance is free to vary unlimitedly, like in your case,
after a long time the two molecules could have a net relative
displacement of say 100 box lengths. They would still be close in space,
but "made whole" the way you want they would be separated by 100 box lengths.

The point is that can not define two unlinked molecules as one entity
over longer timescales than it takes to diffuse apart half a box length.

#11 Updated by David Mobley about 13 years ago

Berk, David:

Berk wrote:

This problem with pbc for two non-interaction molecules CAN NOT be solved.
If the molecules do not interact via an interaction that keeps their
distance within some bounds, there is no way to tell how they should
be "made whole".

OK, but in this case, they are bound to one another.

I am not asking you to come up with a unique way of calculating the dihedral angle between two
molecules that will always give you the same answer regardless of PBC. I'm just asking you to come up
with a consistent way of calculating the dihedral angle between two molecules. That is, if I give you a
trajectory and a topology file and ask you to calculate the dihedral angles between a set of atoms, you
should give me back the same dihedral angles for that trajectory regardless of whether the topology
contains restraints or not.

I well agree that if I give you a trajectory and ask you to do that, and then I assign the same trajectory
to someone David, say, and ask him to do it, too, you two could both come up with different
conventions for whether to wrap the atoms back into the box or not and, depending on how you dealt
with this, you might disagree with one another on what the "right" dihedral angles are. I expect THAT.
I'm just saying that whatever convention you come up with should be applied consistently.

I mean, don't you agree that it makes no sense for GROMACS to, say, return two entirely different
dihedral angles for the same atoms in the same structure, depending on whether or not there is a
restraint (even with spring constant 0!) listed in the topology?

If their distance is free to vary unlimitedly, like in your case,
after a long time the two molecules could have a net relative
displacement of say 100 box lengths. They would still be close in space,
but "made whole" the way you want they would be separated by 100 box lengths.

The ligand is bound to the protein, so it is NOT varying in an unlimited fashion, or even drifting one
box length. It stays put.

The point is that can not define two unlinked molecules as one entity
over longer timescales than it takes to diffuse apart half a box length.

Again, they are bound together. No diffusion.

David wrote:

In formulating these thoughts, it occurs to me that maybe there is a bug
(missing feature) anyway, in the case that you do not specify the extra
interaction between protein and ligand, you may still expect the programs to
compute the dihedral correctly. Basically, intermolecular interactions are
treated different from intramolecular ones, which should not be the case. (here
molecule is defined in the chemical sense, not the topology sense).

I totally agree with this. My big problem here is just consistency: Currently, the dihedral angle is not
being handled in a consistent way.
I agree that, as Berk points out, in general there is no way to define a unique dihedral if, say, two
molecules are drifting away from one another, or even anytime there is a periodic box. I'm not asking
for uniqueness, just consistency in the intermolecular/intramolecular (topology) sense: If I store a
trajectory or set of coordinate files, I should get the same answer for a dihedral angle from that
trajectory regardless of whether I have a restraint line in my topology or not.

#12 Updated by David Mobley about 13 years ago

To reiterate my last comment, let me quote another of David's earlier remarks:

Now the problem is in the y coordinate of atom 2236 (part of your ligand). If
you define the angle restraint, gromacs will treat protein + ligand as one
molecule and move the ligand in the right place (i.e. it will subtract a box
length of 8.3). If you don't define the angle restraint these molecules will be
treated separately.

This, in my mind, is the consistency problem I keep talking about.

I am not sure how to solve it. Maybe just remove the "feature" that things with restraints between them
be mapped into the same box, or add an option to turn it off? (But this might break other things -- I am
not sure). Or, make sure that things that are part of the same "moleule" (in the topology sense) are
ALWAYS mapped into the same box, regardless of whether or not there are restraints between them
(which would mean that adding restraints wouldn't change how this is handled).

David

#13 Updated by David van der Spoel about 13 years ago

David,

can you generate two new tpr files, identical to the previous but with pbc =
full and rerun the analysis? It wouldn't hurt to upload them as well.

#14 Updated by David Mobley about 13 years ago

Created an attachment (id=66)
tpr file with pbc full with and without restraints in topology

#15 Updated by David Mobley about 13 years ago

I still get the same thing. I uploaded the tpr files.

#16 Updated by David van der Spoel about 13 years ago

Fixed the bug by forcing periodicity computation in the analysis of dihedrals
and angles. CVS updated for file src/tools/anadih.c. This fix effects g_angle
and g_chi.

#17 Updated by David Mobley about 13 years ago

OK, David, I just tried out the version of g_angle I got from the release 3.3
patches branch of cvs, and I get 90^o for every dihedral angle I measure. Can
you check it?

Thanks.

#18 Updated by David Mobley about 13 years ago

Guys,

Can you PLEASE take a look at this in the near future? I have a student who is
stalled until we can get this resolved... The problem is twofold:
(a) In the current release 3.3 patches branch, g_angle always gives a dihedral
of 90^o regardless of the actual dihedral
(b) Before that "fix" that caused (a), g_angle gives a dihedral in some cases
that is inconsistent with the value used internally in GROMACS, meaning that if
I use g_angle to measure a dihedral, and then try to use GROMACS to restrain
the dihedral to stay near that same value (in some cases, as in the tprs I
supplied before), I get huge restraining energies. For example, g_angle might
tell me the dihedral angle between four atoms is 66^o; I turn on restraints to
keep it there, but GROMACS internally thinks the dihedral angle starts off as
-70 (or some such) and so the restraining energy to keep it at 66^o is very
large. What I really WANT is to keep it as -70.

I think the fix that is needed is just to make sure that the analysis tool
g_angle computes a dihedral angle in the SAME WAY that it is computed internally
in GROMACS, so that if I compute a dihedral angle with g_angle and then restrain
to that angle, the two will be consistent.

#19 Updated by David van der Spoel about 13 years ago

This gives me a deja-vu feeling... IIRC this was caused by a missing include
file. Is this happening with the same input still? Have you recompiled gromacs
completely or just g_angle?

#20 Updated by David van der Spoel about 13 years ago

By the way, which hardware does this run on?

Did you do a complete cvs update of the source tree?

#21 Updated by David Mobley about 13 years ago

David,

I tried two things. First, a couple weeks ago, I compiled from scratch. Then,
when I tested it, I got this problem. Then yesterday, I did "cvs update" and
compiled again (without doing make clean) and still had the same problem.

This is on a Xeon with Fedora Core 4.

Thanks,
David

#22 Updated by David van der Spoel about 13 years ago

make distclean might help.
can you please confirm that the attached files still produce the problem, or
otherwise provide new ones?

#23 Updated by David Mobley about 13 years ago

David,

Yes, it is still broken, with those input files (I just tested again) or with
any other input files. Regardless of what input file I use, I always get a
dihedral angle of 90^o for ANY dihedral angle.

I tested this on two architectures, now: Xeon and Opteron. Same thing. I
compiled these independently from scratch from the cvs about two weeks ago.

I will try again by doing a cvs update and then compiling after make distclean,
but I assume it won't make any differences.

I AM still following the instructions on the web page for getting the 3.3
patches branch from CVS. I recall you had commented beofre that those might be
incorrect, so if there is some other procedure I should be following to get
these from CVS, please let me know.

#24 Updated by David van der Spoel about 13 years ago

Apparantly there was a left over from the old code, which I didn't take away.
It's fixed in CVS now.

#25 Updated by David Mobley about 13 years ago

OK, thanks. We think this is fixed now, but will keep you posted.

Also available in: Atom PDF