4.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.

4.8.1.1. Generating a density from a MD trajectory

A common use case is to analyze the solvent density around a protein of interest. The density is calculated with density_from_Universe() in the fixed coordinate system of the simulation unit cell. It is therefore necessary to orient and fix the protein with respect to the box coordinate system. In practice this means centering and superimposing the protein, frame by frame, on a reference structure and translating and rotating all other components of the simulation with the protein. In this way, the solvent will appear in the reference frame of the protein.

An input trajectory must

1. have been centered on the protein of interest;

2. have all molecules made whole that have been broken across periodic boundaries 1;

3. have the solvent molecules remapped so that they are closest to the solute (this is important when using triclinic unit cells such as a dodecahedron or a truncated octahedron) 1.

4. have a fixed frame of reference; for instance, by superimposing a protein on a reference structure so that one can study the solvent density around it 2.

To generate the density of water molecules around a protein (assuming that the trajectory is already appropriately treated for periodic boundary artifacts and is suitably superimposed to provide a fixed reference frame) 3

from MDAnalysis.analysis.density import density_from_Universe
u = Universe(TPR, XTC)
D = density_from_Universe(u, delta=1.0, atomselection="name OW")
D.convert_density('TIP4P')
D.export("water.dx", type="double")

The positions of all water oxygens are histogrammed on a grid with spacing delta = 1 Å. Initially the density is measured in \(\text{Å}^{-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 TIP4P water model at ambient conditions (see the values in MDAnalysis.units.water for details). Finally, the density is written as an OpenDX compatible file that can be read in VMD, Chimera, 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.

4.8.1.2. Creating densities

The following functions take trajectory or coordinate data and generate a Density object.

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, verbose=False, interval=1, quiet=None, parameters=None, gridcenter=None, xdim=None, ydim=None, zdim=None)

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.

Parameters

• universe (MDAnalysis.Universe) – MDAnalysis.Universe object with a trajectory

• atomselection (str (optional)) – selection string (MDAnalysis syntax) for the species to be analyzed [“name OH2”]

• delta (float (optional)) – bin size for the density grid in Angstroem (same in x,y,z) [1.0]

• start (int (optional)) –

• stop (int (optional)) –

• step (int (optional)) – Slice the trajectory as trajectory[start:stop:step]; default is to read the whole trajectory.

• metadata (dict. optional) – dict of additional data to be saved with the object; the meta data are passed through as they are.

• padding (float (optional)) – increase histogram dimensions by padding (on top of initial box size) in Angstroem. Padding is ignored when setting a user defined grid. [2.0]

• soluteselection (str (optional)) – MDAnalysis selection for the solute, e.g. “protein” [None]

• cutoff (float (optional)) – With cutoff, select “<atomsel> NOT WITHIN <cutoff> OF <soluteselection>” (Special routines that are faster than the standard AROUND selection); any value that evaluates to False (such as the default 0) disables this special selection.

• update_selection (bool (optional)) –

  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

• verbose (bool (optional)) –

  Print status update to the screen for every interval frame? [True]

  • False: no status updates when a new frame is processed

  • True: status update every frame (including number of atoms processed, which is interesting with update_selection=True)

• interval (int (optional)) – Show status update every interval frame [1]

• parameters (dict (optional)) – dict with some special parameters for Density (see docs)

• gridcenter (numpy ndarray, float32 (optional)) – 3 element numpy array detailing the x, y and z coordinates of the center of a user defined grid box in Angstroem [None]

• xdim (float (optional)) – User defined x dimension box edge in ångström; ignored if gridcenter is None

• ydim (float (optional)) – User defined y dimension box edge in ångström; ignored if gridcenter is None

• zdim (float (optional)) – User defined z dimension box edge in ångström; ignored if gridcenter is None

Returns

A Density instance with the histogrammed data together with associated metadata.

Return type

Density

Notes

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*)”, i.e., the water oxygen atoms that are within 4 Å of the protein heavy atoms) 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 (see Examples below).

Examples

Basic use for creating a water density (just using the water oxygen atoms “OW”):

density = density_from_Universe(universe, delta=1.0, atomselection='name OW')

If you are only interested in water within a certain region, e.g., within a vicinity around a binding site, you can use a selection that updates every step by setting the update_selection keyword argument:

site_density = density_from_Universe(universe, delta=1.0,
                                     atomselection='name OW and around 5 (resid 156 157 305)',
                                     update_selection=True)

A special case for an updating selection is to create the “bulk density”, i.e., the water outside the immediate solvation shell of a protein: Select all water oxygen atoms that are farther away than a given cut-off (say, 4 Å) from the solute (here, heavy atoms of the protein):

bulk = density_from_Universe(universe, delta=1.0, atomselection='name OW',
                             solute="protein and not name H*",
                             cutoff=4)

(Using the special case for the bulk with soluteselection and cutoff improves performance over the simple update_selection approach.)

If you are interested in explicitly setting a grid box of a given edge size and origin, you can use the gridcenter and x/y/zdim arguments. For example to plot the density of waters within 5 Å of a ligand (in this case the ligand has been assigned the residue name “LIG”) in a cubic grid with 20 Å edges which is centered on the centre of mass (COM) of the ligand:

# Create a selection based on the ligand
ligand_selection = universe.select_atoms("resname LIG")

# Extract the COM of the ligand
ligand_COM = ligand_selection.center_of_mass()

# Generate a density of waters on a cubic grid centered on the ligand COM
# In this case, we update the atom selection as shown above.
water_density = density_from_Universe(universe, delta=1.0,
                                      atomselection='name OW around 5 resname LIG',
                                      update_selection=True,
                                      gridcenter=ligand_COM,
                                      xdim=20.0, ydim=20.0, zdim=20.0)

(It should be noted that the `padding` keyword is not used when a user
defined grid is assigned).

Changed in version 0.20.0: ProgressMeter now iterates over the number of frames analysed.

Changed in version 0.19.0: gridcenter, xdim, ydim and zdim keywords added to allow for user defined boxes

Changed in version 0.13.0: update_selection and quiet keywords added

Deprecated since version 0.16: The keyword argument quiet is deprecated in favor of verbose.

MDAnalysis.analysis.density.density_from_PDB(pdb, **kwargs)

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().

Parameters

• pdb (str) – PDB filename (should have the temperatureFactor set); ANISO records are currently not processed

• atomselection (str) – selection string (MDAnalysis syntax) for the species to be analyzed [‘resname HOH and name O’]

• delta (float) – bin size for the density grid in Angstroem (same in x,y,z) [1.0]

• metadata (dict) – dictionary of additional data to be saved with the object [None]

• padding (float) – increase histogram dimensions by padding (on top of initial box size) [1.0]

• sigma (float) – width (in Angstrom) of the gaussians that are used to build up the density; if None then uses B-factors from pdb [None]

Returns

object with a density measured relative to the water density at standard conditions

Return type

Density

Notes

The current implementation is painfully slow.

See also

Bfactor2RMSF()

MDAnalysis.analysis.density.Bfactor2RMSF(B)

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

4.8.1.3. Supporting classes and functions

The main output of the density creation functions is a Density instance, which is derived from a gridData.core.Grid. A Density is essentially, a 3D array with origin and lengths together with associated metadata (which can be used in downstream processing).

class MDAnalysis.analysis.density.Density(*args, **kwargs)

Bases: gridData.core.Grid

Class representing a density on a regular cartesian grid.

Parameters

• grid (array_like) – histogram or density, typically a numpy.ndarray

• edges (list) – list of arrays, the lower and upper bin edges along the axes

• parameters (dict) –

  dictionary of class parameters; saved with Density.save(). The following keys are meaningful to the class. Meaning of the values are listed:

  isDensity

  • False: grid is a histogram with counts [default]

  • True: a density

  Applying Density.make_density`() sets it to True.

• units (dict) –

  A dict with the keys

  • length: physical unit of grid edges (Angstrom or nm) [Angstrom]

  • density: unit of the density if isDensity=True or None otherwise; the default is “Angstrom^{-3}” for densities (meaning \(\text{Å}^{-3}\)).

  (Actually, the default unit is the value of MDAnalysis.core.flags['length_unit']; in most cases this is “Angstrom”.)

• metadata (dict) – a user defined dictionary of arbitrary values associated with the density; the class does not touch Density.metadata but stores it with Density.save()

grid

counts or density

Type

array

edges

The boundaries of each cell in grid along all axes (equivalent to what numpy.histogramdd() returns).

Type

list of 1d-arrays

delta

Cell size in each dimension.

Type

array

origin

Coordinates of the center of the cell at index grid[0, 0, 0, …, 0], which is considered to be the front lower left corner.

Type

array

units

The units for lengths and density; change units with the method convert_length() or convert_density().

Type

dict

Notes

The data (Density.grid) can be manipulated as a standard numpy array. Changes can be saved to a file using the Density.save() method. The grid can be restored using the Density.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 with Density.save().

The Density.export() method always exports a 3D object (written in such a way to be readable in VMD, Chimera, and PyMOL), the rest should work for an array of any dimension. Note that PyMOL only understands DX files with the DX data type “double” in the “array” object (see known issues when writing OpenDX files and issue MDAnalysis/GridDataFormats#35 for details). Using the keyword type="double" for the method Density.export(), the user can ensure that the DX file is written in a format suitable for PyMOL.

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 by Density.grid.sum() or use normed=True right away in histogramdd().

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.

See also

gridData.core.Grid

the base class of Density.

Examples

Typical use:

1. 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()
   

2. 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")
   

3. 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'})
   

4. 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 Ångstrom, 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.

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)

Arguments

grid

histogram or density, defined on numpy nD array, or filename

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]

file_format

file format; only necessary when grid is a filename (see Grid.load()); default is None and the file format is autodetected.

Changed in version 0.5.0: New file_format keyword argument.

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.

1. other is a scalar

2. other is a grid defined on the same edges

Raises

TypeError if not compatible.

convert_density(unit='Angstrom')

Convert the density to the physical units given by unit.

Parameters

unit (str (optional)) –

The target unit that the density should be converted to.

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

Raises

• RuntimeError – If the density does not have a unit associated with it to begin with (i.e., is not a density) then no conversion can take place.

• ValueError – for unknown unit.

Notes

1. This method only works if there is already a length unit associated with the density; otherwise raises RuntimeError

2. Conversions always go back to unity so there can be rounding and floating point artifacts for multiple conversions.

convert_length(unit='Angstrom')

Convert Grid object to the new unit.

Parameters

unit (str (optional)) – unit that the grid should be converted to: one of “Angstrom”, “nm”

Notes

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, type=None, typequote='"')

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 Grid.load() to restore); Grid.save() is simpler than export(format='python').

Parameters

• filename (str) – name of the output file

• file_format ({'dx', 'pickle', None} (optional)) – output file format, the default is “dx”

• type (str (optional)) –

  for DX, set the output DX array type, e.g., “double” or “float”. By default (None), the DX type is determined from the numpy dtype of the array of the grid (and this will typically result in “double”).

  New in version 0.4.0.

• typequote (str (optional)) –

  For DX, set the character used to quote the type string; by default this is a double-quote character, ‘”’. Custom parsers like the one from NAMD-GridForces (backend for MDFF) expect no quotes, and typequote=’’ may be used to appease them.

  New in version 0.5.0.

property interpolated

B-spline function over the data grid(x,y,z).

The interpolated() function allows one to obtain data values for any values of the coordinates:

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

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)

Note

Values are interpolated with a spline function. It is possible that the spline will generate values that would not normally appear in the data. For example, a density is non-negative but a cubic spline interpolation can generate negative values, especially at the boundary between 0 and high values.

property interpolation_spline_order

Order of the B-spline interpolation of the data.

3 = cubic; 4 & 5 are also supported

Only choose values that are acceptable to scipy.ndimage.spline_filter()!

See also

interpolated

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()

Convert the grid (a histogram, counts in a cell) to a density (counts/volume).

This method 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, so calling make_density multiple times does not do any harm.

resample(edges)

Resample data to a new grid with edges edges.

This method creates a new grid with the data from the current grid resampled to a regular grid specified by edges. The order of the interpolation is set by Grid.interpolation_spline_order: change the value before calling resample().

Parameters

edges (tuple of arrays or Grid) – edges of the new grid or a Grid instance that provides Grid.edges

Returns

a new Grid with the data interpolated over the new grid cells

Return type

Grid

Examples

Providing edges (a tuple of three arrays, indicating the boundaries of each grid cell):

g = grid.resample(edges)

As a convenience, one can also supply another Grid as the argument for this method

g = grid.resample(othergrid)

and the edges are taken from Grid.edges.

resample_factor(factor)

Resample to a new regular grid.

Parameters

factor (float) – The number of grid cells are scaled with factor in each dimension, i.e., factor * N_i cells along each dimension i.

Returns

Return type

Grid

See also

resample()

save(filename)

Save a grid object to <filename>.pickle

Internally, this calls Grid.export(filename, format="python"). A grid can be regenerated from the saved data with

g = Grid(filename="grid.pickle")

Note

The pickle format depends on the Python version and therefore it is not guaranteed that a grid saved with, say, Python 2.7 can also be read with Python 3.5. The OpenDX format is a better alternative for portability.

class MDAnalysis.analysis.density.BfactorDensityCreator(pdb, delta=1.0, atomselection='resname HOH and name O', metadata=None, padding=1.0, sigma=None)

Create a density grid from a pdb file using MDAnalysis.

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.

Parameters

• pdb (str) – PDB file or MDAnalysis.Universe;

• atomselection (str) – selection string (MDAnalysis syntax) for the species to be analyzed

• delta (float) – bin size for the density grid in Angstroem (same in x,y,z) [1.0]

• metadata (dict) – dictionary of additional data to be saved with the object

• padding (float) – increase histogram dimensions by padding (on top of initial box size)

• sigma (float) – width (in Angstrom) of the gaussians that are used to build up the density; if None (the default) then uses B-factors from pdb

Notes

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.

Examples

The following creates the density with the B-factors from the pdb file:

DC = BfactorDensityCreator(pdb, delta=1.0, atomselection="name HOH",
                           padding=2, sigma=None)
density = DC.Density()

See also

density_from_PDB()

Density(threshold=None)

Returns a Density object.

MDAnalysis.analysis.density.notwithin_coordinates_factory(universe, sel1, sel2, cutoff, not_within=True, use_kdtree=True, updating_selection=False)

Generate optimized selection for ‘sel1 not within cutoff of sel2’

Parameters

• universe (MDAnalysis.Universe) – Universe object on which to operate

• sel1 (str) – Selection string for the solvent selection (should be the group with the larger number of atoms to make the KD-Tree search more efficient)

• sel2 (str) – Selection string for the solute selection

• cutoff (float) – Distance cutoff

• not_within (bool) –

  • True: selection behaves as “not within” (As described above)

  • False: selection is a “<sel1> WITHIN <cutoff> OF <sel2>”

• use_kdtree (bool) –

  • True: use fast KD-Tree based selections

  • False: use distance matrix approach

• updating_selection (bool) – If True, re-evaluate the selection string each frame.

Notes

• Periodic boundary conditions are not taken into account: the naive minimum image convention employed in the distance check is currently not being applied to remap the coordinates themselves, and hence it would lead to counts in the wrong region.

• With updating_selection=True, the selection is evaluated every turn; do not use distance based selections (such as “AROUND”) in your selection string because it will likely completely negate any gains from using this function factory in the first place.

Examples

notwithin_coordinates_factory() creates an optimized function that, when called, returns the coordinates of the “solvent” selection that are not within a given cut-off distance of the “solute”. Because it is KD-tree based, it is cheap to query the KD-tree with a different cut-off:

notwithin_coordinates = notwithin_coordinates_factory(universe, 'name OH2', 'protein and not name H*', 3.5)
...
coord = notwithin_coordinates()        # get coordinates outside cutoff 3.5 A
coord = notwithin_coordinates(cutoff2) # can use different cut off

For programmatic convenience, the function can also function as a factory for a simple within cutoff query if the keyword not_within=False is set:

within_coordinates = notwithin_coordinates_factory(universe, 'name OH2','protein and not name H*', 3.5,
                                                   not_within=False)
...
coord = within_coordinates()        # get coordinates within cutoff 3.5 A
coord = within_coordinates(cutoff2) # can use different cut off

(Readability is enhanced by properly naming the generated function within_coordinates().)

Footnotes

1(1,2,3)

Making molecules whole can be accomplished with the MDAnalysis.core.groups.AtomGroup.wrap() of Universe.atoms (use compound="fragments").

When using, for instance, the Gromacs command gmx trjconv

gmx trjconv -pbc mol -center -ur compact

one can make the molecules whole -pbc whole, center it on a group (-center), and also pack all molecules in a compact unitcell representation, which can be useful for density generation.

2

Superposition can be performed with MDAnalysis.analysis.align.AlignTraj.

The Gromacs command gmx trjconv

gmx trjconv -fit rot+trans

will also accomplish such a superposition. Note that the fitting has to be done in a separate step from the treatment of the periodic boundaries 1.

3

Note that the trajectory in the example (XTC) is not properly made whole and fitted to a reference structure; these steps were omitted to clearly show the steps necessary for the actual density calculation.