Skip to content

Commit

Permalink
Switch rms.RMSD and rms.RMSF to using the Results class (#3277)
Browse files Browse the repository at this point in the history 
* Fixes #3274
* Towards #3261
* Switch rms.RMSD and rms.RMSF to using the Results class
  • Loading branch information
@IAlibay
IAlibay authored May 7, 2021
1 parent a43c0a0 commit da46e84
Show file tree
Hide file tree
Showing 3 changed files with 96 additions and 33 deletions.
6 changes: 6 additions & 0 deletions package/CHANGELOG
View file
Original file line number Diff line number Diff line change
Expand Up @@ -174,6 +174,8 @@ Enhancements
checking if it can be used in parallel analysis. (Issue #2996, PR #2950)

Changes
* `analysis.rms.RMSD` and `analysis.rms.RMSF` classes now store `rmsd` and
`rmsf` data using the `analysis.base.Results` class (Issues #3274 #3261)
* `analysis.dihedrals` classes now store angle data using the
`results.angles` attribute (Issue #3261)
* `msd.EinsteinMSD` now uses the `analysis.base.Results` class to store
Expand Down Expand Up @@ -233,6 +235,10 @@ Changes
* Added OpenMM coordinate and topology converters (Issue #2863, PR #2917)

Deprecations
* The `analysis.rms.RMSD.rmsd` and `analysis.rms.RMSF.rmsf` attributes are
now deprecated in favour of `analysis.rms.RMSD.results.rmsd` and
`analysis.rms.RMSF.results.rmsf` respectively. They will be removed in
3.0.0 (Issues #3274, #3261)
* The `angles` attribute for the `dihedrals` classes (Dihedral,
Ramachandran, Janin) is now deprecated in favour of `results.angles`. It
will be removed in 3.0.0 (Issue #3261)
Expand Down
76 changes: 59 additions & 17 deletions package/MDAnalysis/analysis/rms.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,7 +73,7 @@
* calculate the backbone RMSD and RMSD for CORE, LID, NMP (backbone atoms)
The trajectory is included with the test data files. The data in
:attr:`RMSD.rmsd` is plotted with :func:`matplotlib.pyplot.plot`::
:attr:`RMSD.results.rmsd` is plotted with :func:`matplotlib.pyplot.plot`::
import MDAnalysis
from MDAnalysis.tests.datafiles import PSF,DCD,CRD
Expand All @@ -91,7 +91,7 @@
R.run()
import matplotlib.pyplot as plt
rmsd = R.rmsd.T # transpose makes it easier for plotting
rmsd = R.results.rmsd.T # transpose makes it easier for plotting
time = rmsd[1]
fig = plt.figure(figsize=(4,4))
ax = fig.add_subplot(111)
Expand All @@ -117,20 +117,41 @@
:members:
:inherited-members:
.. attribute:: rmsd
.. attribute:: results.rmsd
Contains the time series of the RMSD as an N×3 :class:`numpy.ndarray`
array with content ``[[frame, time (ps), RMSD (A)], [...], ...]``.
.. versionadded:: 2.0.0
.. attribute:: rmsd
Alias to the :attr:`results.rmsd` attribute.
.. deprecated:: 2.0.0
Will be removed in MDAnalysis 3.0.0. Please use :attr:`results.rmsd`
instead.
.. autoclass:: RMSF
:members:
:inherited-members:
.. attribute:: rmsf
.. attribute:: results.rmsf
Results are stored in this N-length :class:`numpy.ndarray` array,
giving RMSFs for each of the given atoms.
.. versionadded:: 2.0.0
.. attribute:: rmsf
Alias to the :attr:`results.rmsf` attribute.
.. deprecated:: 2.0.0
Will be removed in MDAnalysis 3.0.0. Please use :attr:`results.rmsf`
instead.
"""
import numpy as np

Expand Down Expand Up @@ -316,11 +337,15 @@ class RMSD(AnalysisBase):
Run the analysis with :meth:`RMSD.run`, which stores the results
in the array :attr:`RMSD.rmsd`.
in the array :attr:`RMSD.results.rmsd`.
.. versionchanged:: 1.0.0
``save()`` method was removed, use ``np.savetxt()`` on
:attr:`RMSD.rmsd` instead.
:attr:`RMSD.results.rmsd` instead.
.. versionchanged:: 2.0.0
:attr:`rmsd` results are now stored in a
:class:`MDAnalysis.analysis.base.Results` instance.
"""
def __init__(self, atomgroup, reference=None, select='all',
Expand Down Expand Up @@ -628,8 +653,8 @@ def _prepare(self):
else:
self._rot = None

self.rmsd = np.zeros((self.n_frames,
3 + len(self._groupselections_atoms)))
self.results.rmsd = np.zeros((self.n_frames,
3 + len(self._groupselections_atoms)))

self._mobile_coordinates64 = self.mobile_atoms.positions.copy().astype(np.float64)

Expand All @@ -638,7 +663,7 @@ def _single_frame(self):
self._mobile_coordinates64[:] = self.mobile_atoms.positions
self._mobile_coordinates64 -= mobile_com

self.rmsd[self._frame_index, :2] = self._ts.frame, self._trajectory.time
self.results.rmsd[self._frame_index, :2] = self._ts.frame, self._trajectory.time

if self._groupselections_atoms:
# superimpose structures: MDAnalysis qcprot needs Nx3 coordinate
Expand All @@ -647,7 +672,7 @@ def _single_frame(self):
# left** so that we can easily use broadcasting and save one
# expensive numpy transposition.

self.rmsd[self._frame_index, 2] = qcp.CalcRMSDRotationalMatrix(
self.results.rmsd[self._frame_index, 2] = qcp.CalcRMSDRotationalMatrix(
self._ref_coordinates64, self._mobile_coordinates64,
self._n_atoms, self._rot, self.weights_select)

Expand All @@ -666,17 +691,25 @@ def _single_frame(self):
for igroup, (refpos, atoms) in enumerate(
zip(self._groupselections_ref_coords64,
self._groupselections_atoms), 3):
self.rmsd[self._frame_index, igroup] = rmsd(
self.results.rmsd[self._frame_index, igroup] = rmsd(
refpos, atoms['mobile'].positions,
weights=self.weights_groupselections[igroup-3],
center=False, superposition=False)
else:
# only calculate RMSD by setting the Rmatrix to None (no need
# to carry out the rotation as we already get the optimum RMSD)
self.rmsd[self._frame_index, 2] = qcp.CalcRMSDRotationalMatrix(
self.results.rmsd[self._frame_index, 2] = qcp.CalcRMSDRotationalMatrix(
self._ref_coordinates64, self._mobile_coordinates64,
self._n_atoms, None, self.weights_select)

@property
def rmsd(self):
wmsg = ("The `rmsd` attribute was deprecated in MDAnalysis 2.0.0 and "
"will be removed in MDAnalysis 3.0.0. Please use "
"`results.rmsd` instead.")
warnings.warn(wmsg, DeprecationWarning)
return self.results.rmsd


class RMSF(AnalysisBase):
r"""Calculate RMSF of given atoms across a trajectory.
Expand All @@ -691,7 +724,7 @@ class RMSF(AnalysisBase):
Run the analysis with :meth:`RMSF.run`, which stores the results
in the array :attr:`RMSF.rmsf`.
in the array :attr:`RMSF.results.rmsf`.
"""
def __init__(self, atomgroup, **kwargs):
Expand Down Expand Up @@ -781,7 +814,8 @@ def __init__(self, atomgroup, **kwargs):
in_memory=True).run()
The trajectory is now fitted to the reference (the RMSD is stored as
`aligner.rmsd` for further inspection). Now we can calculate the RMSF::
`aligner.results.rmsd` for further inspection). Now we can calculate
the RMSF::
from MDAnalysis.analysis.rms import RMSF
Expand All @@ -792,7 +826,7 @@ def __init__(self, atomgroup, **kwargs):
import matplotlib.pyplot as plt
plt.plot(calphas.resnums, rmsfer.rmsf)
plt.plot(calphas.resnums, rmsfer.results.rmsf)
Expand Down Expand Up @@ -829,8 +863,16 @@ def _single_frame(self):

def _conclude(self):
k = self._frame_index
self.rmsf = np.sqrt(self.sumsquares.sum(axis=1) / (k + 1))
self.results.rmsf = np.sqrt(self.sumsquares.sum(axis=1) / (k + 1))

if not (self.rmsf >= 0).all():
if not (self.results.rmsf >= 0).all():
raise ValueError("Some RMSF values negative; overflow " +
"or underflow occurred")

@property
def rmsf(self):
wmsg = ("The `rmsf` attribute was deprecated in MDAnalysis 2.0.0 and "
"will be removed in MDAnalysis 3.0.0. Please use "
"`results.rmsd` instead.")
warnings.warn(wmsg, DeprecationWarning)
return self.results.rmsf
47 changes: 31 additions & 16 deletions testsuite/MDAnalysisTests/analysis/test_rms.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -185,14 +185,14 @@ def correct_values_backbone_group(self):
def test_rmsd(self, universe, correct_values):
RMSD = MDAnalysis.analysis.rms.RMSD(universe, select='name CA')
RMSD.run(step=49)
assert_almost_equal(RMSD.rmsd, correct_values, 4,
assert_almost_equal(RMSD.results.rmsd, correct_values, 4,
err_msg="error: rmsd profile should match" +
"test values")

def test_rmsd_unicode_selection(self, universe, correct_values):
RMSD = MDAnalysis.analysis.rms.RMSD(universe, select=u'name CA')
RMSD.run(step=49)
assert_almost_equal(RMSD.rmsd, correct_values, 4,
assert_almost_equal(RMSD.results.rmsd, correct_values, 4,
err_msg="error: rmsd profile should match" +
"test values")

Expand All @@ -202,13 +202,13 @@ def test_rmsd_atomgroup_selections(self, universe):
select="resid 1-30").run()
R2 = MDAnalysis.analysis.rms.RMSD(universe.atoms.select_atoms("name CA"),
select="resid 1-30").run()
assert not np.allclose(R1.rmsd[:, 2], R2.rmsd[:, 2])
assert not np.allclose(R1.results.rmsd[:, 2], R2.results.rmsd[:, 2])

def test_rmsd_single_frame(self, universe):
RMSD = MDAnalysis.analysis.rms.RMSD(universe, select='name CA',
).run(start=5, stop=6)
single_frame = [[5, 6, 0.91544906]]
assert_almost_equal(RMSD.rmsd, single_frame, 4,
assert_almost_equal(RMSD.results.rmsd, single_frame, 4,
err_msg="error: rmsd profile should match" +
"test values")

Expand All @@ -218,14 +218,14 @@ def test_mass_weighted(self, universe, correct_values):
RMSD = MDAnalysis.analysis.rms.RMSD(universe, select='name CA',
weights='mass').run(step=49)

assert_almost_equal(RMSD.rmsd, correct_values, 4,
assert_almost_equal(RMSD.results.rmsd, correct_values, 4,
err_msg="error: rmsd profile should match"
"test values")

def test_custom_weighted(self, universe, correct_values_mass):
RMSD = MDAnalysis.analysis.rms.RMSD(universe, weights="mass").run(step=49)

assert_almost_equal(RMSD.rmsd, correct_values_mass, 4,
assert_almost_equal(RMSD.results.rmsd, correct_values_mass, 4,
err_msg="error: rmsd profile should match"
"test values")

Expand All @@ -234,15 +234,15 @@ def test_weights_mass_is_mass_weighted(self, universe):
weights="mass").run(step=49)
RMSD_cust = MDAnalysis.analysis.rms.RMSD(universe,
weights=universe.atoms.masses).run(step=49)
assert_almost_equal(RMSD_mass.rmsd, RMSD_cust.rmsd, 4,
assert_almost_equal(RMSD_mass.results.rmsd, RMSD_cust.results.rmsd, 4,
err_msg="error: rmsd profiles should match for 'mass' "
"and universe.atoms.masses")

def test_custom_weighted_list(self, universe, correct_values_mass):
weights = universe.atoms.masses
RMSD = MDAnalysis.analysis.rms.RMSD(universe,
weights=list(weights)).run(step=49)
assert_almost_equal(RMSD.rmsd, correct_values_mass, 4,
assert_almost_equal(RMSD.results.rmsd, correct_values_mass, 4,
err_msg="error: rmsd profile should match" +
"test values")

Expand All @@ -253,7 +253,7 @@ def test_custom_groupselection_weights_applied_1D_array(self, universe):
weights=None,
weights_groupselections=[[1, 0, 0, 0, 0], None]).run(step=49)

assert_almost_equal(RMSD.rmsd.T[3], RMSD.rmsd.T[4], 4,
assert_almost_equal(RMSD.results.rmsd.T[3], RMSD.results.rmsd.T[4], 4,
err_msg="error: rmsd profile should match "
"for applied weight array and selected resid")

Expand All @@ -265,7 +265,7 @@ def test_custom_groupselection_weights_applied_mass(self, universe, correct_valu
weights_groupselections=['mass',
universe.atoms.masses]).run(step=49)

assert_almost_equal(RMSD.rmsd.T[3], RMSD.rmsd.T[4], 4,
assert_almost_equal(RMSD.results.rmsd.T[3], RMSD.results.rmsd.T[4], 4,
err_msg="error: rmsd profile should match "
"between applied mass and universe.atoms.masses")

Expand Down Expand Up @@ -308,7 +308,7 @@ def test_rmsd_group_selections(self, universe, correct_values_group):
RMSD = MDAnalysis.analysis.rms.RMSD(universe,
groupselections=['backbone', 'name CA']
).run(step=49)
assert_almost_equal(RMSD.rmsd, correct_values_group, 4,
assert_almost_equal(RMSD.results.rmsd, correct_values_group, 4,
err_msg="error: rmsd profile should match"
"test values")

Expand All @@ -321,7 +321,7 @@ def test_rmsd_backbone_and_group_selection(self, universe,
groupselections=['backbone and resid 1:10',
'backbone and resid 10:20']).run(step=49)
assert_almost_equal(
RMSD.rmsd, correct_values_backbone_group, 4,
RMSD.results.rmsd, correct_values_backbone_group, 4,
err_msg="error: rmsd profile should match test values")

def test_ref_length_unequal_len(self, universe):
Expand All @@ -347,7 +347,7 @@ def test_ref_mobile_mass_mismapped(self, universe,correct_values_mass_add_ten):
weights='mass',
tol_mass=100)
RMSD.run(step=49)
assert_almost_equal(RMSD.rmsd, correct_values_mass_add_ten, 4,
assert_almost_equal(RMSD.results.rmsd, correct_values_mass_add_ten, 4,
err_msg="error: rmsd profile should match "
"between true values and calculated values")

Expand All @@ -359,6 +359,14 @@ def test_group_selections_unequal_len(self, universe):
reference=reference,
groupselections=['resname MET', 'type NH3'])

def test_rmsd_attr_warning(self, universe):
RMSD = MDAnalysis.analysis.rms.RMSD(
universe, select='name CA').run(stop=2)

wmsg = "The `rmsd` attribute was deprecated in MDAnalysis 2.0.0"
with pytest.warns(DeprecationWarning, match=wmsg):
assert_equal(RMSD.rmsd, RMSD.results.rmsd)


class TestRMSF(object):
@pytest.fixture()
Expand All @@ -370,14 +378,14 @@ def test_rmsf(self, universe):
rmsfs.run()
test_rmsfs = np.load(rmsfArray)

assert_almost_equal(rmsfs.rmsf, test_rmsfs, 5,
assert_almost_equal(rmsfs.results.rmsf, test_rmsfs, 5,
err_msg="error: rmsf profile should match test "
"values")

def test_rmsf_single_frame(self, universe):
rmsfs = rms.RMSF(universe.select_atoms('name CA')).run(start=5, stop=6)

assert_almost_equal(rmsfs.rmsf, 0, 5,
assert_almost_equal(rmsfs.results.rmsf, 0, 5,
err_msg="error: rmsfs should all be zero")

def test_rmsf_identical_frames(self, universe, tmpdir):
Expand All @@ -392,5 +400,12 @@ def test_rmsf_identical_frames(self, universe, tmpdir):
universe = mda.Universe(GRO, outfile)
rmsfs = rms.RMSF(universe.select_atoms('name CA'))
rmsfs.run()
assert_almost_equal(rmsfs.rmsf, 0, 5,
assert_almost_equal(rmsfs.results.rmsf, 0, 5,
err_msg="error: rmsfs should all be 0")

def test_rmsf_attr_warning(self, universe):
rmsfs = rms.RMSF(universe.select_atoms('name CA')).run(stop=2)

wmsg = "The `rmsf` attribute was deprecated in MDAnalysis 2.0.0"
with pytest.warns(DeprecationWarning, match=wmsg):
assert_equal(rmsfs.rmsf, rmsfs.results.rmsf)

0 comments on commit da46e84

Please sign in to comment.