From b6dfa8809d1a006a26d290c9870b62f22b74dd60 Mon Sep 17 00:00:00 2001 From: mattihappy Date: Wed, 27 Jan 2016 13:40:53 +0100 Subject: [PATCH 01/10] Added lineardensity module --- package/MDAnalysis/analysis/lineardensity.py | 146 +++++++++++++++++++ 1 file changed, 146 insertions(+) create mode 100644 package/MDAnalysis/analysis/lineardensity.py diff --git a/package/MDAnalysis/analysis/lineardensity.py b/package/MDAnalysis/analysis/lineardensity.py new file mode 100644 index 00000000000..c100685bdbb --- /dev/null +++ b/package/MDAnalysis/analysis/lineardensity.py @@ -0,0 +1,146 @@ +# -*- coding: utf-8 -*- +# pylint: disable=E1101 +from __future__ import division +import os.path as path +from MDAnalysis.analysis.base import AnalysisBase +import numpy as np + +class LinearDensity(AnalysisBase): + def __init__(self, universe, selection, description="", grouping='atoms', binsize=0.25, + start=None, stop=None, step=None): + self._ags = [selection] # allows use of run(parallel=True) + self._setup_frames(universe, start, stop, step) + + self.binsize = binsize + + # group of atoms on which to compute the COM (same as used in + # AtomGroup.wrap()) + self.grouping = grouping + + # Take root of trajectory filename for output file naming + self.trajname = path.splitext(path.basename(universe.trajectory.dcdfilename))[0] + # additional string for naming the output file + self.description = description+"_"+str(grouping) + + # Dictionary containing results + self.results = {'x': {'dim': 0}, 'y': {'dim': 1}, 'z': {'dim': 2}} + # Box sides + self.dimensions = [side for side in self._universe.dimensions[0:3]] + self.volume = np.prod(self.dimensions) + bins = [int(side//self.binsize) for side in self.dimensions] # 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 = max(bins) + slices_vol = [self.volume/n for n in 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 = [elem.total_mass() for elem in group] + self.charges = [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) + + #@profile + def _single_frame(self, timestep): + 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] = np.sum([self.results[dim][key], hist], axis=0) + self.results[dim][key_std] = np.sum([self.results[dim][key_std], + np.square(hist)], axis=0) + + 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] = np.sum([self.results[dim][key], hist], axis=0) + self.results[dim][key_std] = np.sum([self.results[dim][key_std], np.square(hist)], + axis=0) + + def _conclude(self): + k = 6.022e-1 # divide by avodagro and convert from A3 to cm3 + bins = np.linspace(0.0, max(self.dimensions), num=self.nbins) + + # 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 + # Computed 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'])) + + + # Create list of results which will be output + output = [bins] + + for dim in ['x', 'y', 'z']: + output.append(self.results[dim]['pos']/(self.results[dim]['slice volume']*k)) + output.append(self.results[dim]['pos_std']/(self.results[dim]['slice volume']*k)) + + for dim in ['x', 'y', 'z']: + output.append(self.results[dim]['char']/(self.results[dim]['slice volume']*k)) + output.append(self.results[dim]['char_std']/(self.results[dim]['slice volume']*k)) + + # Define filename and write to output + filename = self.trajname+"."+self.description+".ldens" + 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), header=header) + + 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] = np.sum([results[dim][key], other[dim][key]], axis=0) + results[dim][key_std] = np.sum([results[dim][key_std], other[dim][key_std]], axis=0) + + key = 'char' + key_std = 'char_std' + results[dim][key] = np.sum([results[dim][key], other[dim][key]], axis=0) + results[dim][key_std] = np.sum([results[dim][key_std], other[dim][key_std]], axis=0) From eadb51af0d4b7bc60930b17143f4335e3e8ea144 Mon Sep 17 00:00:00 2001 From: mattihappy Date: Wed, 27 Jan 2016 19:26:58 +0100 Subject: [PATCH 02/10] Fixed: np arrays instead of lists for charges and masses, correct import order, added save method, use of in-place sum of np arrays --- package/MDAnalysis/analysis/lineardensity.py | 90 +++++++++++++------- 1 file changed, 58 insertions(+), 32 deletions(-) diff --git a/package/MDAnalysis/analysis/lineardensity.py b/package/MDAnalysis/analysis/lineardensity.py index c100685bdbb..066b13364e6 100644 --- a/package/MDAnalysis/analysis/lineardensity.py +++ b/package/MDAnalysis/analysis/lineardensity.py @@ -1,27 +1,24 @@ # -*- coding: utf-8 -*- # pylint: disable=E1101 from __future__ import division + import os.path as path -from MDAnalysis.analysis.base import AnalysisBase import numpy as np +from MDAnalysis.analysis.base import AnalysisBase + class LinearDensity(AnalysisBase): - def __init__(self, universe, selection, description="", grouping='atoms', binsize=0.25, + def __init__(self, universe, selection, grouping='atoms', binsize=0.25, start=None, stop=None, step=None): self._ags = [selection] # allows use of run(parallel=True) self._setup_frames(universe, start, stop, step) self.binsize = binsize - + # group of atoms on which to compute the COM (same as used in # AtomGroup.wrap()) self.grouping = grouping - # Take root of trajectory filename for output file naming - self.trajname = path.splitext(path.basename(universe.trajectory.dcdfilename))[0] - # additional string for naming the output file - self.description = description+"_"+str(grouping) - # Dictionary containing results self.results = {'x': {'dim': 0}, 'y': {'dim': 1}, 'z': {'dim': 2}} # Box sides @@ -55,11 +52,11 @@ def _prepare(self): # Get masses and charges for the selection try: # in case it's not an atom - self.masses = [elem.total_mass() for elem in group] - self.charges = [elem.total_charge() for elem in group] + 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.masses = np.array(self._ags[0].masses) + self.charges = np.array(self._ags[0].charges) self.totalmass = np.sum(self.masses) @@ -83,9 +80,8 @@ def _single_frame(self, timestep): hist, _ = np.histogram(positions[:, idx], weights=self.masses, bins=self.nbins, range=(0.0, max(self.dimensions))) - self.results[dim][key] = np.sum([self.results[dim][key], hist], axis=0) - self.results[dim][key_std] = np.sum([self.results[dim][key_std], - np.square(hist)], axis=0) + self.results[dim][key] += hist + self.results[dim][key_std] += np.square(hist) key = 'char' key_std = 'char_std' @@ -93,43 +89,73 @@ def _single_frame(self, timestep): hist, _ = np.histogram(positions[:, idx], weights=self.charges, bins=self.nbins, range=(0.0, max(self.dimensions))) - self.results[dim][key] = np.sum([self.results[dim][key], hist], axis=0) - self.results[dim][key_std] = np.sum([self.results[dim][key_std], np.square(hist)], - axis=0) + 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 - bins = np.linspace(0.0, max(self.dimensions), num=self.nbins) # 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 - # Computed standard deviation for the error + # 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'): + # Take root of trajectory filename for output file naming + trajname = path.splitext(path.basename(self._universe.trajectory.dcdfilename))[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']/(self.results[dim]['slice volume']*k)) - output.append(self.results[dim]['pos_std']/(self.results[dim]['slice volume']*k)) + 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']/(self.results[dim]['slice volume']*k)) - output.append(self.results[dim]['char_std']/(self.results[dim]['slice volume']*k)) + output.append(self.results[dim]['char']) + output.append(self.results[dim]['char_std']) - # Define filename and write to output - filename = self.trajname+"."+self.description+".ldens" 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), header=header) + + 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 @@ -137,10 +163,10 @@ def _add_other_results(self, other): for dim in ['x', 'y', 'z']: key = 'pos' key_std = 'pos_std' - results[dim][key] = np.sum([results[dim][key], other[dim][key]], axis=0) - results[dim][key_std] = np.sum([results[dim][key_std], other[dim][key_std]], axis=0) + results[dim][key] += other[dim][key] + results[dim][key_std] += other[dim][key_std] key = 'char' key_std = 'char_std' - results[dim][key] = np.sum([results[dim][key], other[dim][key]], axis=0) - results[dim][key_std] = np.sum([results[dim][key_std], other[dim][key_std]], axis=0) + results[dim][key] += other[dim][key] + results[dim][key_std] += other[dim][key_std] From f38bc451a3c3c23575d14ba03aa99e1c364ac208 Mon Sep 17 00:00:00 2001 From: mattihappy Date: Sun, 31 Jan 2016 02:50:22 +0100 Subject: [PATCH 03/10] Fixed improper use of list comprehension; Added docstrings --- package/MDAnalysis/analysis/lineardensity.py | 91 ++++++++++++++++++-- 1 file changed, 83 insertions(+), 8 deletions(-) diff --git a/package/MDAnalysis/analysis/lineardensity.py b/package/MDAnalysis/analysis/lineardensity.py index 066b13364e6..6387a61f68d 100644 --- a/package/MDAnalysis/analysis/lineardensity.py +++ b/package/MDAnalysis/analysis/lineardensity.py @@ -1,5 +1,14 @@ # -*- coding: utf-8 -*- # pylint: disable=E1101 + +""" +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 @@ -8,6 +17,52 @@ from MDAnalysis.analysis.base import AnalysisBase class LinearDensity(AnalysisBase): + """Linear density profile + LinearDensity(universe, selection, grouping='atoms', binsize=0.25) + + Parameters + ---------- + universe : `Universe` object + universe whose density profiles will be computed + selection : `AtomGroup` object + Any atomgroup belonging to universe + + 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 universe and + a selection, then use the `run` method + ldens = LinearDensity(universe, selection) + ldens.run() + Analysis can be run in parallel on multicore machines by supplying + the argument `parallel=True` to the `run` method: + ldens.run(parallel=True, nthreads=4) + + Density proiles are output to file through the `save` method: + ldens.save(description='mydensprof', form='txt') + which will write the density profiles in a file named + .mydensprof_.ldens + Results can be output in npz format by specifying `form='npz'` + + .. versionadded:: 0.14.0 + """ def __init__(self, universe, selection, grouping='atoms', binsize=0.25, start=None, stop=None, step=None): self._ags = [selection] # allows use of run(parallel=True) @@ -22,14 +77,14 @@ def __init__(self, universe, selection, grouping='atoms', binsize=0.25, # Dictionary containing results self.results = {'x': {'dim': 0}, 'y': {'dim': 1}, 'z': {'dim': 2}} # Box sides - self.dimensions = [side for side in self._universe.dimensions[0:3]] + self.dimensions = self._universe.dimensions[:3] self.volume = np.prod(self.dimensions) - bins = [int(side//self.binsize) for side in self.dimensions] # number of bins + 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 = max(bins) - slices_vol = [self.volume/n for n in bins] + self.nbins = bins.max() + slices_vol = self.volume/bins self.keys = ['pos', 'pos_std', 'char', 'char_std'] @@ -55,13 +110,12 @@ def _prepare(self): 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 = np.array(self._ags[0].masses) - self.charges = np.array(self._ags[0].charges) + self.masses = self._ags[0].masses + self.charges = self._ags[0].charges self.totalmass = np.sum(self.masses) - #@profile - def _single_frame(self, timestep): + def _single_frame(self, _): self.group = getattr(self._ags[0], self.grouping) self._ags[0].wrap(compound=self.grouping) @@ -112,6 +166,27 @@ def _conclude(self): 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 + .mydensprof_.ldens + """ # Take root of trajectory filename for output file naming trajname = path.splitext(path.basename(self._universe.trajectory.dcdfilename))[0] # additional string for naming the output file From a6771157b071c50f2275c3be8650fba915d8f549 Mon Sep 17 00:00:00 2001 From: mattihappy Date: Wed, 3 Feb 2016 00:23:35 +0100 Subject: [PATCH 04/10] Added standard header --- package/MDAnalysis/analysis/lineardensity.py | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/package/MDAnalysis/analysis/lineardensity.py b/package/MDAnalysis/analysis/lineardensity.py index 6387a61f68d..f68dd77bc9d 100644 --- a/package/MDAnalysis/analysis/lineardensity.py +++ b/package/MDAnalysis/analysis/lineardensity.py @@ -1,5 +1,18 @@ -# -*- coding: utf-8 -*- -# pylint: disable=E1101 +# -*- 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` From 5605394a5f70afd8665c8d2c6bad38979369be7c Mon Sep 17 00:00:00 2001 From: mattihappy Date: Wed, 3 Feb 2016 00:28:24 +0100 Subject: [PATCH 05/10] Created stub rst for the analysis module and added link under Volumetric Analysis --- .../source/documentation_pages/analysis/lineardensity.rst | 2 ++ .../doc/sphinx/source/documentation_pages/analysis_modules.rst | 1 + 2 files changed, 3 insertions(+) create mode 100644 package/doc/sphinx/source/documentation_pages/analysis/lineardensity.rst diff --git a/package/doc/sphinx/source/documentation_pages/analysis/lineardensity.rst b/package/doc/sphinx/source/documentation_pages/analysis/lineardensity.rst new file mode 100644 index 00000000000..d46d87dff90 --- /dev/null +++ b/package/doc/sphinx/source/documentation_pages/analysis/lineardensity.rst @@ -0,0 +1,2 @@ +.. automodule:: MDAnalysis.analysis.lineardensity + diff --git a/package/doc/sphinx/source/documentation_pages/analysis_modules.rst b/package/doc/sphinx/source/documentation_pages/analysis_modules.rst index 7cd3965c4fb..df3f66c98ad 100644 --- a/package/doc/sphinx/source/documentation_pages/analysis_modules.rst +++ b/package/doc/sphinx/source/documentation_pages/analysis_modules.rst @@ -106,6 +106,7 @@ Volumetric analysis :maxdepth: 1 analysis/density + analysis/lineardensity analysis/waterdynamics From a3a9550f77a726e856ccb9a33a1f4bcbf46ba0e1 Mon Sep 17 00:00:00 2001 From: mattihappy Date: Wed, 3 Feb 2016 00:40:23 +0100 Subject: [PATCH 06/10] Added myself to authors --- package/AUTHORS | 1 + package/doc/sphinx/source/conf.py | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/package/AUTHORS b/package/AUTHORS index a1b550c40ca..d042c7d7cec 100644 --- a/package/AUTHORS +++ b/package/AUTHORS @@ -66,6 +66,7 @@ Chronological list of authors - Hai Nguyen 2016 - Balasubramanian + - Mattia F. Palermo External code ------------- diff --git a/package/doc/sphinx/source/conf.py b/package/doc/sphinx/source/conf.py index a63641c803b..bfb6a470edc 100644 --- a/package/doc/sphinx/source/conf.py +++ b/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""" From 38d09a272d43465b1ea1b2e0cef695174dd6d86e Mon Sep 17 00:00:00 2001 From: mattihappy Date: Wed, 3 Feb 2016 00:49:03 +0100 Subject: [PATCH 07/10] Modified changelog --- package/CHANGELOG | 3 +++ 1 file changed, 3 insertions(+) diff --git a/package/CHANGELOG b/package/CHANGELOG index 485294be9a8..6797e5167c4 100644 --- a/package/CHANGELOG +++ b/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 From 90c574cd6b5169a1d5b72c3dbeeda80533a9115d Mon Sep 17 00:00:00 2001 From: mattihappy Date: Wed, 3 Feb 2016 22:21:01 +0100 Subject: [PATCH 08/10] Lineardensity does not take universe as argument anymore and is compatible with AnalysisBase of MDAnalysis 0.13.0 --- package/MDAnalysis/analysis/lineardensity.py | 30 +++++++++----------- 1 file changed, 13 insertions(+), 17 deletions(-) diff --git a/package/MDAnalysis/analysis/lineardensity.py b/package/MDAnalysis/analysis/lineardensity.py index f68dd77bc9d..700c98b2b29 100644 --- a/package/MDAnalysis/analysis/lineardensity.py +++ b/package/MDAnalysis/analysis/lineardensity.py @@ -31,14 +31,12 @@ class LinearDensity(AnalysisBase): """Linear density profile - LinearDensity(universe, selection, grouping='atoms', binsize=0.25) + LinearDensity(selection, grouping='atoms', binsize=0.25) Parameters ---------- - universe : `Universe` object - universe whose density profiles will be computed selection : `AtomGroup` object - Any atomgroup belonging to universe + Any atomgroup Keywords -------- @@ -60,26 +58,24 @@ class LinearDensity(AnalysisBase): Example ------- - First create a LinearDensity object by supplying a universe and - a selection, then use the `run` method - ldens = LinearDensity(universe, selection) + First create a LinearDensity object by supplying a selection, + then use the `run` method: + ldens = LinearDensity(selection) ldens.run() - Analysis can be run in parallel on multicore machines by supplying - the argument `parallel=True` to the `run` method: - ldens.run(parallel=True, nthreads=4) - - Density proiles are output to file through the `save` method: + + Density profiles can be written to file through the `save` method: ldens.save(description='mydensprof', form='txt') - which will write the density profiles in a file named + which will output the density profiles in a file named .mydensprof_.ldens - Results can be output in npz format by specifying `form='npz'` + Results can be saved in npz format by specifying `form='npz'` .. versionadded:: 0.14.0 """ - def __init__(self, universe, selection, grouping='atoms', binsize=0.25, + 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._setup_frames(universe, start, stop, step) + self._universe = selection.universe + self._setup_frames(self._universe.trajectory, start, stop, step) self.binsize = binsize @@ -128,7 +124,7 @@ def _prepare(self): self.totalmass = np.sum(self.masses) - def _single_frame(self, _): + def _single_frame(self): self.group = getattr(self._ags[0], self.grouping) self._ags[0].wrap(compound=self.grouping) From fd44b72854108e76a053177a664e6ad0f7e52232 Mon Sep 17 00:00:00 2001 From: mattihappy Date: Wed, 3 Feb 2016 22:21:32 +0100 Subject: [PATCH 09/10] Test added --- .../analysis/test_lineardensity.py | 44 +++++++++++++++++++ 1 file changed, 44 insertions(+) create mode 100644 testsuite/MDAnalysisTests/analysis/test_lineardensity.py diff --git a/testsuite/MDAnalysisTests/analysis/test_lineardensity.py b/testsuite/MDAnalysisTests/analysis/test_lineardensity.py new file mode 100644 index 00000000000..5f964d6731c --- /dev/null +++ b/testsuite/MDAnalysisTests/analysis/test_lineardensity.py @@ -0,0 +1,44 @@ +# -*- 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() + print self.xpos + print ld.results['x']['pos'] + 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']) + + From 0f0aefac3a9987cfc20dc2454a28276f05c2ddaa Mon Sep 17 00:00:00 2001 From: mattihappy Date: Sat, 6 Feb 2016 22:22:37 +0100 Subject: [PATCH 10/10] Trajectory name taken from variable filename and not dcdfilename, removed print statements from test --- package/MDAnalysis/analysis/lineardensity.py | 2 +- testsuite/MDAnalysisTests/analysis/test_lineardensity.py | 2 -- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/package/MDAnalysis/analysis/lineardensity.py b/package/MDAnalysis/analysis/lineardensity.py index 700c98b2b29..511aa13838e 100644 --- a/package/MDAnalysis/analysis/lineardensity.py +++ b/package/MDAnalysis/analysis/lineardensity.py @@ -197,7 +197,7 @@ def save(self, description='', form='txt'): .mydensprof_.ldens """ # Take root of trajectory filename for output file naming - trajname = path.splitext(path.basename(self._universe.trajectory.dcdfilename))[0] + 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" diff --git a/testsuite/MDAnalysisTests/analysis/test_lineardensity.py b/testsuite/MDAnalysisTests/analysis/test_lineardensity.py index 5f964d6731c..6c861518d95 100644 --- a/testsuite/MDAnalysisTests/analysis/test_lineardensity.py +++ b/testsuite/MDAnalysisTests/analysis/test_lineardensity.py @@ -32,8 +32,6 @@ def setUp(self): def test_serial(self): ld = LinearDensity(self.selection, binsize=5) ld.run() - print self.xpos - print ld.results['x']['pos'] assert_allclose(self.xpos, ld.results['x']['pos'], rtol=1e-6, atol=0) # def test_parallel(self):