Source code for MDAnalysis.analysis.rdf

# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*-
# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
#
# MDAnalysis --- http://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.
#
# 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
#

"""
Radial Distribution Functions --- :mod:`MDAnalysis.analysis.rdf`
================================================================

Tools for calculating pair distribution functions ("radial
distribution functions" or "RDF").

.. Not Implemented yet:
.. - Structure factor?
.. - Coordination number

"""
from __future__ import division, absolute_import
import numpy as np

from ..lib.util import blocks_of
from ..lib import distances
from .base import AnalysisBase


[docs]class InterRDF(AnalysisBase): """Intermolecular pair distribution function InterRDF(g1, g2, nbins=75, range=(0.0, 15.0)) Arguments --------- g1 : AtomGroup First AtomGroup g2 : AtomGroup Second AtomGroup nbins : int (optional) Number of bins in the histogram [75] range : tuple or list (optional) The size of the RDF [0.0, 15.0] exclusion_block : tuple (optional) A tuple representing the tile to exclude from the distance array. [None] start : int (optional) The frame to start at (default is first) stop : int (optional) The frame to end at (default is last) step : int (optional) The step size through the trajectory in frames (default is every frame) Example ------- First create the :class:`InterRDF` object, by supplying two AtomGroups then use the :meth:`run` method :: rdf = InterRDF(ag1, ag2) rdf.run() Results are available through the :attr:`bins` and :attr:`rdf` attributes:: plt.plot(rdf.bins, rdf.rdf) The `exclusion_block` keyword allows the masking of pairs from within the same molecule. For example, if there are 7 of each atom in each molecule, the exclusion mask `(7, 7)` can be used. .. versionadded:: 0.13.0 """ def __init__(self, g1, g2, nbins=75, range=(0.0, 15.0), exclusion_block=None, **kwargs): super(InterRDF, self).__init__(g1.universe.trajectory, **kwargs) self.g1 = g1 self.g2 = g2 self.u = g1.universe self.rdf_settings = {'bins': nbins, 'range': range} self._exclusion_block = exclusion_block def _prepare(self): # Empty histogram to store the RDF count, edges = np.histogram([-1], **self.rdf_settings) count = count.astype(np.float64) count *= 0.0 self.count = count self.edges = edges self.bins = 0.5 * (edges[:-1] + edges[1:]) # Need to know average volume self.volume = 0.0 # Allocate a results array which we will reuse self._result = np.zeros((len(self.g1), len(self.g2)), dtype=np.float64) # If provided exclusions, create a mask of _result which # lets us take these out if self._exclusion_block is not None: self._exclusion_mask = blocks_of(self._result, *self._exclusion_block) self._maxrange = self.rdf_settings['range'][1] + 1.0 else: self._exclusion_mask = None def _single_frame(self): distances.distance_array(self.g1.positions, self.g2.positions, box=self.u.dimensions, result=self._result) # Maybe exclude same molecule distances if self._exclusion_mask is not None: self._exclusion_mask[:] = self._maxrange count = np.histogram(self._result, **self.rdf_settings)[0] self.count += count self.volume += self._ts.volume def _conclude(self): # Number of each selection nA = len(self.g1) nB = len(self.g2) N = nA * nB # If we had exclusions, take these into account if self._exclusion_block: xA, xB = self._exclusion_block nblocks = nA / xA N -= xA * xB * nblocks # Volume in each radial shell vol = np.power(self.edges[1:], 3) - np.power(self.edges[:-1], 3) vol *= 4/3.0 * np.pi # Average number density box_vol = self.volume / self.n_frames density = N / box_vol rdf = self.count / (density * vol * self.n_frames) self.rdf = rdf