Solvation Shell

doc/sphinx/source/analysis.txt

@@ -17,4 +17,5 @@ are for users who wish to construct their own analyses.
 
    analysis/ensemble
    analysis/ensemble_analysis
-   analysis/dihedral
+   analysis/solvation
+   analysis/dihedral

doc/sphinx/source/analysis/solvation.txt

@@ -0,0 +1,12 @@
+========================
+Solvation Shell Analysis
+========================
+
+Analyzes the number of solvent molecules within given distances of the solute.
+
+.. versionadded:: 0.8.0
+
+.. autoclass:: mdpow.analysis.solvation.SolvationAnalysis
+    :members:
+
+    .. automethod:: run

mdpow/analysis/solvation.py

@@ -0,0 +1,85 @@
+# MDPOW: solvation.py
+# 2021 Alia Lescoulie
+
+from typing import List
+
+import numpy as np
+import pandas as pd
+
+import MDAnalysis as mda
+from MDAnalysis.lib.distances import capped_distance
+
+from .ensemble import EnsembleAnalysis, Ensemble, EnsembleAtomGroup
+
+import logging
+
+logger = logging.getLogger('mdpow.dihedral')
+
+
+class SolvationAnalysis(EnsembleAnalysis):
+    """Measures the number of solvent molecules withing the given distances
+    in an :class:`~mdpow.analysis.ensemble.Ensemble` .
+
+    :keyword:
+
+    *solute*
+        An :class:`~mdpow.analysis.ensemble.EnsembleAtom` containing the solute
+        used to measure distance.
+
+    *solvent*
+        An :class:`~mdpow.analysis.ensemble.EnsembleAtom` containing the solvents
+        counted in by the distance measurement. Each solvent atom is counted by the
+        distance calculation.
+
+
+    *distances*
+        The cutoff distances around the solute measured in Angstroms.
+
+    The data is returned in a :class:`pandas.DataFrame` with observations sorted by
+    distance, solvent, interaction, lambda, time.
+
+    .. ruberic:: Example
+
+    Typical Workflow::
+
+        ens = Ensemble(dirname='Mol')
+        solvent = ens.select_atoms('resname SOL and name OW')
+        solute = ens.select_atoms('not resname SOL')
+
+        solv_dist = SolvationAnalysis(solute, solvent, [1.2, 2.4]).run(start=0, stop=10, step=1)
+
+    """
+    def __init__(self, solute: EnsembleAtomGroup, solvent: EnsembleAtomGroup, distances: List[float]):
+        self.check_groups_from_common_ensemble([solute, solvent])
+        super(SolvationAnalysis, self).__init__(solute.ensemble)
+        self._solute = solute
+        self._solvent = solvent
+        self._dists = distances
+
+    def _prepare_ensemble(self):
+        self._sel = 'name '
+        keys = [k for k in self._solute.keys()]
+        for n in self._solute[keys[0]].names:
+            self._sel += f' {n}'
+        self._col = ['distance', 'solvent', 'interaction',
+                     'lambda', 'time', 'N_solvent']
+        self.results = pd.DataFrame(columns=self._col)
+        self._res_dict = {key: [] for key in self._col}
+
+    def _single_frame(self):
+        solute = self._solute[self._key]
+        solvent = self._solvent[self._key]
+        pairs, distaces = capped_distance(solute.positions, solvent.positions,
+                                          max(self._dists), box=self._ts.dimensions)
+        solute_i, solvent_j = np.transpose(pairs)
+        for d in range(len(self._dists)):
+            close_solv_atoms = solvent[solvent_j[distaces < self._dists[d]]]
+            n = len(close_solv_atoms)
+            result = [self._dists[d], self._key[0], self._key[1],
+                      self._key[2], self._ts.time, n]
+            for i in range(len(self._col)):
+                self._res_dict[self._col[i]].append(result[i])
+
+    def _conclude_ensemble(self):
+        for k in self._col:
+            self.results[k] = self._res_dict[k]

mdpow/tests/test_solv_shell.py

@@ -0,0 +1,55 @@
+from __future__ import absolute_import
+
+import numpy as np
+
+from . import tempdir as td
+
+import py.path
+
+import pybol
+import pytest
+
+from numpy.testing import assert_almost_equal
+from scipy.stats import variation
+
+from ..analysis.ensemble import Ensemble, EnsembleAnalysis, EnsembleAtomGroup
+
+from ..analysis.solvation import SolvationAnalysis
+
+from pkg_resources import resource_filename
+
+RESOURCES = py.path.local(resource_filename(__name__, 'testing_resources'))
+MANIFEST = RESOURCES.join("manifest.yml")
+
+
+class TestSolvShell(object):
+    mean = 2654.0
+    std = 2654.465059103246
+
+    def setup(self):
+        self.tmpdir = td.TempDir()
+        self.m = pybol.Manifest(str(RESOURCES / 'manifest.yml'))
+        self.m.assemble('example_FEP', self.tmpdir.name)
+        self.ens = Ensemble(dirname=self.tmpdir.name, solvents=['water'])
+        self.solute = self.ens.select_atoms('not resname SOL')
+        self.solvent = self.ens.select_atoms('resname SOL and name OW')
+
+    def teardown(self):
+        self.tmpdir.dissolve()
+
+    def test_dataframe(self):
+        solv = SolvationAnalysis(self.solute, self.solvent, [1.2]).run(start=0, stop=4, step=1)
+
+        for d in solv.results['distance']:
+            assert d == 1.2
+        for s in solv.results['solvent']:
+            assert s == 'water'
+        for i in solv.results['interaction'][:12]:
+            assert i == 'Coulomb'
+
+    def test_selection(self):
+        solv = SolvationAnalysis(self.solute, self.solvent, [2, 10]).run(start=0, stop=4, step=1)
+        mean = np.mean(solv.results['N_solvent'])
+        std = np.std(solv.results['N_solvent'])
+        assert_almost_equal(mean, self.mean, 6)
+        assert_almost_equal(std, self.std, 6)