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

[Mark Abraham]

On 5/07/2012 10:11 AM, Jan Domanski wrote:
> 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.)

See manual D.1 - the -debug flag takes an argument.

>
> 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,

This is expected. See 
http://www.gromacs.org/Documentation/Floating_Point_Arithmetic. Mass 
weighting changes the manner in which round-off error accumulates. 
Double precision is less affected by this. Your g_rms is single 
precision, and apparently MDAnalysis is double.

Mark