[gmx-users] g_rms on CA atoms only, with and without mass weighting -nomw

[Jan Domanski]

Hi,

I'm using gromacs 4.5.4 and I've got a detailed question on how the mass
weighting works.

Given a trajectory and a pdb from
http://code.google.com/p/mdanalysis/source/browse/testsuite/MDAnalysisTests/data/adk_oplsaa.pdb
http://code.google.com/p/mdanalysis/source/browse/testsuite/MDAnalysisTests/data/adk_oplsaa.trr

I called the following
structure="adk_oplsaa.pdb"
trajectory="adk_oplsaa.trr"
g_rms -s $structure -f $trajectory -o d_rmsd.xvg << EOF
3
3
EOF

g_rms -nomw -s $structure -f $trajectory -o d_rmsd_nomw.xvg << EOF
3
3
EOF

I was expecting the RMSD values to be identical: all the CA
atoms have identical weights with the -mw flag on, which (to my mind)
should be yield same results to when -nomw is specified and the same set
of atoms is selected. The question: is this an unreasonable expectation?

sdiff d_rmsd.xvg  d_rmsd_nomw.xvg  | tail # consistent only for 5 sig figs

  0.0000000    0.0000369 |       0.0000000    0.0000368
100.0000076    0.9818456 |     100.0000076    0.9818469
200.0000153    0.8223789 |     200.0000153    0.8223798
300.0000000    0.6450028 |     300.0000000    0.6450034
400.0000305    0.8285408 |     400.0000305    0.8285414
500.0000305    1.9386379 |     500.0000305    1.9386411
600.0000000    1.6888058 |     600.0000000    1.6888076
700.0000610    1.6885630 |     700.0000610    1.6885644
800.0000610    1.8211102 |     800.0000610    1.8211110
900.0000610    2.1306074 |     900.0000610    2.1306095

(BTW, the g_rms -h mentions something about a '-debug flag' but it
seems not to be working.)

As a check, I've used MDAnalysis 0.7.6-devel
(http://mdanalysis.googlecode.com) and the following code to get RMSD on
the same data, which was consistent over 14 sig figs.

from MDAnalysis import *
from MDAnalysis.analysis.align import *
conf = "adk_oplsaa.pdb"
traj = "adk_oplsaa.trr"
ref = Universe(conf)
mob = Universe(conf, traj)
rms_fit_trj(mob, ref, select="name CA", rmsdfile="rmsd.dat",
mass_weighted=True)
rms_fit_trj(mob, ref, select="name CA", rmsdfile="rmsd_nomw.dat",
mass_weighted=False)

sdiff rmsd.dat  rmsd_nomw.dat # consistent over 14 sig figs
3.686669177327545608e-04 |	3.686671021799712840e-04
9.818467126090068220e+00 |	9.818467126090039798e+00
8.223796870415890581e+00	8.223796870415890581e+00
6.450033517100112412e+00 |	6.450033517100091096e+00
8.285414822442923821e+00 |	8.285414822442906058e+00
1.938640271585943964e+01 |	1.938640271585941832e+01
1.688807649241958941e+01 |	1.688807649241956454e+01
1.688564785098214571e+01 |	1.688564785098212795e+01
1.821111456215330549e+01 |	1.821111456215327706e+01
2.130609827021595137e+01 |	2.130609827021593006e+01

Thanks for helping me figure it out guys,


Jan