Generating Densities from Trajectories — pmda.density

This module contains parallel versions of analysis tasks in MDAnalysis.analysis.density.

class pmda.density.DensityAnalysis(atomgroup, delta=1.0, atomselection=None, metadata=None, padding=2.0, updating=False, parameters=None, gridcenter=None, xdim=None, ydim=None, zdim=None)[source]

Parallel density analysis.

The trajectory is read, frame by frame, and the atoms selected with atomselection are histogrammed on a grid with spacing delta.

Parameters

  • atomgroup (AtomGroup) – Group of atoms (such as all the water oxygen atoms) being analyzed. If this is an updating AtomGroup then you need to set ‘atomselection’ and also ‘updating=True’.

  • 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 ångström (same in x,y,z). 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 ångström. Padding is ignored when setting a user defined grid.

  • updating (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

  • 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 ångström.

  • 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”.

See also

MDAnalysis.analysis.density.density_from_Universe

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 updating=True.

For more details about density calculations, refer to [Awtrey2019].

Examples

A common use case is to analyze the solvent density around a protein of interest. The density is calculated with DensityAnalysis 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, first create the DensityAnalysis object by supplying an AtomGroup, then use the run() method:

from pmda.density import DensityAnalysis
U = Universe(TPR, XTC)
ow = U.select_atoms("name OW")
D = DensityAnalysis(ow, delta=1.0)
D.run()
D.convert_density('TIP4P')

The positions of all water oxygens are histogrammed on a grid with spacing delta = 1 Å. Initially the density is measured in inverse cubic angstroms. 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). In particular, the density is stored as a NumPy array in grid, which can be processed in any manner. Results are available through the density attribute, which has the grid attribute that contains the histogrammed density data. DensityAnalysis.density is a gridData.core.Grid object. In particular, its contents can be exported to different formats. For example, to write a DX file water.dx that can be read with VMD, PyMOL, or Chimera:

D.density.export("water.dx", type="double")

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

D = DensityAnalysis(universe.atoms, 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 updating keyword argument:

atomselection = 'name OW and around 5 (resid 156 157 305)'
D_site = DensityAnalysis(universe.atoms, atomselection=atomselection,
                         updating=True)

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

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

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

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.

New in version 0.3.0.

Parameters

  • Universe (Universe) – a MDAnalysis.core.groups.Universe (the atomgroups must belong to this Universe)

  • atomgroups (tuple of AtomGroup) – atomgroups that are iterated in parallel

_results

The raw data from each process are stored as a list of lists, with each sublist containing the return values from pmda.parallel.ParallelAnalysisBase._single_frame().

Type

list

static current_coordinates(atomgroup, atomselection, updating)[source]

Retrieves the current coordinates of all atoms in the chosen atom selection. Note: currently required to allow for updating selections

readonly_attributes()

Set the attributes of this class to be read only

Useful to avoid the class being modified when passing it around.

To be used as a context manager:

with analysis.readonly_attributes():
    some_function(analysis)
run(start=None, stop=None, step=None, n_jobs=1, n_blocks=None)

Perform the calculation

Parameters

  • start (int, optional) – start frame of analysis

  • stop (int, optional) – stop frame of analysis

  • step (int, optional) – number of frames to skip between each analysed frame

  • n_jobs (int, optional) – number of jobs to start, if -1 use number of logical cpu cores. This argument will be ignored when the distributed scheduler is used

  • n_blocks (int, optional) – number of blocks to divide trajectory into. If None set equal to n_jobs or number of available workers in scheduler.

See also

MDAnalysis.analysis.density.density_from_Universe