4.1.4. Generating the hopping trajectory — hop.trajectory

Based on a definition of grid sites, convert a molecular dynamics trajectory into a trajectory of site hops.

You will also need the following modules to create the input for HoppingTraj: hop.sitemap.

4.1.4.1. Classes

class hop.trajectory.HoppingTrajectory(trajectory=None, group=None, density=None, filename=None, hopdcd=None, hoppsf=None, fixtrajectory=None, verbosity=3)

Provides a time-sequence of sites visited by individual molecules, called a ‘hopping trajectory’ because the molecules hop between sites. Their coordinates are mapped to site labels, which have been defined on a grid previously (using hop.sitemap).

Output format:

For simplicity and code reusal this is again a dcd with the site as the x-coordinate; the y coordinate is set to the ‘orbit site’, i.e. it records the site the particle was last at for as long as it does not enter a new site. It describes the site in whose ‘basin of attraction’ the particle orbits. Note, however, that the transition to a new site is still counted as belonging to the previous site (which is arguably incorrect); the hop.graph module, however, does a proper analysis, which is cannot be done here for efficieny reasons. The z field is unused at the moment and set to 0.

Attributes:

ts MDAnalysis.Timestep object n_frames number of frames in hopping trajectory group AtomGroup of atoms that are tracked

Methods:

## [start:stop] object can be used as an iterator over the ## hopping trajectory (disabled du to problems when doing random ## access on large dcds; either a bug in DCDReader or python) next() advances time step in the hopping trajectory map_dcd() iterator that updates the ts and maps the trajectory

coordinates to site labels

_map_next_timestep() map next coordinate trajectory step to hopping time step _read_next_timestep() read next timestep from hopping trajectory

write() write the hopping trajectory to a dcd file + psf write_psf() write a dummy psf for visualization

Converts a trajectory into a hopping trajectory, using a sitemap as an index for sites.

>>> h = HoppingTrajectory(trajectory=DCDReader,group=AtomGroup,density=Density,
                          fixtrajectory=<dict>,verbosity=3)
>>> h = HoppingTrajectory(filename=<name>)

Create from a coordinate trajectory of a group of atoms and a site map:

u = MDAnalysis.Universe(psf,dcd) water = u.select_atoms(‘name OH2’) h = HoppingTrajectory(trajectory=u.trajectory,group=water,density=water_density)

Load from a saved hopping trajectory (in dcd format with dummy psf)

h = HoppingTrajectory(hopdcd=’hops.trajectory’,hoppsf=’hops.psf’)

Arguments:

trajectory MDAnalysis.trajectory trajectory instance group MDAnalysis.group instance density grid3Dc.Grid instance with sitemap set

hopdcd dcd written by write() hoppsf psf written by write() (or write_psf()) filename or simply provide one filename prefix for psf and dcd

fixtrajectory dictionary with attributes of a dcd object and new

values; used to provide correct values after using a catdcd-generated trajectory (hack!), e.g. fixtrajectory = {‘delta’:10.22741474887299}

verbosity show status messages for >= 3

filename(filename=None, ext=None, set_default=False, use_my_ext=False)

Supply a file name for the object.

fn = filename() —> <default_filename> fn = filename(‘name.ext’) —> ‘name’ fn = filename(ext=’pickle’) —> <default_filename>’.pickle’ fn = filename(‘name.inp’,’pdf’) –> ‘name.pdf’ fn = filename(‘foo.pdf’,ext=’png’,use_my_ext=True) –> ‘foo.pdf’

The returned filename is stripped of the extension (use_my_ext=False) and if provided, another extension is appended. Chooses a default if no filename is given. Raises a ValueError exception if no default file name is known.

If set_default=True then the default filename is also set.

use_my_ext=True lets the suffix of a provided filename take priority over a default ext(tension).

map_dcd(start=None, stop=None, step=None)

Generator to read the trajectory from start to stop and map positions to sites.

ts = map_dcd(**kwargs)

Arguments: start starting frame number (None means first) stop last frame to read (exclusive) (None means last)

(Those are arguments to dcd[start:stop].)

Iterator Returns: ts hopping trajectory timestep object (iterator)

next()

Provides the next time step of a hopping trajectory.

ts = next()

If a hopping trajectory file exists then this is used. Otherwise, the coordinate trajectory is mapped on the fly (which is computationally more expensive).

ts

Timestep of the hoptraj

write(filename, start=None, step=None, delta=None, load=True)

Write hopping trajectory as standard dcd file, together with a minimal psf.

write(‘hop’)

Arguments:

load = True Immediately loads the trajectory so that further

calls to next() will use the computed trajectory and don’t use expensive mapping.

Ignore the other options and leave them at the defaults. Currently, only the whole trajectory is written. For visualization one also needs the dummy psf of the group.

Results:

filename.trajectory and filename.psf

Note that it is your responsibility to load the hopping trajectory and the appropriate psf together as there is very limited information stored in the dcd itself.

write_psf(filename)

Write a dummy psf just for the atoms in the selected group so that one can visualize the hopping trajectory.

write_psf(filename)

The psf is NOT a fully functional psf. It only contains the header and the ATOMS section. It is sufficient to display the hopping trajectory in VMD and can be read in by the MDAnalysis tools in order to store the atom numbers for the hopping trajectory.

Format from psffres.src

CHEQ: II,LSEGID,LRESID,LRES,TYPE(I),IAC(I),CG(I),AMASS(I),IMOVE(I),ECH(I),EHA(I)

standard format:

(I8,1X,A4,1X,A4,1X,A4,1X,A4,1X,I4,1X,2G14.6,I8,2G14.6) (I8,1X,A4,1X,A4,1X,A4,1X,A4,1X,A4,1X,2G14.6,I8,2G14.6) XPLOR

expanded format EXT:

(I10,1X,A8,1X,A8,1X,A8,1X,A8,1X,I4,1X,2G14.6,I8,2G14.6) (I10,1X,A8,1X,A8,1X,A8,1X,A8,1X,A4,1X,2G14.6,I8,2G14.6) XPLOR

no CHEQ: II,LSEGID,LRESID,LRES,TYPE(I),IAC(I),CG(I),AMASS(I),IMOVE(I)

standard format:

(I8,1X,A4,1X,A4,1X,A4,1X,A4,1X,I4,1X,2G14.6,I8) (I8,1X,A4,1X,A4,1X,A4,1X,A4,1X,A4,1X,2G14.6,I8) XPLOR

expanded format EXT:

(I10,1X,A8,1X,A8,1X,A8,1X,A8,1X,I4,1X,2G14.6,I8) (I10,1X,A8,1X,A8,1X,A8,1X,A8,1X,A4,1X,2G14.6,I8) XPLOR

class hop.trajectory.TAPtrajectory(trajectory=None, group=None, TAPradius=2.8, TAPsteps=3, filename=None, dcd=None, psf=None, fixtrajectory=None, verbosity=3)

Provides a Time-Averaged Position (TAP) version of the input trajectory.

The method is described in Henchman and McCammon, J Comp Chem 23 (2002), 861 doi:10.1002/jcc.10074

Attributes:

ts MDAnalysis.Timestep object n_frames number of frames in TAP trajectory group AtomGroup of atoms that are tracked

Methods:

## [start:stop] object can be used as an iterator over the ## hopping trajectory (disabled due to dcdreader bug) next() advances time step in the hopping trajectory map_dcd() iterator that updates the ts and maps the trajectory

coordinates to site labels

_map_next_timestep() map next coordinate trajectory step to hopping time step _read_next_timestep() read next timestep from hopping trajectory

write() write the hopping trajectory to a dcd file + psf

A TAP trajectory object converts a trajectory into a TAP trajectory.

Create from a coordinate trajectory of a group of water residues:

u = MDAnalysis.Universe(psf,dcd) water = u.select_atoms(‘resname TIP*’) # see NOTE below!! water = u.select_atoms(‘name OH2’) # better, see NOTE below!! h = TAPtrajectory(trajectory=u.trajectory,group=water)

Load from a saved hopping trajectory (in dcd format with dummy psf)

h = TAPtrajectory(dcd=’TAP.trajectory’,psf=’TAP.psf’)

The given atom group is filtered according to the Time-Averaged Positon algorithm (Henchman and McCammon, J Comp Chem 23 (2002), 861). Original positions are replaced by their TAPs: A particles last position (TAP) is retained unless it has moved farther than TAPradius from its TAP measured by its root mean square distance over the last TAPsteps frames.

One can use a TAP filtered trajectory ‘on-the-fly’ to build the density:

u = Universe(psf,dcd) oxy = u.select_atoms(‘name OH2’) TAP = TAPtrajectory(u.trajectory,oxy) u.trajectory = TAP.trajectory # <— replace orig dcd with TAP !! dens = hop.density.density_from_Universe(u,atomselection=’name OH2’)

NOTE: In the current implementation residues are often ripped apart because all coordinates are processed independently. It is recommended to only do TAP on the water oxygens (for speed). This will create a trajectory in which hydrogens are always ripped from the oxygen but this trajectory is ONLY being used for creating a density from those oxygen using hop.sitemap.build_density().

(This could be fixed at the cost of speed; in this case TAP would be done on the centre of mass and the whole residue would be translated.)

Arguments:

trajectory MDAnalysis.trajectory trajectory instance group MDAnalysis.group instance (from the same Universe as trajectory) TAPradius particles are considered to be on the TAP as long as they

haven’t moved farther than TAPradius over the last TAPsteps frames

TAPsteps RMS distance of particle from TAP over TAPsteps is compared

to TAPradius

dcd dcd written by write() psf psf written by write() (or write_psf()) filename or simply provide one filename prefix for psf and dcd

fixtrajectory dictionary with attributes of a dcd object and new

values; used to provide correct values after using a catdcd-generated trajectory (hack!), e.g. fixtrajectory = {‘delta’:10.22741474887299}

verbosity show status messages for >= 3

filename(filename=None, ext=None, set_default=False, use_my_ext=False)

Supply a file name for the object.

fn = filename() —> <default_filename> fn = filename(‘name.ext’) —> ‘name’ fn = filename(ext=’pickle’) —> <default_filename>’.pickle’ fn = filename(‘name.inp’,’pdf’) –> ‘name.pdf’ fn = filename(‘foo.pdf’,ext=’png’,use_my_ext=True) –> ‘foo.pdf’

The returned filename is stripped of the extension (use_my_ext=False) and if provided, another extension is appended. Chooses a default if no filename is given. Raises a ValueError exception if no default file name is known.

If set_default=True then the default filename is also set.

use_my_ext=True lets the suffix of a provided filename take priority over a default ext(tension).

map_dcd(start=None, stop=None, skip=1)

Generator to read the trajectory from start to stop and map positions to TAP sites.

ts = map_dcd(**kwargs)

Arguments: start starting frame number (None means first) stop last frame to read (exclusive) (None means last)

(Those are arguments to dcd[start:stop].)

Iterator Returns: ts hopping trajectory timestep object (iterator)

next()

Provides the next time step of a TAP trajectory.

ts = next()

If a TAP trajectory file exists then this is used. Otherwise, the coordinate trajectory is mapped on the fly (which is computationally more expensive).

write(filename, start=None, step=None, delta=None, load=True)

Write hopping trajectory as standard dcd file.

write(‘TAP’)

Arguments:

load = True Immediately loads the trajectory so that further

calls to next() will use the computed trajectory and don’t use expensive mapping.

Ignore the other options and leave them at the defaults. Currently, only the whole trajectory is written. All atoms in the original trajectory are written to the output so you should be able to use your original psf file.

NOTE: Fixed atoms are possibly not accounted for properly.

Note that it is your responsibility to load the TAP trajectory and the appropriate psf together as there is very limited information stored in the dcd itself.

class hop.trajectory.ThinDCDReader(datafeeder)

DCD-like object that supports a subsection of the DCDReader interface such as iteration over frames and most attributes. The important part is that the __iter__() method is overriden to provide data from another source. This allows a filter architecture for trajectories.