# Centering a trajectory in the box¶

Here we use MDAnalysis transformations to make a protein whole, center it in the box, and then wrap the water back into the box. We then look at how to do this on-the-fly.

Last updated: July 3, 2020 with MDAnalysis 1.0.0

Minimum version of MDAnalysis: 1.0.0

Packages required:

[1]:

import MDAnalysis as mda
from MDAnalysis.tests.datafiles import TPR, XTC
import numpy as np
import nglview as nv


The test files we will be working with here are trajectories of a adenylate kinase (AdK), a phosophotransferase enzyme. ([BDPW09])

For the step-by-step transformations, we need to load the trajectory into memory so that our changes to the coordinates persist. If your trajectory is too large for that, see the on-the-fly transformation section for how to do this out-of-memory.

[2]:

u = mda.Universe(TPR, XTC, in_memory=True)


## Before transformation¶

If you have NGLView installed, you can view the trajectory as it currently is.

[3]:

view = nv.show_mdanalysis(u)
view.center()
view


Otherwise, we embed a gif of it below.

For easier analysis and nicer visualisation, we want to center this protein in the box.

## Unwrapping the protein¶

The first step is to “unwrap” the protein from the border of the box, to make the protein whole. MDAnalysis provides the AtomGroup.unwrap function to do this easily. Note that this function requires your universe to have bonds in it.

We loop over the trajectory to unwrap for each frame.

[4]:

protein = u.select_atoms('protein')

for ts in u.trajectory:
protein.unwrap(compound='fragments')


As you can see, the protein is now whole, but not centered.

[ ]:

unwrapped = nv.show_mdanalysis(u)
unwrapped.center()
unwrapped


Over the course of the trajectory it leaves the box.

## Centering in the box¶

The next step is to center the protein in the box. We calculate the center-of-mass of the protein and the center of the box for each timestep. We then move all the atoms so that the protein center-of-mass is in the center of the box.

If you don’t have masses in your trajectory, try using the center_of_geometry.

[5]:

for ts in u.trajectory:
protein_center = protein.center_of_mass(pbc=True)
dim = ts.triclinic_dimensions
box_center = np.sum(dim, axis=0) / 2
u.atoms.translate(box_center - protein_center)


The protein is now in the center of the box, but the solvent is likely outside it, as we have just moved all the atoms.

[ ]:

centered = nv.show_mdanalysis(u)
centered.center()
centered


## Wrapping the solvent back into the box¶

Luckily, MDAnalysis also has AtomGroup.wrap to wrap molecules back into the box. Our trajectory has dimensions defined, which the function will find automatically. If your trajectory does not, or you wish to use a differently sized box, you can pass in a box with dimensions in the form [lx, ly, lz, alpha, beta, gamma].

[6]:

not_protein = u.select_atoms('not protein')

for ts in u.trajectory:
not_protein.wrap(compound='residues')


And now it is centered!

[ ]:

wrapped = nv.show_mdanalysis(u)
wrapped.center()
wrapped


## Doing all this on-the-fly¶

Running all the transformations above can be difficult if your trajectory is large, or your computational resources are limited. Use on-the-fly transformations to keep your data out-of-memory.

Some common transformations are defined in MDAnalysis.transformations.

[7]:

import MDAnalysis.transformations as trans


[8]:

u2 = mda.Universe(TPR, XTC)

protein2 = u2.select_atoms('protein')
not_protein2 = u2.select_atoms('not protein')


From version 1.0.0 onwards, the MDAnalysis.transformations module contains wrap and unwrap functions that correspond with the AtomGroup versions above. You can only use add_transformations once, so pass them all at the same time.

[9]:

transforms = [trans.unwrap(protein2),
trans.center_in_box(protein2, wrap=True),
trans.wrap(not_protein2)]


[11]:

otf = nv.show_mdanalysis(u2)