GSoC Report

4 minute read

Published:

Introduction

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

Residue 1

if PBC are applied to this and atoms are wrapped into simulation cell the atoms would appear as below

Residue wrapped 1

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[2]

Whether a molecule splits across the periodic boundary or not, effects calculations of various properties.

Projection1

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, select_wrap1

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.

select_unwrap 1

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.

The Objectives

  1. Extend Implementation of Unwrap to molecules which are not bonded
  2. Add unwrap keyword to different methods of AtomGroup
  3. Adds ASV benchmarks to track benchmarks

Part 1: Extend Implementation of Unwrap

Pull Requests: #23322, #22682

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[2]
group.atoms.arrange_closest(universe.dimensions[:3]/2, reference='cog', inplace="True")

The resultant positions are shown below.

Periodic Image Closest to the center of cell 1

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

Pull Requests: #2275, #2279, #2286, #2288, #2290, #22952, #23052, #22982, #22962

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

Pull Requests: #22972

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.

TODO

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.

  1. Image created using OVITO  2 3 4 5 6

  2. Pull Requests yet to be merged  2 3 4 5 6 7 8