Analysis - Linear Density module

package/AUTHORS

@@ -66,6 +66,7 @@ Chronological list of authors
   - Hai Nguyen
 2016
   - Balasubramanian
+  - Mattia F. Palermo
 
 External code
 -------------

package/CHANGELOG

@@ -31,6 +31,9 @@ Enhancement
   * select_atoms now implicitly ORs multiple values after a keyword for
     many types of selections (Issue #345)
   * Performance improvements for the PDBReader of about 10%
+  * LinearDensity analysis module added, which allows to compute linear mass
+    and charge density profiles along the three cartesian axes of the cell.
+    Works for orthorombic, fixed volume cells only. (Issue #670)
 
 Changes
 

package/MDAnalysis/analysis/lineardensity.py

@@ -0,0 +1,256 @@
+# -*- 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-2015 Naveen Michaud-Agrawal, Elizabeth J. Denning, Oliver Beckstein
+# and contributors (see AUTHORS for the full list)
+#
+# Released under the GNU Public Licence, v2 or any higher version
+#
+# Please cite your use of MDAnalysis in published work:
+#
+# 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
+#
+
+"""
+Linear Density --- :mod:`MDAnalysis.analysis.lineardensity`
+===========================================================
+
+A tool to compute mass and charge density profiles along the three
+cartesian axes of the simulation cell. Works only for orthorombic,
+fixed volume cells (thus for simulations in canonical NVT ensemble).
+"""
+from __future__ import division
+
+import os.path as path
+import numpy as np
+
+from MDAnalysis.analysis.base import AnalysisBase
+
+class LinearDensity(AnalysisBase):
+    """Linear density profile
+    LinearDensity(selection, grouping='atoms', binsize=0.25)
+
+    Parameters
+    ----------
+    selection : `AtomGroup` object
+      Any atomgroup
+
+    Keywords
+    --------
+    grouping : str {'atoms', 'residues', 'segments', 'fragments'}
+          Density profiles will be computed on the center of geometry
+          of a selected group of atoms ['atoms']
+
+    binsize : float
+          Bin width in Angstrom used to build linear density
+          histograms. Defines the resolution of the resulting density 
+          profile (smaller --> higher resolution) [0.25]
+
+    start : int
+          The frame to start at [0]
+    stop : int
+          The frame to end at [-1]
+    step : int
+          The step size through the trajectory in frames [0]
+
+    Example
+    -------
+    First create a LinearDensity object by supplying a selection, 
+    then use the `run` method:
+      ldens = LinearDensity(selection)
+      ldens.run()
+
+    Density profiles can be written to file through the `save` method:
+      ldens.save(description='mydensprof', form='txt')
+    which will output the density profiles in a file named
+    <trajectory_filename>.mydensprof_<grouping>.ldens
+    Results can be saved in npz format by specifying `form='npz'`
+
+    .. versionadded:: 0.14.0
+    """
+    def __init__(self, selection, grouping='atoms', binsize=0.25,
+                 start=None, stop=None, step=None):
+        self._ags = [selection] # allows use of run(parallel=True)
+        self._universe = selection.universe
+        self._setup_frames(self._universe.trajectory, start, stop, step)
+
+        self.binsize = binsize
+
+        # group of atoms on which to compute the COM (same as used in
+        # AtomGroup.wrap())
+        self.grouping = grouping
+
+        # Dictionary containing results
+        self.results = {'x': {'dim': 0}, 'y': {'dim': 1}, 'z': {'dim': 2}}
+        # Box sides
+        self.dimensions = self._universe.dimensions[:3]
+        self.volume = np.prod(self.dimensions)
+        bins = (self.dimensions // self.binsize).astype(int) # number of bins
+
+        # Here we choose a number of bins of the largest cell side so that
+        # x, y and z values can use the same "coord" column in the output file
+        self.nbins = bins.max()
+        slices_vol = self.volume/bins
+
+        self.keys = ['pos', 'pos_std', 'char', 'char_std']
+
+        # Initialize results array with zeros
+        for dim in self.results:
+            idx = self.results[dim]['dim']
+            self.results[dim].update({'slice volume': slices_vol[idx]})
+            for key in self.keys:
+                self.results[dim].update({key: np.zeros(self.nbins)})
+
+        # Variables later defined in _prepare() method
+        self.masses = None
+        self.charges = None
+        self.totalmass = None
+
+    def _prepare(self):
+        # group must be a local variable, otherwise there will be
+        # issues with parallelization
+        group = getattr(self._ags[0], self.grouping)
+
+        # Get masses and charges for the selection
+        try: # in case it's not an atom
+            self.masses = np.array([elem.total_mass() for elem in group])
+            self.charges = np.array([elem.total_charge() for elem in group])
+        except AttributeError: # much much faster for atoms
+            self.masses = self._ags[0].masses
+            self.charges = self._ags[0].charges
+
+        self.totalmass = np.sum(self.masses)
+
+    def _single_frame(self):
+        self.group = getattr(self._ags[0], self.grouping)
+        self._ags[0].wrap(compound=self.grouping)
+
+        # Find position of atom/group of atoms
+        if self.grouping == 'atoms':
+            positions = self._ags[0].positions # faster for atoms
+        else:
+            positions = np.array([elem.centroid() for elem in self.group]) # COM for res/frag/etc
+
+        for dim in ['x', 'y', 'z']:
+            idx = self.results[dim]['dim']
+
+            key = 'pos'
+            key_std = 'pos_std'
+            # histogram for positions weighted on masses
+            hist, _ = np.histogram(positions[:, idx], weights=self.masses,
+                                   bins=self.nbins, range=(0.0, max(self.dimensions)))
+
+            self.results[dim][key] += hist
+            self.results[dim][key_std] += np.square(hist)
+
+            key = 'char'
+            key_std = 'char_std'
+            # histogram for positions weighted on charges
+            hist, _ = np.histogram(positions[:, idx], weights=self.charges,
+                                   bins=self.nbins, range=(0.0, max(self.dimensions)))
+
+            self.results[dim][key] += hist
+            self.results[dim][key_std] += np.square(hist)
+
+    def _conclude(self):
+        k = 6.022e-1 # divide by avodagro and convert from A3 to cm3
+
+        # Average results over the  number of configurations
+        for dim in ['x', 'y', 'z']:
+            for key in ['pos', 'pos_std', 'char', 'char_std']:
+                self.results[dim][key] /= self.nframes
+            # Compute standard deviation for the error
+            self.results[dim]['pos_std'] = np.sqrt(self.results[dim]['pos_std']
+                                                   - np.square(self.results[dim]['pos']))
+            self.results[dim]['char_std'] = np.sqrt(self.results[dim]['char_std']
+                                                    - np.square(self.results[dim]['char']))
+
+        for dim in ['x', 'y', 'z']:
+            self.results[dim]['pos'] /= (self.results[dim]['slice volume']*k)
+            self.results[dim]['char'] /= (self.results[dim]['slice volume']*k)
+            self.results[dim]['pos_std'] /= (self.results[dim]['slice volume']*k)
+            self.results[dim]['char_std'] /= (self.results[dim]['slice volume']*k)
+
+    def save(self, description='', form='txt'):
+        """ Save density profile to file
+        Allows to save the density profile to either a ASCII txt file or a
+        binary numpy npz file. Output file has extension 'ldens' and begins
+        with the name of the trajectory file.
+
+        Keywords
+        --------
+        description : str
+          An arbitrary description added to the output filename. Can be useful
+        form : str {'txt', 'npz'}
+          Format of the output. 'txt' will generate an ASCII text file while 'npz'
+          will produce a numpy binary file.
+
+        Example
+        -------
+        After initializing and running a `LinearDensity` object, results can be written
+        to file as follows:
+          ldens.save(description='mydensprof', form='txt')
+        which will output the linear density profiles in a file named
+        <trajectory_filename>.mydensprof_<grouping>.ldens
+        """
+        # Take root of trajectory filename for output file naming
+        trajname = path.splitext(path.basename(self._universe.trajectory.filename))[0]
+        # additional string for naming the output file
+        description = description + "_" + str(self.grouping)
+        filename = trajname+"." + description + ".ldens"
+
+        if form is 'txt':
+            self._savetxt(filename)
+        elif form is 'npz':
+            self._savez(filename)
+        else:
+            raise ValueError('form argument must be either txt or npz')
+
+    def _savetxt(self, filename):
+        bins = np.linspace(0.0, max(self.dimensions), num=self.nbins)
+
+        # Create list of results which will be output
+        output = [bins]
+
+        for dim in ['x', 'y', 'z']:
+            output.append(self.results[dim]['pos'])
+            output.append(self.results[dim]['pos_std'])
+
+        for dim in ['x', 'y', 'z']:
+            output.append(self.results[dim]['char'])
+            output.append(self.results[dim]['char_std'])
+
+        density = self.totalmass/self.volume
+        header = "1 coord [Ang] 2-7 mass density (x,sx,y,sz,z,sz) [g/cm^3]" + \
+            "8-13 charge density (x,sx,y,sz,z,sz) [e/A^3]\n Average density: "\
+            + str(density) + " g/cm3"
+        np.savetxt(filename, np.column_stack(output), fmt='%10.5f', header=header)
+
+    def _savez(self, filename):
+        bins = np.linspace(0.0, max(self.dimensions), num=self.nbins)
+        dictionary = {'bins': bins}
+
+        for dim in self.results:
+            self.results[dim].pop('dim')
+            self.results[dim].pop('slice volume')
+            for key in self.results[dim]:
+                dictionary[dim+"_"+key] = self.results[dim][key]
+
+        np.savez(filename, **dictionary)
+
+    def _add_other_results(self, other):
+        # For parallel analysis
+        results = self.results
+        for dim in ['x', 'y', 'z']:
+            key = 'pos'
+            key_std = 'pos_std'
+            results[dim][key] += other[dim][key]
+            results[dim][key_std] += other[dim][key_std]
+
+            key = 'char'
+            key_std = 'char_std'
+            results[dim][key] += other[dim][key]
+            results[dim][key_std] += other[dim][key_std]

package/doc/sphinx/source/conf.py

@@ -85,7 +85,7 @@ def __getattr__(cls, name):
     Feltz, Philip Fowler, Joseph Goose, Richard J. Gowers, Lukas
     Grossar, Benjamin Hall, Joe Jordan, Max Linke, Jinju Lu, Robert
     McGibbon, Alex Nesterenko, Manuel Nuno Melo, Hai Nguyen,
-    Caio S. Souza, Danny Parton, Joshua L. Phillips, Tyler Reddy,
+    Caio S. Souza, Mattia F. Palermo, Danny Parton, Joshua L. Phillips, Tyler Reddy,
     Paul Rigor, Sean L. Seyler, Andy Somogyi, Lukas Stelzl,
     Gorman Stock, Isaac Virshup, Zhuyi Xue, Carlos Yáñez S.,
     and Oliver Beckstein"""

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

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

package/doc/sphinx/source/documentation_pages/analysis_modules.rst

@@ -106,6 +106,7 @@ Volumetric analysis
    :maxdepth: 1
 
    analysis/density
+   analysis/lineardensity
    analysis/waterdynamics
 
 

testsuite/MDAnalysisTests/analysis/test_lineardensity.py

@@ -0,0 +1,42 @@
+# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*-
+# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 fileencoding=utf-8
+#
+# MDAnalysis --- http://www.MDAnalysis.org
+# Copyright (c) 2006-2015 Naveen Michaud-Agrawal, Elizabeth J. Denning, Oliver Beckstein
+# and contributors (see AUTHORS for the full list)
+#
+# Released under the GNU Public Licence, v2 or any higher version
+#
+# Please cite your use of MDAnalysis in published work:
+#
+# 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 MDAnalysis as mda
+import numpy as np
+
+from MDAnalysisTests.datafiles import waterPSF, waterDCD
+from MDAnalysis.analysis.lineardensity import LinearDensity
+from numpy.testing import TestCase, assert_allclose
+
+class TestLinearDensity(TestCase):
+    def setUp(self):
+        self.universe = mda.Universe(waterPSF, waterDCD)
+        self.sel_string = 'all'
+        self.selection = self.universe.select_atoms(self.sel_string)
+
+        self.xpos = np.array([ 0., 0., 0., 0.0072334, 0.00473299, 0.,
+                               0., 0., 0., 0. ])
+
+    def test_serial(self):
+        ld = LinearDensity(self.selection, binsize=5)
+        ld.run()
+        assert_allclose(self.xpos, ld.results['x']['pos'], rtol=1e-6, atol=0)
+
+#    def test_parallel(self):
+#        ld = LinearDensity(self.universe, self.selection, binsize=5)
+#        ld.run(parallel=True)
+#        assert_equal(self.xpos, ld.results['x']['pos'])
+
+