Fix and cover RMSD

package/MDAnalysis/analysis/rms.py

@@ -338,9 +338,9 @@ def __init__(self, traj, reference=None, select='all',
             self.reference = self.universe
         else:
             self.reference = reference
-        self.select = _process_selection(select)
+        self.select = process_selection(select)
         if groupselections is not None:
-            self.groupselections = [_process_selection(s) for s in groupselections]
+            self.groupselections = [process_selection(s) for s in groupselections]
         else:
             self.groupselections = []
         self.mass_weighted = mass_weighted
@@ -351,12 +351,12 @@ def __init__(self, traj, reference=None, select='all',
         self.ref_atoms = self.reference.select_atoms(*self.select['reference'])
         self.traj_atoms = self.universe.select_atoms(*self.select['mobile'])
         if len(self.ref_atoms) != len(self.traj_atoms):
-            logger.exception()
-            raise SelectionError("Reference and trajectory atom selections do "
-                                 "not contain the same number of atoms: "
-                                 "N_ref={0:d}, N_traj={1:d}".format(
-                                     self.ref_atoms.n_atoms,
-                                     self.traj_atoms.n_atoms))
+            errmsg = ("Reference and trajectory atom selections do "+
+                      "not contain the same number of atoms: "+
+                      "N_ref={0:d}, N_traj={1:d}".format(self.ref_atoms.n_atoms,
+                                                        self.traj_atoms.n_atoms))
+            logger.exception(errmsg)
+            raise SelectionError(errmsg)
         logger.info("RMS calculation for {0:d} atoms.".format(len(self.ref_atoms)))
         mass_mismatches = (np.absolute(self.ref_atoms.masses - self.traj_atoms.masses) > self.tol_mass)
         if np.any(mass_mismatches):
@@ -366,8 +366,9 @@ def __init__(self, traj, reference=None, select='all',
                 if ar.name != at.name:
                     logger.error("{0!s:>4} {1:3d} {2!s:>3} {3!s:>3} {4:6.3f}  |  {5!s:>4} {6:3d} {7!s:>3} {8!s:>3} {9:6.3f}".format(ar.segid, ar.resid, ar.resname, ar.name, ar.mass,
                                  at.segid, at.resid, at.resname, at.name, at.mass))
-            errmsg = "Inconsistent selections, masses differ by more than {0:f}; mis-matching atoms are shown above.".format( \
-                     self.tol_mass)
+            errmsg = ("Inconsistent selections, masses differ by more than"
+                     + "{0:f}; mis-matching atoms are shown above.".format(
+                     self.tol_mass))
             logger.error(errmsg)
             raise SelectionError(errmsg)
         del mass_mismatches
@@ -386,7 +387,7 @@ def __init__(self, traj, reference=None, select='all',
         for igroup, (sel, atoms) in enumerate(zip(self.groupselections,
                                                   self.groupselections_atoms)):
             if len(atoms['mobile']) != len(atoms['reference']):
-                logger.exception()
+                logger.exception('SelectionError: Group Selection')
                 raise SelectionError(
                     "Group selection {0}: {1} | {2}: Reference and trajectory "
                     "atom selections do not contain the same number of atoms: "

testsuite/MDAnalysisTests/analysis/test_rms.py

@@ -21,16 +21,19 @@
 import MDAnalysis as mda
 from MDAnalysis.analysis import rms, align
 
-from numpy.testing import TestCase, assert_almost_equal, raises, assert_
+from numpy.testing import (TestCase, assert_almost_equal, raises, assert_,
+                           assert_array_almost_equal)
+
 import numpy as np
 
 import os
 
+from MDAnalysis.exceptions import SelectionError, NoDataError
 from MDAnalysisTests.datafiles import GRO, XTC, rmsfArray, PSF, DCD
 from MDAnalysisTests import tempdir
 
 
-class TestRMSD(object):
+class Testrmsd(object):
     def __init__(self):
         shape = (5, 3)
         # vectors with length one
@@ -128,6 +131,82 @@ def test_with_superposition_equal(self):
         rmsd_superposition = rms.rmsd(A, B, center=True, superposition=True)
         assert_almost_equal(rmsd, rmsd_superposition)
 
+class TestRMSD(object):
+    def setUp(self):
+        self.universe = MDAnalysis.Universe(PSF, DCD)
+        self.tempdir = tempdir.TempDir()
+        self.outfile = os.path.join(self.tempdir.name, 'rmsd.txt')
+        self.correct_values = [[0, 0, 0], [49, 48.9999, 4.68953]]
+        self.correct_values_group = [[0 , 0, 0, 0, 0],
+                                     [49, 48.9999, 4.7857, 4.7002,
+                                      4.68981]]
+
+    def tearDown(self):
+        del self.universe
+        del self.tempdir
+
+    def test_rmsd(self):
+        RMSD = MDAnalysis.analysis.rms.RMSD(self.universe, select='name CA')
+        RMSD.run(step=49)
+        assert_array_almost_equal(RMSD.rmsd, self.correct_values, 4,
+                                  err_msg="error: rmsd profile should match" +
+                                  "test values")
+
+    def test_rmsd_single_frame(self):
+        RMSD = MDAnalysis.analysis.rms.RMSD(self.universe, select='name CA')
+        RMSD.run(start=5, stop=6)
+        single_frame = [[ 5, 5, 0.91544906]]
+        assert_array_almost_equal(RMSD.rmsd, single_frame, 4,
+                                  err_msg="error: rmsd profile should match" +
+                                  "test values")
+
+    def test_mass_weighted_and_save(self):
+        RMSD = MDAnalysis.analysis.rms.RMSD(self.universe, select='name CA')
+        RMSD.run(step=49, mass_weighted= True)
+        RMSD.save(self.outfile)
+        saved = np.loadtxt(self.outfile)
+        assert_array_almost_equal(RMSD.rmsd, saved, 4,
+                                  err_msg="error: rmsd profile should match " +
+                                          "test values")
+
+    def test_rmsd_group_selections(self):
+        RMSD = MDAnalysis.analysis.rms.RMSD(self.universe,
+                                            groupselections=
+                                            ['backbone','name CA'])
+        RMSD.run(step=49)
+
+        assert_array_almost_equal(RMSD.rmsd, self.correct_values_group, 4,
+                                  err_msg="error: rmsd profile should match" +
+                                  "test values")
+
+
+    @raises(SelectionError)
+    def test_ref_length_unequal_len(self):
+        reference = MDAnalysis.Universe(PSF, DCD)
+        reference.atoms = reference.atoms[:-1]
+        RMSD = MDAnalysis.analysis.rms.RMSD(self.universe,
+                                            reference=reference)
+
+    @raises(SelectionError)
+    def test_mass_mismatches(self):
+        reference = MDAnalysis.Universe(PSF, DCD)
+        reference.atoms.masses = 10
+        RMSD = MDAnalysis.analysis.rms.RMSD(self.universe,
+                                            reference=reference)
+
+
+    @raises(SelectionError)
+    def test_group_selections_unequal_len(self):
+        reference = MDAnalysis.Universe(PSF, DCD)
+        reference.atoms[0].resname='NOTMET'
+        RMSD = MDAnalysis.analysis.rms.RMSD(self.universe,
+                                            reference=reference,
+                                            groupselections=
+                                            ['resname MET','type NH3'])
+    @raises(NoDataError)
+    def test_save_before_run(self):
+        RMSD = MDAnalysis.analysis.rms.RMSD(self.universe)
+        RMSD.save('blah')
 
 class TestRMSF(TestCase):
     def setUp(self):