diff --git a/waterdynamics/tests/test_waterdynamics.py b/waterdynamics/tests/test_waterdynamics.py
index af68982..eb50ac0 100644
--- a/waterdynamics/tests/test_waterdynamics.py
+++ b/waterdynamics/tests/test_waterdynamics.py
@@ -43,31 +43,48 @@ def universe():
 
 def test_WaterOrientationalRelaxation(universe):
     wor = waterdynamics.WaterOrientationalRelaxation(
-        universe, SELECTION1, 0, 5, 2, order=2)
-    wor.run()
-    assert_almost_equal(wor.timeseries[1][2], 0.35887,
+        universe.atoms, ["name OH2", "name OH2 ", "name H1 or name H2"], 2, 1.2, update_selections=False, order=1)
+    wor.run(start=0, stop=4)
+    assert_almost_equal(wor.results["correlation"][2],  0.7148607552,
                         decimal=5)
 
 
-def test_WaterOrientationalRelaxation_zeroMolecules(universe):
-    wor_zero = waterdynamics.WaterOrientationalRelaxation(
-        universe, SELECTION2, 0, 5, 2, order=2)
-    wor_zero.run()
-    assert_almost_equal(wor_zero.timeseries[1], (0.0, 0.0, 0.0))
-
-def test_WaterOrientationalRelaxation_order_1(universe):
+def test_WaterOrientationalRelaxationRegion(universe):
     wor = waterdynamics.WaterOrientationalRelaxation(
-        universe, SELECTION1, 0, 5, 2, order=1)
-    wor.run()
-    assert_almost_equal(wor.timeseries[1][2], 0.71486,
+        universe.atoms, ["name OH2", "name OH2 and prop z <5 ", "name H1 or name H2"], 2, 1.2, order=1)
+    wor.run(start=0, stop=4)
+    assert_almost_equal(wor.results["correlation"][2],  0.810989,
                         decimal=5)
 
 
-def test_WaterOrientationalRelaxation_order_1_zeroMolecules(universe):
-    wor_zero = waterdynamics.WaterOrientationalRelaxation(
-        universe, SELECTION2, 0, 5, 2, order=1)
-    wor_zero.run()
-    assert_almost_equal(wor_zero.timeseries[1], (0.0, 0.0, 0.0))
+# def test_WaterOrientationalRelaxationDeprecated(universe):
+#     wor = waterdynamics.WaterOrientationalRelaxation(
+#         universe, SELECTION1, 0, 5, 2, order=2)
+#     wor.run()
+#     assert_almost_equal(wor.timeseries[1][2], 0.35887,
+#                         decimal=5)
+
+
+# def test_WaterOrientationalRelaxation_zeroMoleculesDeprecated(universe):
+#     wor_zero = waterdynamics.WaterOrientationalRelaxation(
+#         universe, SELECTION2, 0, 5, 2, order=2)
+#     wor_zero.run()
+#     assert_almost_equal(wor_zero.timeseries[1], (0.0, 0.0, 0.0))
+
+# def test_WaterOrientationalRelaxation_order_1Deprecated(universe):
+#     wor = waterdynamics.WaterOrientationalRelaxation(
+#         universe, SELECTION1, 0, 5, 2, order=1)
+#     wor.run()
+#     print(wor.timeseries)
+#     assert_almost_equal(wor.timeseries[1][2], 0.71486,
+#                         decimal=5)
+
+
+# def test_WaterOrientationalRelaxation_order_1_zeroMoleculesDeprecated(universe):
+#     wor_zero = waterdynamics.WaterOrientationalRelaxation(
+#         universe, SELECTION2, 0, 5, 2, order=1)
+#     wor_zero.run()
+#     assert_almost_equal(wor_zero.timeseries[1], (0.0, 0.0, 0.0))
 
 
 def test_AngularDistribution(universe):
@@ -99,7 +116,8 @@ def test_SurvivalProbability_intermittency1and2(universe):
     """
     with patch.object(universe, 'select_atoms') as select_atoms_mock:
         ids = [(9, 8), (), (8,), (9,), (8,), (), (9, 8), (), (8,), (9, 8)]
-        select_atoms_mock.side_effect = lambda selection: Mock(ids=ids.pop())   # atom IDs fed set by set
+        select_atoms_mock.side_effect = lambda selection: Mock(
+            ids=ids.pop())   # atom IDs fed set by set
         sp = waterdynamics.SurvivalProbability(universe, "")
         sp.run(tau_max=3, stop=10, verbose=True, intermittency=2)
         assert all(x == {9, 8} for x in sp._intermittent_selected_ids)
@@ -113,7 +131,8 @@ def test_SurvivalProbability_intermittency2lacking(universe):
     """
     with patch.object(universe, 'select_atoms') as select_atoms_mock:
         ids = [(9,), (), (), (), (9,), (), (), (), (9,)]
-        select_atoms_mock.side_effect = lambda selection: Mock(ids=ids.pop())   # atom IDs fed set by set
+        select_atoms_mock.side_effect = lambda selection: Mock(
+            ids=ids.pop())   # atom IDs fed set by set
         sp = waterdynamics.SurvivalProbability(universe, "")
         sp.run(tau_max=3, stop=9, verbose=True, intermittency=2)
         assert_almost_equal(sp.sp_timeseries, [1, 0, 0, 0])
@@ -125,8 +144,10 @@ def test_SurvivalProbability_intermittency1_step5_noSkipping(universe):
     No frames should be skipped.
     """
     with patch.object(universe, 'select_atoms') as select_atoms_mock:
-        ids = [(2, 3), (3,), (2, 3), (3,), (2,), (3,), (2, 3), (3,), (2, 3), (2, 3)]
-        select_atoms_mock.side_effect = lambda selection: Mock(ids=ids.pop())   # atom IDs fed set by set
+        ids = [(2, 3), (3,), (2, 3), (3,), (2,),
+               (3,), (2, 3), (3,), (2, 3), (2, 3)]
+        select_atoms_mock.side_effect = lambda selection: Mock(
+            ids=ids.pop())   # atom IDs fed set by set
         sp = waterdynamics.SurvivalProbability(universe, "")
         sp.run(tau_max=2, stop=10, verbose=True, intermittency=1, step=5)
         assert all((x == {2, 3} for x in sp._intermittent_selected_ids))
@@ -141,7 +162,8 @@ def test_SurvivalProbability_intermittency1_step5_Skipping(universe):
     with patch.object(universe, 'select_atoms') as select_atoms_mock:
         ids = [(1,), (), (1,), (), (1,), (), (1,), (), (1,), (1,)]
         beforepopsing = len(ids) - 2
-        select_atoms_mock.side_effect = lambda selection: Mock(ids=ids.pop())   # atom IDs fed set by set
+        select_atoms_mock.side_effect = lambda selection: Mock(
+            ids=ids.pop())   # atom IDs fed set by set
         sp = waterdynamics.SurvivalProbability(universe, "")
         sp.run(tau_max=1, stop=10, verbose=True, intermittency=1, step=5)
         assert all((x == {1} for x in sp._intermittent_selected_ids))
@@ -151,44 +173,50 @@ def test_SurvivalProbability_intermittency1_step5_Skipping(universe):
 
 def test_intermittency_none():
     # No changes asked - returns the same data
-    input_ids = [{1}, {1}, {1}, set(), set(), {1}, {1}, {1}, set(), set(), {1}, set(), set(), {1}]
+    input_ids = [{1}, {1}, {1}, set(), set(), {1}, {1}, {1}, set(),
+                 set(), {1}, set(), set(), {1}]
     corrected = correct_intermittency(input_ids, intermittency=0)
-    assert all(x == y for x,y in zip(input_ids, corrected))
+    assert all(x == y for x, y in zip(input_ids, corrected))
 
 
 def test_intermittency_1and2():
     # The maximum gap in the dataset is 2, so the IDs are always present after correction
-    input_ids = [{9, 8}, set(), {8, }, {9, }, {8, }, set(), {9, 8}, set(), {8, }, {9, 8, }]
+    input_ids = [{9, 8}, set(), {8, }, {9, }, {8, }, set(),
+                 {9, 8}, set(), {8, }, {9, 8, }]
     corrected = correct_intermittency(input_ids, intermittency=2)
     assert all((x == {9, 8} for x in corrected))
 
 
 def test_intermittency_2tooShort():
-    #The IDs are abscent for too long/
-    input_ids = [{9,}, {}, {}, {}, {9,}, {}, {}, {}, {9,}]
+    # The IDs are abscent for too long/
+    input_ids = [{9, }, {}, {}, {}, {9, }, {}, {}, {}, {9, }]
     corrected = correct_intermittency(input_ids, intermittency=2)
     assert all(x == y for x, y in zip(input_ids, corrected))
 
 
 def test_intermittency_setsOfSets():
     # Verificaiton for the case of hydrogen bonds (sets of sets)
-    input_ids = [{frozenset({1,2}), frozenset({3, 4})},set(), set(),
+    input_ids = [{frozenset({1, 2}), frozenset({3, 4})}, set(), set(),
                  {frozenset({1, 2}), frozenset({3, 4})}, set(), set(),
                  {frozenset({1, 2}), frozenset({3, 4})}, set(), set(),
                  {frozenset({1, 2}), frozenset({3, 4})}]
     corrected = correct_intermittency(input_ids, intermittency=2)
-    assert all((x == {frozenset({1, 2}), frozenset({3, 4})} for x in corrected))
+    assert all((x == {frozenset({1, 2}), frozenset({3, 4})}
+               for x in corrected))
 
 
 def test_autocorrelation_alwaysPresent():
     input = [{1, 2}, {1, 2}, {1, 2}, {1, 2}, {1, 2}, {1, 2}, {1, 2}]
-    tau_timeseries, sp_timeseries, sp_timeseries_data = autocorrelation(input, tau_max=3)
+    tau_timeseries, sp_timeseries, sp_timeseries_data = autocorrelation(
+        input, tau_max=3)
     assert all(np.equal(sp_timeseries, 1))
 
 
 def test_autocorrelation_definedTaus():
-    input_ids = [{9, 8, 7}, {8, 7, 6}, {7, 6, 5}, {6, 5, 4}, {5, 4, 3}, {4, 3, 2}, {3, 2, 1}]
-    tau_timeseries, sp_timeseries, sp_timeseries_data = autocorrelation(input_ids, tau_max=3)
+    input_ids = [{9, 8, 7}, {8, 7, 6}, {7, 6, 5}, {
+        6, 5, 4}, {5, 4, 3}, {4, 3, 2}, {3, 2, 1}]
+    tau_timeseries, sp_timeseries, sp_timeseries_data = autocorrelation(
+        input_ids, tau_max=3)
     assert_almost_equal(sp_timeseries, [1, 2/3., 1/3., 0])
 
 
@@ -197,7 +225,8 @@ def test_autocorrelation_intermittency1_windowJump_intermittencyAll():
     Step leads to skipping frames if (tau_max + 1) + (intermittency * 2) < step.
     No frames should be skipped so intermittency should be applied to all.
     """
-    input_ids = [{2, 3}, {3,}, {2, 3}, {3,}, {2,}, {3,}, {2, 3}, {3,}, {2, 3}, {2, 3}]
+    input_ids = [{2, 3}, {3, }, {2, 3}, {3, }, {
+        2, }, {3, }, {2, 3}, {3, }, {2, 3}, {2, 3}]
     corrected = correct_intermittency(input_ids, intermittency=1)
     tau_timeseries, sp_timeseries, sp_timeseries_data = autocorrelation(corrected, tau_max=2,
                                                                         window_step=5)
@@ -206,21 +235,26 @@ def test_autocorrelation_intermittency1_windowJump_intermittencyAll():
 
 
 def test_autocorrelation_windowBigJump():
-    #The empty sets are ignored (no intermittency)
-    input_ids = [{1}, {1}, {1}, set(), set(), {1}, {1}, {1}, set(), set(), {1}, {1}, {1}]
-    tau_timeseries, sp_timeseries, sp_timeseries_data = autocorrelation(input_ids, tau_max=2, window_step=5)
+    # The empty sets are ignored (no intermittency)
+    input_ids = [{1}, {1}, {1}, set(), set(), {1}, {1}, {1},
+                 set(), set(), {1}, {1}, {1}]
+    tau_timeseries, sp_timeseries, sp_timeseries_data = autocorrelation(
+        input_ids, tau_max=2, window_step=5)
     assert_almost_equal(sp_timeseries, [1, 1, 1])
 
 
 def test_autocorrelation_windowBigJump_absence():
     # In the last frame the molecules are absent
-    input_ids = [{1}, {1}, {1}, set(), set(), {1}, {1}, {1}, set(), set(), {1}, set(), set()]
-    tau_timeseries, sp_timeseries, sp_timeseries_data = autocorrelation(input_ids, tau_max=2, window_step=5)
+    input_ids = [{1}, {1}, {1}, set(), set(), {1}, {1}, {1},
+                 set(), set(), {1}, set(), set()]
+    tau_timeseries, sp_timeseries, sp_timeseries_data = autocorrelation(
+        input_ids, tau_max=2, window_step=5)
     assert_almost_equal(sp_timeseries, [1, 2/3., 2/3.])
 
 
 def test_autocorrelation_intermittency1_many():
-    input_ids = [{1}, set(), {1}, set(), {1}, set(), {1}, set(), {1}, set(), {1}, set(), {1}, set(), {1}]
+    input_ids = [{1}, set(), {1}, set(), {1}, set(), {1}, set(), {
+        1}, set(), {1}, set(), {1}, set(), {1}]
     corrected = correct_intermittency(input_ids, intermittency=1)
     tau_timeseries, sp_timeseries, sp_timeseries_data = autocorrelation(corrected, tau_max=14,
                                                                         window_step=5)
@@ -229,7 +263,8 @@ def test_autocorrelation_intermittency1_many():
 
 def test_autocorrelation_intermittency2_windowBigJump():
     # The intermittency corrects the last frame
-    input_ids = [{1}, {1}, {1}, set(), set(), {1}, {1}, {1}, set(), set(), {1}, set(), set(), {1}]
+    input_ids = [{1}, {1}, {1}, set(), set(), {1}, {1}, {1}, set(),
+                 set(), {1}, set(), set(), {1}]
     corrected = correct_intermittency(input_ids, intermittency=2)
     tau_timeseries, sp_timeseries, sp_timeseries_data = autocorrelation(corrected, tau_max=2,
                                                                         window_step=5)
@@ -238,8 +273,10 @@ def test_autocorrelation_intermittency2_windowBigJump():
 
 def test_SurvivalProbability_t0tf(universe):
     with patch.object(universe, 'select_atoms') as select_atoms_mock:
-        ids = [(0, ), (0, ), (7, 6, 5), (6, 5, 4), (5, 4, 3), (4, 3, 2), (3, 2, 1), (0, )]
-        select_atoms_mock.side_effect = lambda selection: Mock(ids=ids.pop(2))   # atom IDs fed set by set
+        ids = [(0, ), (0, ), (7, 6, 5), (6, 5, 4),
+               (5, 4, 3), (4, 3, 2), (3, 2, 1), (0, )]
+        select_atoms_mock.side_effect = lambda selection: Mock(
+            ids=ids.pop(2))   # atom IDs fed set by set
         sp = waterdynamics.SurvivalProbability(universe, "")
         sp.run(tau_max=3, start=2, stop=7)
         assert_almost_equal(sp.sp_timeseries, [1, 2 / 3.0, 1 / 3.0, 0])
@@ -247,8 +284,10 @@ def test_SurvivalProbability_t0tf(universe):
 
 def test_SurvivalProbability_definedTaus(universe):
     with patch.object(universe, 'select_atoms') as select_atoms_mock:
-        ids = [(9, 8, 7), (8, 7, 6), (7, 6, 5), (6, 5, 4), (5, 4, 3), (4, 3, 2), (3, 2, 1)]
-        select_atoms_mock.side_effect = lambda selection: Mock(ids=ids.pop())   # atom IDs fed set by set
+        ids = [(9, 8, 7), (8, 7, 6), (7, 6, 5), (6, 5, 4),
+               (5, 4, 3), (4, 3, 2), (3, 2, 1)]
+        select_atoms_mock.side_effect = lambda selection: Mock(
+            ids=ids.pop())   # atom IDs fed set by set
         sp = waterdynamics.SurvivalProbability(universe, "")
         sp.run(tau_max=3, start=0, stop=7, verbose=True)
         assert_almost_equal(sp.sp_timeseries, [1, 2 / 3.0, 1 / 3.0, 0])
diff --git a/waterdynamics/waterdynamics.py b/waterdynamics/waterdynamics.py
index 4e9f75a..efee77d 100644
--- a/waterdynamics/waterdynamics.py
+++ b/waterdynamics/waterdynamics.py
@@ -39,7 +39,10 @@
 
 .. footbibliography::
 """
+from MDAnalysis.lib.log import ProgressBar
 from MDAnalysis.lib.correlations import autocorrelation, correct_intermittency
+from MDAnalysis.lib.distances import minimize_vectors
+from MDAnalysis.analysis.base import AnalysisBase
 from itertools import zip_longest
 import logging
 import warnings
@@ -47,10 +50,179 @@
 
 
 logger = logging.getLogger('waterdynamics')
-from MDAnalysis.lib.log import ProgressBar
 
 
-class WaterOrientationalRelaxation(object):
+class WaterOrientationalRelaxation(AnalysisBase):
+    r"""Water orientation relaxation analysis
+
+    Function to evaluate the Water Orientational Relaxation proposed by Yu-ling
+    Yeh and Chung-Yuan Mou :footcite:p:`Yeh1999`. WaterOrientationalRelaxation
+    indicates "how fast" water molecules are rotating or changing direction.
+    This is a time correlation function given by:
+
+    .. math::
+        C_{\hat u}(\tau)=\langle \mathit{P}_2[\mathbf{\hat{u}}(t_0)\cdot\mathbf{\hat{u}}(t_0+\tau)]\rangle
+
+    where :math:`P_2=(3x^2-1)/2` is the second-order Legendre polynomial and :math:`\hat{u}` is
+    a unit vector along HH, OH or dipole vector. Another option is to select the first-order Legendre
+    polynomial, :math:`P_1=x`.
+
+
+    Parameters
+    ----------
+    atomgroup : AtomGroup
+      AtomGroup object
+    selection : str
+      Selection is a list of three strings,
+      The first one is the selection of all relevant oxygen atoms that are possible to be the atoms of interest
+      The second one is the selection of the oxygen atoms of interest, for example, atoms in a region.
+      The third one is the selection of the hydrogen atoms that are possible to attach to the oxygen atoms of interest.
+      Example: ["name OH2", "name OH2 and prop z <5 " ,"name H1 or name H2"]
+    lag_max : int
+      Maximum lag time for correaltion
+    OH_cutoff : float
+        Cutoff distance for the hydrogen atoms to be considered as attached to the oxygen atoms of interest.
+    update_selections: bool
+      If True, the selection of the oxygen atoms of interest and
+      the hydrogen atoms that are possible to attach to the oxygen atoms of interest will be updated at each frame.
+    order : 1 or 2 (default)
+      first- or second-order Legendre polynomial
+    """
+
+    def __init__(self, atomgroup, selection, lag_max, OH_cutoff,
+                 nproc=1, order=2, update_selections=False, **kwargs):
+        super(WaterOrientationalRelaxation, self).__init__(atomgroup.universe.trajectory,
+                                                           **kwargs)
+        # selection is a list of two strings, one for all O, and for subsect,
+        # one for H
+
+        self._ag = atomgroup
+        # selection of all relevant oxygen atoms that are possible to be the
+        # atoms of interest
+        self.selection_oxygen = selection[0]
+        # selection of the oxygen atoms of interest, for example, atoms in a region.
+        # The oxygen atoms can flow in and out of this region.
+        self.selection_oxygen_subset = selection[1]
+        # selection of the hydrogen atoms that are possible to attach to the
+        # oxygen atoms of interest
+        self.selection_hydrogen = selection[2]
+        self.update_selections = update_selections
+
+        self.OH_cutoff = OH_cutoff
+        self.lag_max = lag_max
+
+        # parallel not yet implemented
+        self.nproc = nproc
+
+        if order != 1 and order != 2:
+            raise ValueError(
+                f"order = {order} but only first- or second-order Legendre polynomial is allowed.")
+        else:
+            self.order = order
+
+    def _prepare(self):
+        # Called before iteration on the trajectory has begun.
+        # TODO: parallel computation.
+        # use lagseries instead of timeseries to avoid confusion with the
+        # quantity over timestep.
+        self.results = {}
+        self.cellpar = self._ag.dimensions
+
+        self._ag_oxygen = self._ag.select_atoms(self.selection_oxygen)
+        self._ag_oxygen_subset = self._ag.select_atoms(
+            self.selection_oxygen_subset)
+        self.n_oxygen = len(self._ag_oxygen)
+
+        self.get_water_list()
+        # matrix to store the dipole vectors, n_frames x n_oxygen x 3
+        self.dipole_matrix = np.empty((self.n_frames, self.n_oxygen, 3))
+        # np.nan is a placeholder for missing data
+        self.dipole_matrix.fill(np.nan)
+
+    def _single_frame(self):
+        # loop over frames
+        if self.update_selections:
+            self.get_water_list()
+
+        for count, idxs_water in enumerate(self.idxs_water_list):
+            if len(idxs_water) == 0:
+                continue
+            pos_O = self._ag[idxs_water[0]].position
+            pos_H1 = self._ag[idxs_water[1]].position
+            pos_H2 = self._ag[idxs_water[2]].position
+
+            # compute dipole vector
+            OHVector1 = pos_H1 - pos_O
+            OHVector2 = pos_H2 - pos_O
+            # consider the periodic boundary conditions
+            vectors = minimize_vectors(
+                np.array([OHVector1, OHVector2]), box=self.cellpar)
+            dipVector = np.mean(vectors, axis=0)
+
+            # normalize the dipole vector
+            normdipVector = np.linalg.norm(dipVector)
+            unitdipVector = dipVector / normdipVector
+
+            self.dipole_matrix[self._frame_index, count] = unitdipVector
+
+    def _conclude(self):
+
+        # lagseries stores the correlation (lag) time
+        # correlation stores the corresponding correlation value
+        self.results["lagseries"] = np.zeros((self.lag_max + 1))
+        self.results["lagseries"][0] = 0  # first lag time is zero.
+        self.results["correlation"] = np.zeros((self.lag_max + 1))
+        self.results["correlation"][0] = 1.0  # first correlation must be 1.0
+
+        # compute the correlation
+        # diole matrix has the shape of (n_frames, n_oxygen, 3)
+        # print(self.dipole_matrix)
+        for lag in range(1, self.lag_max + 1):
+            # first we slice the dipole matrix according to the lag time
+            reduced_dipole_matrix = self.dipole_matrix[::lag]
+            # then we compute the correlation by shifting one frame
+            correlation = np.sum(
+                reduced_dipole_matrix[1:] * reduced_dipole_matrix[:-1], axis=2)
+            if self.order == 1:
+                correlation = self.lg1(correlation)
+            elif self.order == 2:
+                correlation = self.lg2(correlation)
+            # nan left in the array, use nan mean to exclude them.
+            correlation = np.nanmean(correlation)
+
+            self.results["lagseries"][lag] = lag
+            self.results["correlation"][lag] = correlation
+
+    def get_water_list(self):
+        self.idxs_water_list = []
+        # print("List of found water molecules")
+        # print("Oxygen atom, Hydrogen atom 1, Hydrogen atom 2")
+        for i in range(self.n_oxygen):
+            if self._ag_oxygen[i].id in self._ag_oxygen_subset.ids:
+                self._tmp_ag_H = self._ag.select_atoms(
+                    f" ({self.selection_hydrogen}) and around {self.OH_cutoff} id {self._ag_oxygen[i].id}")
+                # if not a water molecule, append place holder [] and skip
+                if len(self._tmp_ag_H) != 2:
+                    self.idxs_water_list.append([])
+                    continue
+                idxs_water = [idx for idx in self._tmp_ag_H.ids]
+                idxs_water.insert(0, self._ag_oxygen[i].id)
+                self.idxs_water_list.append(idxs_water)
+            else:
+                self.idxs_water_list.append([])
+
+    @staticmethod
+    def lg1(x):
+        """First Legendre polynomial"""
+        return x
+
+    @staticmethod
+    def lg2(x):
+        """Second Legendre polynomial"""
+        return (3 * x * x - 1) / 2
+
+
+class WaterOrientationalRelaxationDeprecated(object):
     r"""Water orientation relaxation analysis
 
     Function to evaluate the Water Orientational Relaxation proposed by Yu-ling
@@ -90,7 +262,8 @@ def __init__(self, universe, select, t0, tf, dtmax, nproc=1, order=2):
         self.dtmax = dtmax
         self.nproc = nproc
         if order != 1 and order != 2:
-            raise ValueError(f"order = {order} but only first- or second-order Legendre polynomial is allowed.")
+            raise ValueError(
+                f"order = {order} but only first- or second-order Legendre polynomial is allowed.")
         else:
             self.order = order
         self.timeseries = None
@@ -177,8 +350,7 @@ def _getOneDeltaPoint(self, universe, repInd, i, t0, dt):
                 valHH += self.lg2(np.dot(unitHHVector0, unitHHVectorp))
                 valdip += self.lg2(np.dot(unitdipVector0, unitdipVectorp))
             n += 1
-        return  (valOH/n, valHH/n, valdip/n) if n > 0 else (0, 0, 0)
-
+        return (valOH / n, valHH / n, valdip / n) if n > 0 else (0, 0, 0)
 
     def _getMeanOnePoint(self, universe, selection1, selection_str, dt,
                          totalFrames):
@@ -203,7 +375,8 @@ def _getMeanOnePoint(self, universe, selection1, selection_str, dt,
 
         # if no water molecules remain in selection, there is nothing to get
         # the mean, so n = 0.
-        return (sumDeltaOH / n, sumDeltaHH / n, sumDeltadip / n) if n > 0 else (0, 0, 0)
+        return (sumDeltaOH / n, sumDeltaHH / n,
+                sumDeltadip / n) if n > 0 else (0, 0, 0)
 
     def _sameMolecTandDT(self, selection, t0d, tf):
         """
@@ -233,7 +406,7 @@ def lg1(x):
     @staticmethod
     def lg2(x):
         """Second Legendre polynomial"""
-        return (3*x*x - 1)/2
+        return (3 * x * x - 1) / 2
 
     def run(self, **kwargs):
         """Analyze trajectory and produce timeseries"""
@@ -486,7 +659,7 @@ def _getOneDeltaPoint(self, universe, repInd, i, t0, dt):
 
         # if no water molecules remain in selection, there is nothing to get
         # the mean, so n = 0.
-        return valO/n if n > 0 else 0
+        return valO / n if n > 0 else 0
 
     def _getMeanOnePoint(self, universe, selection1, selection_str, dt,
                          totalFrames):
@@ -509,7 +682,7 @@ def _getMeanOnePoint(self, universe, selection1, selection_str, dt,
 
         # if no water molecules remain in selection, there is nothing to get
         # the mean, so n = 0.
-        return sumDeltaO/n if n > 0 else 0
+        return sumDeltaO / n if n > 0 else 0
 
     def _sameMolecTandDT(self, selection, t0d, tf):
         """
@@ -558,9 +731,9 @@ class SurvivalProbability(object):
     .. math::
         P(\tau) = \langle \frac{ N(t, t + \tau )} { N(t) }\rangle
 
-    where :math:`\tau` is the timestep, :math:`N(t)` the number of particles at time 
-    :math:`t`, and :math:`N(t, t+\tau)` is the number of particles at every frame from 
-    :math:`t` to :math:`t + \tau`. The angular brackets represent an average over all time 
+    where :math:`\tau` is the timestep, :math:`N(t)` the number of particles at time
+    :math:`t`, and :math:`N(t, t+\tau)` is the number of particles at every frame from
+    :math:`t` to :math:`t + \tau`. The angular brackets represent an average over all time
     origins, :math:`t`. See :func:`MDAnalysis.lib.correlations.autocorrelation` for
     technical details.
 
@@ -571,7 +744,7 @@ class SurvivalProbability(object):
       Universe object
     select : str
       Selection string; any selection is allowed. With this selection you
-      define the region/zone where to analyze, e.g.: "resname SOL and around 5 
+      define the region/zone where to analyze, e.g.: "resname SOL and around 5
       (resid 10)".
     verbose : Boolean, optional
       When True, prints progress and comments to the console.
@@ -698,13 +871,15 @@ def run(self, tau_max=20, start=None, stop=None, step=None, residues=False,
         # and for extra frames that were loaded (intermittency)
         window_jump = step - num_frames_to_skip
 
-        self._intermittent_selected_ids = correct_intermittency(self._selected_ids, intermittency=intermittency)
+        self._intermittent_selected_ids = correct_intermittency(
+            self._selected_ids, intermittency=intermittency)
         tau_timeseries, sp_timeseries, sp_timeseries_data = autocorrelation(self._intermittent_selected_ids,
                                                                             tau_max, window_jump)
 
         # warn the user if the NaN are found
         if all(np.isnan(sp_timeseries[1:])):
-            logger.warning('NaN Error: Most likely data was not found. Check your atom selections. ')
+            logger.warning(
+                'NaN Error: Most likely data was not found. Check your atom selections. ')
 
         # user can investigate the distribution and sample size
         self.sp_timeseries_data = sp_timeseries_data