This post entails a description of the work done by me during GSoC.
The Problem Statement (Why is unwrap needed?)
Periodic Boundary Conditions(PBC) are extensively used in Molecular Dynamics to approximate a larger system using a small lattice. For example see Fig below
if PBC are applied to this and atoms are wrapped into simulation cell the atoms would appear as below
Often atoms positions are represented as above in data files.
Note the above atom positions can be obtained from the following code.
from MDAnalysisTests.datafiles import ( TRZ_psf, TRZ, ) universe = mda.Universe(TRZ_psf, TRZ) group = universe.residues
Whether a molecule splits across the periodic boundary or not, effects calculations of various properties.
The figure above shows four periodic images projected on the X-Y plane.
If the atoms are wrapped about the PBC, the atoms in red shown below are selected, 1
Now if we calculate the Moment of Interia(MoI), for the wrapped atoms we get the following value
[[2177.3019906, -237.02156344, -554.46003421] [-237.02156344, 2548.66172155, -381.71081487] [-554.46003421, -381.71081487, 2557.54815044]]
Which is wrong because it includes atoms 3 different periodic images.
However, if we unwrap the atoms, the atoms shown below are selected.
Now considering the selected positions for calculations of MoI we get the following values,
[[2820.12756093, 297.26274974, -143.17261246] [ 297.26274974, 2390.56824468, 1030.83862229] [-143.17261246, 1030.83862229, 1860.72757184]]
These values are correct as they consider an entire molecule for the calculation of MOI.
Hence we need to unwrap atoms before calculating the properties.
- Extend Implementation of Unwrap to molecules which are not bonded
- Add unwrap keyword to different methods of AtomGroup
- Adds ASV benchmarks to track benchmarks
Part 1: Extend Implementation of Unwrap
Group.unwrap was implemented by Johannes Zeman in PR #2189. This method work for molecules which are connected with bonds but not for non-bonded molecules like DNA. Part 1 of the project was to extend the current working of unwrap to work for non-bonded molecules. This is done by first unwrapping each molecule and then moving them to their own periodic image which is closer to a reference point.
PR #2332 implements arrange_closest method, which moves a molecule to the periodic image which is closest to the target_point. For example, to move the original positions to the periodic image which closest to the center of box we use.
universe = mda.Universe(TRZ_psf, TRZ) group = universe.residues group.atoms.arrange_closest(universe.dimensions[:3]/2, reference='cog', inplace="True")
The resultant positions are shown below.
PR #2268 extends the implementation of unwrap by moving molecules to the closest periodic image to the reference point after unwrapping them.
Part 2: Add unwrap keyword
The keyword unwrap allows the user to unwrap the molecules before calculating to property. I have added unwrap keyword to various methods (Each PR corresponds to a different method). An example of the usage of unwrap is shown below.
from MDAnalysisTests.core.util import UnWrapUniverse u = UnWrapUniverse(is_triclinic=False) group = u.atoms[31:39] # get center with unwrap center1 = group.center(weights=None, pbc=False, compound='segments', unwrap=True) # center1 is [[10.1000002, 7.5, 7.]] # get center without unwrap center2 = group.center(weights=None, pbc=False, compound='segments', unwrap=False) # center1=2 is [[5.1, 7.5, 7. ]]
Part 3: ASV Benchmarks
Airspeed Velocity Benchmarks allow us to track memory usage and performance of the code over different versions.
They are straightforward to implement. The only problem I faced was to handle functions which were not implemented in previous versions.
For this, I moved all the new methods to a different class and then returning an NotImplementedError error for the setup of that class. Look here to see how this is done.
Some of the pull requests (marked with 2) are yet to merged. Once they are merged, relevant benchmarks will be added in air speed velocity.