Update 1

3 minute read

Published:

Introduction

I will cover my work updates and things I have learned during the project in this post.

Updates

  • PR #2268 adds reference_point keyword in the unwrap method. The tests for cubic systems have been updated, and I am currently updating the tests for hexagonal systems.
  • PR #2275 adds unwrap keyword to center. It is currently being reviewed.

The Unwrap Universe

UnWrapUniverse was created by Johannes Zeman to test the unwrap method. An approximate projection onto the x-z plane of box is shown below

         :                   :
         :                   :
    6    :        8          :                   :
  - + - -:+-------+----------:- - - - - - - - - -: -
         |5       8          |                   :
         |                   |                   :
         |         (10)      |                   :
         |   o-o-o-o-x-x-x-x |                   :
         +x-x-x-x     o-o-o-o+                   :
         |            (11)   |                   :
         |                   |                   :
         |         1         |       2           : 3
         |                   |                   :
         |          9        |                   :
         |   4      !       7|                   :
         |   !      9       !|                   :
         |   4-4            7+7                  :
         |                   |                   :
        5+5                  |                   :
  - + - -:+------------------:- - - - - - - - - -: -
    6-6  :                   :                   :
         :                   :
         :                   :
         :                   :
         :             (12)  :
         :        x-x-x-x-o-o+o-o
         :                   :

UnWrapUniverse contains 47 atoms in 12 molecules. In this post, we will consider molecule 10 and molecule 11. Molecule 10 lies entire in the unit cell while molecule 11 is split across the surface due to periodic boundary conditions.

Both of them have 8 atoms, and their positions are shown below

molecule_10_positions =   [[ 2.,  7.,  8.],
                           [ 3.,  7.,  8.],
                           [ 4.,  7.,  8.],
                           [ 5.,  7.,  8.],
                           [ 6.,  7.,  8.],
                           [ 7.,  7.,  8.],
                           [ 8.,  7.,  8.],
                           [ 9.,  7.,  8.]]

molecule_11_positions =   [[ 6.6,  7.5,  7. ],
                           [ 7.6,  7.5,  7. ],
                           [ 8.6,  7.5,  7. ],
                           [ 9.6,  7.5,  7. ],
                           [ 0.6,  7.5,  7. ],
                           [ 1.6,  7.5,  7. ],
                           [ 2.6,  7.5,  7. ],
                           [ 3.6,  7.5,  7. ]]

Now, if we were to find the center of each molecule by taking the average the positions of their atoms, we would get the following values.

molecule_10_center = [4.75, 7, 8]

molecule_11_center = [5.1, 7.5, 7]

MDAnalysis has a method center to find the above result. While the value of the center is correct for molecule 10, the value for molecule 11 is wrong. In order to get correct result for center the positions of atoms of molecule 11 need to be unwrapped before averaging them. This is done by adding a keyword “unwrap”. When it is set to True molecules will be unwrapped before calculating center. The working on unwrap keyword is show in the next section.

Current Working of center

from MDAnalysisTests.core.util import UnWrapUniverse

u = UnWrapUniverse(is_triclinic=False)
group = u.atoms[31:39]  # molecule 11


# 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. ]]

Lessons

  • Decorators

    Thanks to Richard Gowers, I ended up implementing decorators for filtering the values for pbc and unwrap keywords. For compound type ‘group’, both ‘unwrap’ and ‘pbc’ cannot be True. This decorator will be used will adding unwrap to other functions, hence will help in avoiding repetition in the code.

  • Pytest Pytest has been quite a useful testing tool. @pytest.mark.parametrize especially helps in adding multiple combinations of tests. For example, just adding reference_point at Line 233 tests different combinations of reference_point with varying values of compound and level.

  • Working with cython

    This was the first time I had to code in cython. One of the major error I made was not passing dtype in numpy array. I had to spend a lot of time figuring out why I was getting random large numbers in my array. A good lesson to learn.