Re: Trouble using rotateReference, centerRefrence and refPositionsGroup in the colvar module

From: Giacomo Fiorin (
Date: Thu Feb 02 2012 - 17:02:08 CST

Hi Ajasja, here are the comments.

The gradients of the dihedral with respect of the atoms of group1 are
always projected correctly from the rotated frame to the laboratory, while
the gradients of the dihedral with respect to the backbone atoms (which
enter the dihedral's formula through the rotation matrix) are neglected.
This is a very good approximation for as long as the internal structure of
the group you're using to do the fitting (i.e. backbone atoms) doesn't

Of course, this is not true for rmsd and eigenvector, where you're supposed
to change the same structure you're fitting, which is why the exception.

What the second warning says is that if you're applying a torque on the
group in the rotated frame, you won't see the effects of that torque in the
rotated frame, because you always fit the group to the reference
coordinates. But the forces from the torque are still applied, in the
laboratory frame. If the torque comes from a restraint of some sort, and
that restraint is defined in the rotated frame, you may end up having a
constant force.

In your case you're 100% fine on this because you're only applying a force
to the center of geometry of group1, so no torque.


On Thu, Feb 2, 2012 at 10:06 AM, Ajasja Ljubetič

> Hum, this is a very peculiar situation... I wish you had defined the
>> variable in terms of atoms that move together with the peptide (such as
>> certain side chains) instead of dummy atoms if you were planning on
>> allowing the protein to move later.
> I was thinking about this too, but I didn't want to be limited only to
> alpha-helices, for example I would like to preform the same calculation on
> a beta-sheet structure. So I was not sure which side chains atoms would
> always be in the correct location.
> Is there any other (theoretical or practical) consideration why defining
> the colvar in this way is a bad idea? The one oneSiteSystemForce options
> seems to work nicely, and the simulation seems stable.
> In any case, I don't believe you need to rotate the dummy atoms. The
>> dummy atoms won't drift or rotate with the protein. Instead, it's the
>> protein atoms that you should fit using the initial structure to go back to
>> the frame where the protein didn't move. So you only apply the fitting
>> options to group1.
> Yes, you are absolutely right! I completely forgot that I can (actually
> have to) transform only one group of the colvar component.
>> Regarding why the PDB file doesn't work, you need to match the number of
>> coordinates of the atoms that you use for fitting (refPositionsGroup), not
>> the number of atoms of group1. The problem in providing a file with the
>> whole system is that the PDB reader will by default read the coordinates of
>> group1, not those of refPositionsGroup.
> Aha, I thought that it would take the indices of the refPositionsGroup and
> only load those from the PDB.
>> If you write a PDB file that just contains the five positions of the
>> refPositionsGroup atoms, it should work.
> Thank you, this works! For completeness I'm pasting the working version
> below. (It was quicker for me to just paste the positions into
> refPositions).
> dihedral {
> oneSiteSystemForce
> group1 { atomnumbers 173 174 175 176 177 181 185 189 193 194
> rotateReference on
> centerReference on
> refPositionsGroup {atomnumbers 200 202 204 208 209}
> refpositions
> (0.9959999918937683,-1.2680000066757202,-0.7979999780654907)
> (2.0460000038146973,-0.8029999732971191,0.07900000363588333)
> (3.1630001068115234,-0.061000000685453415,-0.6990000009536743)
> (1.49399995803833,0.12099999934434891,1.1390000581741333)
> (1.8849999904632568,0.041999999433755875,2.3010001182556152) } ;
> group2 { dummyAtom
> (3.1630001068115234,-0.061000000685453415,-0.6990000009536743) } ;
> group3 { dummyAtom
> (3.1630001068115234,-0.061000000685453415,9.300999999046326) } ;
> group4 { dummyAtom
> (13.163000106811523,-0.061000000685453415,9.300999999046326) } ;
> }
> }
> Update:
> I analysed the first 50 ns of the simulation with the new coordinate
> system. The result does not match the previous one (the PMF seems
> distorted). I'm almost 100% I made some mistake in the configuration file,
> but just to be sure, could you comment on the following two warnings?
> - From the manual under rotateReference/centerRefrence:
> *"Note*:*the derivatives of the colvars with respect to the
> rotation/translation are usually neglected (except by *rmsd* and *
> eigenvector*)*."
> What exactly does this mean? Does this come into play when
> the Cartesian biasing force is calculated (eq 4<>)
> from the colvar derivative? (BTW there is still a small typo in that
> equation - the gradient should be along Cartesian coordinates and not
> along colvars space).
> - From the log file:
> "Warning: atom group "group1" is set to be rotated to a reference
> orientation: a torque different than zero on this group could make the
> simulation unstable. If this happens, set "disableForces" to yes for this
> group."
> If I understand rotateReference/centerRefrence correctly, the
> atoms in the group are first rotated and translated to
> best superimpose with the reference positions, then the value of the colvar
> is calculated (and the derivative at that point) and the Cartesian forces
> are estimated from the average forces on the colvars in this bin. Then
> this Cartesian forces are rotated back into the original reference frame.
> In my case this almost certainly result in a non zero torque, but the
> simulation still seems stable. Should I be worried about this
> non-zero torque?
> Again, thank you very much & best regards,
> Ajasja

This archive was generated by hypermail 2.1.6 : Mon Dec 31 2012 - 23:21:10 CST