Solvant Accesible surface area

package/AUTHORS

@@ -206,6 +206,7 @@ Chronological list of authors
   - Christian Pfaendner
   - Pratham Chauhan
   - Meet Brijwani
+  - Pegerto Fernandez
 
 
 External code

package/CHANGELOG

@@ -31,6 +31,8 @@ Enhancements
   * AuxReaders are now pickle-able and copy-able (Issue #1785, PR #3887)
   * Add pickling support for Atom, Residue, Segment, ResidueGroup 
     and SegmentGroup. (PR #3953)
+  * Add calculation for the solvant accessible area, using Shrake-Rupley
+    algorithm (PR #4025)
 
 Changes
   * As per NEP29 the minimum supported NumPy version has been raised to 1.21

package/MDAnalysis/analysis/accesibility.py

@@ -0,0 +1,126 @@
+# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*-
+# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
+#
+# MDAnalysis --- https://www.mdanalysis.org
+# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors
+# (see the file AUTHORS for the full list of names)
+#
+# Released under the GNU Public Licence, v2 or any higher version
+#
+# Please cite your use of MDAnalysis in published work:
+#
+# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler,
+# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein.
+# MDAnalysis: A Python package for the rapid analysis of molecular dynamics
+# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th
+# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy.
+# doi: 10.25080/majora-629e541a-00e
+#
+# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein.
+# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations.
+# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
+#
+
+import numpy as np
+
+from scipy.spatial import KDTree
+from numpy import arange, pi, sin, cos, arccos
+from .base import AnalysisBase
+
+
+ATOMIC_RADI = {
+    "H": 1.200,
+    "HE": 1.400,
+    "C": 1.700,
+    "N": 1.550,
+    "O": 1.520,
+    "F": 1.470,
+    "NA": 2.270,
+    "MG": 1.730,
+    "P": 1.800,
+    "S": 1.800,
+    "CL": 1.750,
+    "K": 2.750,
+    "CA": 2.310,
+    "NI": 1.630,
+    "CU": 1.400,
+    "ZN": 1.390,
+    "SE": 1.900,
+    "BR": 1.850,
+    "CD": 1.580,
+    "I": 1.980,
+    "HG": 1.550,
+}
+
+
+def atomic_radi(e):
+    """van der Waals radii"""
+    DEFAULT_ATOMIC_RADI = 2.0
+    return ATOMIC_RADI[e] if e in ATOMIC_RADI else DEFAULT_ATOMIC_RADI
+
+
+class SASA(AnalysisBase):
+    r"""Calculation of solvent accesible surface areas using Shrake-Rupley
+    'rolling ball' algorithm
+    """
+
+    def __init__(self, ag, n_points=100, probe_radius=1.40, **kwargs):
+        """Parameters
+        ----------
+        ag : AtomGroup
+            Atom group used for the calculation.
+        n_points: int (optional)
+            Number of points used for the estimation of the atom surface area
+        probe_radius: double ( optional )
+            Radius of the probe atom, by default watter radious
+        """
+        super(SASA, self).__init__(ag.universe.trajectory, **kwargs)
+        self._ag = ag
+        self._atom_neighbours = KDTree(ag.positions, 10)
+        self._n_points = n_points
+        self._sphere = self._compute_sphere()
+        self._twice_max_radi = (max(list(ATOMIC_RADI.values())) + probe_radius) * 2
+        self._probe_radius = probe_radius
+
+    def _compute_sphere(self):
+        """Fibonacci lattice to evently distribute points in the sphere."""
+        golden_ratio = (1 + 5**0.5) / 2
+        i = arange(0, self._n_points)
+        theta = 2 * pi * i / golden_ratio
+        phi = arccos(1 - 2 * (i + 0.5) / self._n_points)
+        x_points, y_points, z_points = (
+            cos(theta) * sin(phi),
+            sin(theta) * sin(phi),
+            cos(phi),
+        )
+        return np.transpose(np.array([x_points, y_points, z_points]))
+
+    def _compute(self):
+        asas = np.zeros(len(self._ag.atoms), dtype=np.double)
+
+        for atom in self._ag.atoms:
+            r_w = atomic_radi(atom.type)
+            raddi = r_w + self._probe_radius
+            # translate and transform the sphere to the centroid
+            sphere = np.array(self._sphere, copy=True) * raddi + atom.position
+            neighbours = self._atom_neighbours.query_ball_point(
+                atom.position, self._twice_max_radi
+            )
+            kdt_sphere = KDTree(sphere, 10)
+            intersect = set()
+            for n in neighbours:
+                if n == atom.index:
+                    continue
+                n_raddi = atomic_radi(self._ag.atoms[n].type) + self._probe_radius
+                cut = kdt_sphere.query_ball_point(self._ag.atoms[n].position, n_raddi)
+                intersect |= set(cut)
+
+            points = self._n_points - len(intersect)
+            total_surface_area = raddi * raddi * 4 * np.pi
+            area_per_point = total_surface_area / self._n_points
+            asas[atom.index] = points * area_per_point
+        return asas
+
+    def atoms(self):
+        """returns accessible surface area in A**2"""
+        return self._compute()

package/doc/sphinx/source/documentation_pages/analysis/accesibility.rst

@@ -0,0 +1,2 @@
+.. automodule:: MDAnalysis.analysis.accesibility
+

testsuite/MDAnalysisTests/analysis/test_accesibility.py

@@ -0,0 +1,43 @@
+# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*-
+# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 fileencoding=utf-8
+#
+# MDAnalysis --- https://www.mdanalysis.org
+# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors
+# (see the file AUTHORS for the full list of names)
+#
+# Released under the GNU Public Licence, v2 or any higher version
+#
+# Please cite your use of MDAnalysis in published work:
+#
+# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler,
+# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein.
+# MDAnalysis: A Python package for the rapid analysis of molecular dynamics
+# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th
+# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy.
+# doi: 10.25080/majora-629e541a-00e
+#
+# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein.
+# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations.
+# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
+#
+import pytest
+import MDAnalysis as mda
+
+from MDAnalysisTests.datafiles import PSF, DCD
+from MDAnalysis.analysis.accesibility import SASA
+
+
+class TestSASA(object):
+    @pytest.fixture()
+    def universe(self):
+        return mda.Universe(PSF, DCD)
+
+    def test_len_result_equal_natoms(self, universe):
+        sasa = SASA(universe.atoms)
+        atoms_accesible_area = sasa.atoms()
+        assert len(atoms_accesible_area) == len(universe.atoms)
+
+    def test_len_atom_area_gt_0(self, universe):
+        sasa = SASA(universe.atoms)
+        atoms_accesible_area = sasa.atoms()
+        assert all(map(lambda x: x >= 0, atoms_accesible_area))