diff --git a/package/AUTHORS b/package/AUTHORS index 5b347d08539..3788d51045f 100644 --- a/package/AUTHORS +++ b/package/AUTHORS @@ -206,6 +206,7 @@ Chronological list of authors - Christian Pfaendner - Pratham Chauhan - Meet Brijwani + - Pegerto Fernandez External code diff --git a/package/CHANGELOG b/package/CHANGELOG index 0a34834cf4e..3fe5594b0b1 100644 --- a/package/CHANGELOG +++ b/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 diff --git a/package/MDAnalysis/analysis/accesibility.py b/package/MDAnalysis/analysis/accesibility.py new file mode 100644 index 00000000000..0b38192d98c --- /dev/null +++ b/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() diff --git a/package/doc/sphinx/source/documentation_pages/analysis/accesibility.rst b/package/doc/sphinx/source/documentation_pages/analysis/accesibility.rst new file mode 100644 index 00000000000..41a9aa36077 --- /dev/null +++ b/package/doc/sphinx/source/documentation_pages/analysis/accesibility.rst @@ -0,0 +1,2 @@ +.. automodule:: MDAnalysis.analysis.accesibility + diff --git a/testsuite/MDAnalysisTests/analysis/test_accesibility.py b/testsuite/MDAnalysisTests/analysis/test_accesibility.py new file mode 100644 index 00000000000..f0972787e8c --- /dev/null +++ b/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))