3.8.1. Generating densities from trajectories — MDAnalysis.analysis.density
¶
Author: | Oliver Beckstein |
---|---|
Year: | 2011 |
Copyright: | GNU Public License v3 |
The module provides classes and functions to generate and represent volumetric data, in particular densities.
3.8.1.1. Generating a density from a MD trajectory¶
An input trajectory is required that
- Has been centered on the protein of interest.
- Has all molecules made whole that have been broken across periodic boundaries.
- Has the solvent molecules remapped so that they are closest to the solute (this is important when using funky unit cells such as a dodecahedron or a truncated octahedron).
To generate the density of water molecules around a protein:
from MDAnalysis.analysis.density import density_from_Universe
u = Universe(PSF,DCD)
D = density_from_Universe(u, delta=1.0, atomselection="name OH2")
D.convert_density('TIP3P')
D.export("water.dx")
The positions of all water oxygens are histogrammed on a grid with spacing
delta = 1 A. Initially the density is measured in 1/A**3. With the
Density.convert_density()
method, the units of measurement are
changed. In the example we are now measuring the density relative to the
literature value of the TIP3P water model at ambient conditions (see the values
in MDAnalysis.units.water
for details). Finally, the density is
writte as an OpenDX compatible file that can be read in VMD or PyMOL.
See Density
for details. In particular, the density is stored
as a NumPy array in Density.grid
, which can be processed in
any manner.
3.8.1.2. Classes and Functions¶
-
class
MDAnalysis.analysis.density.
Density
(*args, **kwargs)[source]¶ Bases:
gridData.core.Grid
Class representing a density on a regular cartesian grid.
The data (
Density.grid
) can be manipulated as a standard numpy array. Changes can be saved to a file using theDensity.save()
method. The grid can be restored using theDensity.load()
method or by supplying the filename to the constructor.The attribute
Density.metadata
holds a user-defined dictionary that can be used to annotate the data. It is also saved withDensity.save()
.The
Density.export()
method always exports a 3D object (written in such a way to be readable in VMD and PyMOL), the rest should work for an array of any dimension.If the input histogram consists of counts per cell then the
Density.make_density()
method converts the grid to a physical density. For a probability density, divide it byDensity.grid.sum()
or usenormed=True
right away inhistogramdd()
.The user should set the parameters keyword (see docs for the constructor); in particular, if the data are already a density, one must set isDensity ==
True
because there is no reliable way to detect if data represent counts or a density. As a special convenience, if data are read from a file and the user has not set isDensity then it is assumed that the data are in fact a density.Typical use:
From a histogram (i.e. counts on a grid):
h,edges = numpy.histogramdd(...) D = Density(h, edges, parameters={'isDensity': False}, units={'length': 'A'}) D.make_density()
From a saved density file (e.g. in OpenDX format), where the lengths are in Angstrom and the density in 1/A**3:
D = Density("density.dx")
From a saved density file (e.g. in OpenDX format), where the lengths are in Angstrom and the density is measured relative to the density of water at ambient conditions:
D = Density("density.dx", units={'density': 'water'})
From a saved histogram (less common, but in order to demonstrate the parameters keyword) where the lengths are in nm:
D = Density("counts.dx", parameters={'isDensity': False}, units={'length': 'nm'}) D.make_density() D.convert_length('Angstrom^{-3}') D.convert_density('water')
After the final step,
D
will contain a density on a grid measured in Angstrom, with the density values itself measured relative to the density of water.
Density
objects can be algebraically manipulated (added, subtracted, multiplied, ...) but there are no sanity checks in place to make sure that units, metadata, etc are compatible!Note
It is suggested to construct the Grid object from a histogram, to supply the appropriate length unit, and to use
Density.make_density()
to obtain a density. This ensures that the length- and the density unit correspond to each other.See also
Grid
which is the base class ofDensity
. (Grid
has been imported fromgridData.Grid
which is part of GridDataFormats).Create a
Density
from data.Parameters: - *grid* – histogram or density, typically a
numpy.ndarray
- *edges* – list of arrays, the lower and upper bin edges along the axes
- *parameters* –
dictionary of class parameters; saved with
Density.save()
. The following keys are meaningful to the class. Meaning of the values are listed:isDensityFalse
: grid is a histogram with counts [default]True
: a density
Applying
Density.make_density`()
sets it toTrue
. - *units* –
A dict with the keys
- length: physical unit of grid edges (Angstrom or nm) [Angstrom]
- density: unit of the density if
isDensity == True
orNone
otherwise; the default is “Angstrom^{-3}” for densities (meaning A^-3).
(Actually, the default unit is the value of
MDAnalysis.core.flags['length_unit']
; in most cases this is “Angstrom”.) - *metadata* – a user defined dictionary of arbitrary values associated with the
density; the class does not touch
Density.metadata
but stores it withDensity.save()
-
__delattr__
¶ x.__delattr__(‘name’) <==> del x.name
-
__format__
()¶ default object formatter
-
__getattribute__
¶ x.__getattribute__(‘name’) <==> x.name
-
__hash__
¶
-
__reduce__
()¶ helper for pickle
-
__reduce_ex__
()¶ helper for pickle
-
__setattr__
¶ x.__setattr__(‘name’, value) <==> x.name = value
-
__sizeof__
() → int¶ size of object in memory, in bytes
-
__str__
¶
-
centers
()¶ Returns the coordinates of the centers of all grid cells as an iterator.
-
check_compatible
(other)¶ Check if other can be used in an arithmetic operation.
- other is a scalar
- other is a grid defined on the same edges
Raises: TypeError
if not compatible.
-
convert_density
(unit='Angstrom')[source]¶ Convert the density to the physical units given by unit.
Grid.convert_to(unit)unit can be one of the following:
name description of the unit Angstrom^{-3} particles/A**3 nm^{-3} particles/nm**3 SPC density of SPC water at standard conditions TIP3P ... see MDAnalysis.units.water
TIP4P ... see MDAnalysis.units.water
water density of real water at standard conditions (0.997 g/cm**3) Molar mol/l Note
- This only works if the initial length unit is provided.
- Conversions always go back to unity so there can be rounding and floating point artifacts for multiple conversions.
There may be some undesirable cross-interactions with
convert_length()
...
-
convert_length
(unit='Angstrom')[source]¶ Convert Grid object to the new unit.
Grid.convert_length(<unit>)Keywords: - unit
Angstrom, nm
This changes the edges but will not change the density; it is the user’s responsibility to supply the appropriate unit if the Grid object is constructed from a density. It is suggested to start from a histogram and a length unit and use
make_density()
.
-
export
(filename, file_format=None)¶ export density to file using the given format.
The format can also be deduced from the suffix of the filename though the format keyword takes precedence.
The default format for export() is ‘dx’. Use ‘dx’ for visualization.
Implemented formats:
- dx
OpenDX
- pickle
- pickle (use :meth:
Grid.load` to restore); :meth:`Grid.save` is simpler than ``export(format='python')
.
-
interpolated
¶ B-spline function over the data grid(x,y,z).
interpolated([x1,x2,...],[y1,y2,...],[z1,z2,...]) -> F[x1,y1,z1],F[x2,y2,z2],...The interpolation order is set in
Grid.interpolation_spline_order
.The interpolated function is computed once and is cached for better performance. Whenever
interpolation_spline_order
is modified,Grid.interpolated()
is recomputed.The value for unknown data is set in
Grid.interpolation_cval
(TODO: also recompute when interpolation_cval value is changed.)- Example usage for resampling::
>>> XX,YY,ZZ = numpy.mgrid[40:75:0.5, 96:150:0.5, 20:50:0.5] >>> FF = interpolated(XX,YY,ZZ)
-
load
(filename, file_format=None)¶ Load saved (pickled or dx) grid and edges from <filename>.pickle
Grid.load(<filename>.pickle) Grid.load(<filename>.dx)The load() method calls the class’s constructor method and completely resets all values, based on the loaded data.
-
make_density
()[source]¶ Convert the grid (a histogram, counts in a cell) to a density (counts/volume).
make_density()- This changes the grid irrevocably.
- For a probability density, manually divide by grid.sum().
If this is already a density, then a warning is issued and nothing is done.
-
resample
(edges)¶ Resample data to a new grid with edges edges.
resample(edges) –> Gridor
resample(otherGrid) –> GridThe order of the interpolation is set by
Grid.interpolation_spline_order
.
-
resample_factor
(factor)¶ Resample to a new regular grid with factor*oldN cells along each dimension.
-
save
(filename)¶ Save a grid object to <filename>.pickle
Grid.save(filename)Internally, this calls Grid.export(filename,format=”python”). A grid can be regenerated from the saved data with
g = Grid(filename=<filename>)
-
MDAnalysis.analysis.density.
density_from_Universe
(universe, delta=1.0, atomselection='name OH2', start=None, stop=None, step=None, metadata=None, padding=2.0, cutoff=0, soluteselection=None, use_kdtree=True, update_selection=False, quiet=False, interval=1, **kwargs)[source]¶ Create a density grid from a
MDAnalysis.Universe
object.The trajectory is read, frame by frame, and the atoms selected with atomselection are histogrammed on a grid with spacing delta:
density_from_Universe(universe, delta=1.0, atomselection='name OH2', ...) --> density
Note
By default, the atomselection is static, i.e., atoms are only selected once at the beginning. If you want dynamically changing selections (such as “name OW and around 4.0 (protein and not name H*)”) then set
update_selection=True
. For the special case of calculating a density of the “bulk” solvent away from a solute use the optimized selections with keywords cutoff and soluteselection.Parameters: universe –
MDAnalysis.Universe
object with a trajectoryKeywords: - atomselection
selection string (MDAnalysis syntax) for the species to be analyzed [“name OH2”]
- delta
bin size for the density grid in Angstroem (same in x,y,z) [1.0]
- start, stop, step
Slice the trajectory as
trajectory[start"stop:step]
; default is to read the whole trajectory.- metadata
dictionary of additional data to be saved with the object
- padding
increase histogram dimensions by padding (on top of initial box size) in Angstroem [2.0]
- soluteselection
MDAnalysis selection for the solute, e.g. “protein” [
None
]- cutoff
With cutoff, select “<atomsel> NOT WITHIN <cutoff> OF <soluteselection>” (Special routines that are faster than the standard
AROUND
selection) [0]- update_selection
Should the selection of atoms be updated for every step? [
False
] -True
: atom selection is updated for each frame, can be slow -False
: atoms are only selected at the beginning- quiet
Print status update to the screen for every interval frame? [
False
] -True
: no status updates when a new frame is processed -False
: status update every frame (including number of atomsprocessed, which is interesting with
update_selection=True
)- interval
Show status update every interval frame [1]
- parameters
dict with some special parameters for
Density
(see doc)- kwargs
metadata, parameters are modified and passed on to
Density
Returns: Changed in version 0.13.0: update_selection and quite keywords added
-
MDAnalysis.analysis.density.
density_from_PDB
(pdb, **kwargs)[source]¶ Create a density from a single frame PDB.
Typical use is to make a density from the crystal water molecules. The density is created from isotropic gaussians centered at each selected atoms. If B-factors are present in the file then they are used to calculate the width of the gaussian.
Using the sigma keyword, one can override this choice and prescribe a gaussian width for all atoms (in Angstrom), which is calculated as described for
Bfactor2RMSF()
.Note
The current implementation is painfully slow.
See also
Parameters: *pdb* – PDB file (should have the temperatureFactor set); ANISO records are currently not processed
Keywords: - atomselection
selection string (MDAnalysis syntax) for the species to be analyzed [‘resname HOH and name O’]
- delta
bin size for the density grid in Angstroem (same in x,y,z) [1.0]
- metadata
dictionary of additional data to be saved with the object [
None
]- padding
increase histogram dimensions by padding (on top of initial box size) [1.0]
- sigma
width (in Angstrom) of the gaussians that are used to build up the density; if
None
then uses B-factors from pdb [None
]
Returns: a
Density
object with a density measured relative to the water density at standard conditions
-
MDAnalysis.analysis.density.
Bfactor2RMSF
(B)[source]¶ Atomic root mean square fluctuation (in Angstrom) from the crystallographic B-factor
RMSF and B-factor are related by [Willis1975]
\[B = \frac{8\pi^2}{3} \rm{RMSF}^2\]and this function returns
\[\rm{RMSF} = \sqrt{\frac{3 B}{8\pi^2}}\]References
[Willis1975] BTM Willis and AW Pryor. Thermal vibrations in crystallography. Cambridge Univ. Press, 1975
-
class
MDAnalysis.analysis.density.
BfactorDensityCreator
(pdb, delta=1.0, atomselection='resname HOH and name O', metadata=None, padding=1.0, sigma=None)[source]¶ Create a density grid from a pdb file using MDAnalysis.
dens = BfactorDensityCreator(pdb,...).Density()The main purpose of this function is to convert crystal waters in an X-ray structure into a density so that one can compare the experimental density with the one from molecular dynamics trajectories. Because a pdb is a single snapshot, the density is estimated by placing Gaussians of width sigma at the position of all selected atoms.
Sigma can be fixed or taken from the B-factor field, in which case sigma is taken as sqrt(3.*B/8.)/pi (see
BFactor2RMSF()
).Construct the density from psf and pdb and the atomselection.
- DC = BfactorDensityCreator(pdb, delta=<delta>, atomselection=<MDAnalysis selection>,
- metadata=<dict>, padding=2, sigma=None)
density = DC.Density()
Parameters: - pdb – PDB file or
MDAnalysis.Universe
; a PDB is read with the simpl PDB reader. If the Bio.PDB reader is required, either set the permissive_pdb_reader flag toFalse
inMDAnalysis.core.flags
or supply a Universe that was created with the permissive =False
keyword. - atomselection – selection string (MDAnalysis syntax) for the species to be analyzed
- delta – bin size for the density grid in Angstroem (same in x,y,z) [1.0]
- metadata – dictionary of additional data to be saved with the object
- padding – increase histogram dimensions by padding (on top of initial box size)
- sigma – width (in Angstrom) of the gaussians that are used to build up the density; if None then uses B-factors from pdb
For assigning X-ray waters to MD densities one might have to use a sigma of about 0.5 A to obtain a well-defined and resolved x-ray water density that can be easily matched to a broader density distribution.
-
class
MDAnalysis.analysis.density.
Grid
(grid=None, edges=None, origin=None, delta=None, metadata={}, interpolation_spline_order=3)¶ Class to manage a multidimensional grid object.
The export(format=’dx’) method always exports a 3D object, the rest should work for an array of any dimension.
The grid (Grid.grid) can be manipulated as a standard numpy array.
The attribute Grid.metadata holds a user-defined dictionary that can be used to annotate the data. It is saved with save().
Create a Grid object from data.
- From a numpy.histogramdd()::
- grid,edges = numpy.histogramdd(...) g = Grid(grid,edges=edges)
- From an arbitrary grid::
- g = Grid(grid,origin=origin,delta=delta)
- From a saved file::
- g = Grid(filename)
- or
- g = Grid() g.load(filename)
Parameters: - grid – histogram or density, defined on numpy nD array
- edges – list of arrays, the lower and upper bin edges along the axes (both are output by numpy.histogramdd())
- origin – cartesian coordinates of the center of grid[0,0,...,0]
- delta – Either n x n array containing the cell lengths in each dimension, or n x 1 array for rectangular arrays.
- metadata – a user defined dictionary of arbitrary values associated with the density; the class does not touch metadata[] but stores it with save()
- interpolation_spline_order – order of interpolation function for resampling; cubic splines = 3 [3]
-
__delattr__
¶ x.__delattr__(‘name’) <==> del x.name
-
__format__
()¶ default object formatter
-
__getattribute__
¶ x.__getattribute__(‘name’) <==> x.name
-
__hash__
¶
-
__reduce__
()¶ helper for pickle
-
__reduce_ex__
()¶ helper for pickle
-
__setattr__
¶ x.__setattr__(‘name’, value) <==> x.name = value
-
__sizeof__
() → int¶ size of object in memory, in bytes
-
__str__
¶
-
centers
()¶ Returns the coordinates of the centers of all grid cells as an iterator.
-
check_compatible
(other)¶ Check if other can be used in an arithmetic operation.
- other is a scalar
- other is a grid defined on the same edges
Raises: TypeError
if not compatible.
-
export
(filename, file_format=None)¶ export density to file using the given format.
The format can also be deduced from the suffix of the filename though the format keyword takes precedence.
The default format for export() is ‘dx’. Use ‘dx’ for visualization.
Implemented formats:
- dx
OpenDX
- pickle
- pickle (use :meth:
Grid.load` to restore); :meth:`Grid.save` is simpler than ``export(format='python')
.
-
interpolated
¶ B-spline function over the data grid(x,y,z).
interpolated([x1,x2,...],[y1,y2,...],[z1,z2,...]) -> F[x1,y1,z1],F[x2,y2,z2],...The interpolation order is set in
Grid.interpolation_spline_order
.The interpolated function is computed once and is cached for better performance. Whenever
interpolation_spline_order
is modified,Grid.interpolated()
is recomputed.The value for unknown data is set in
Grid.interpolation_cval
(TODO: also recompute when interpolation_cval value is changed.)- Example usage for resampling::
>>> XX,YY,ZZ = numpy.mgrid[40:75:0.5, 96:150:0.5, 20:50:0.5] >>> FF = interpolated(XX,YY,ZZ)
-
load
(filename, file_format=None)¶ Load saved (pickled or dx) grid and edges from <filename>.pickle
Grid.load(<filename>.pickle) Grid.load(<filename>.dx)The load() method calls the class’s constructor method and completely resets all values, based on the loaded data.
-
resample
(edges)¶ Resample data to a new grid with edges edges.
resample(edges) –> Gridor
resample(otherGrid) –> GridThe order of the interpolation is set by
Grid.interpolation_spline_order
.
-
resample_factor
(factor)¶ Resample to a new regular grid with factor*oldN cells along each dimension.
-
save
(filename)¶ Save a grid object to <filename>.pickle
Grid.save(filename)Internally, this calls Grid.export(filename,format=”python”). A grid can be regenerated from the saved data with
g = Grid(filename=<filename>)