MDAnalysis comparing position of ion in trajectory to previous ts

Question

So I have a system where I need to be able to determine the exact position of my ions and run an equation on the average position of that ion. I found my ion positions were inconsistent due to some ions wrapping across the periodic boundary and severely changing the position for that one window. Leading me to have an average of say +20 when the ion just shuffled between +40 and -40.

I was wanting to correct that by implementing a way to unwrap my wrapped coordinates for ions on the edge of my box.

Essentially I was thinking that for each frame in my trajectory, MDAnalysis would check the position of ION 1 in frame 1. Then in frame 2 it would check the same ion once more and compare it to the previous position. If it for example goes from + coordinates to - coordinates then I would have a count that adds +1 meaning that it wrapped once. If it goes from - to + I would have it subtract 1. Then by the end of all of the frames I would have a number that could help me identify how I could perform my analysis.

However my coding skills are less than lackluster and I wanted to know how I would go about implementing this? I have essentially gotten the count down, but the comparison between frames is where I am confused. How would I do this comparison?

Thanks in advance

Answer 1

There are a few ways to answer this question. Firstly,

Essentially I was thinking that for each frame in my trajectory, MDAnalysis would check the position of ION 1 in frame 1. Then in frame 2 it would check the same ion once more and compare it to the previous position. If it for example goes from + coordinates to - coordinates then I would have a count that adds +1 meaning that it wrapped once. If it goes from - to + I would have it subtract 1. Then by the end of all of the frames I would have a number that could help me identify how I could perform my analysis.

You could write your own analysis class. One untested way to do it is prototyped below -- the tutorial goes more into what each method (_prepare, _conclude, etc) does.

from MDAnalysis.analysis.base import AnalysisBase
import numpy as np

class CountWrappings(AnalysisBase):

    def __init__(self, universe, select="name NA"):
        super().__init__(universe.universe.trajectory)
        # these are your selected ions
        self.atomgroup = universe.select_atoms(select)
        self.n_atoms = len(self.atomgroup)

    def _prepare(self):
        # self.results is a dictionary of results
        self.results.wrapping_per_frame = np.zeros((self.n_frames, self.n_atoms), dtype=bool)
        self._last_positions = self.atomgroup.positions
    
    def _single_frame(self):
        # does sign change for any element in 2D array?
        compare_signs = np.sign(self.atomgroup.positions) == np.sign(self._last_positions)
        sign_changes_any_axis = np.any(compare_signs, axis=1)
        # _frame_index is the relative index of the frame being currently analyzed
        self.results.wrapping_per_frame[self._frame_index] = sign_changes_any_axis
        self._last_positions = self.atomgroup.positions
    
    def _conclude(self):
        self.results.n_wraps = self.results.wrapping_per_frame.sum(axis=0)


n_wraps = CountWrappings(my_universe, select="name NA CL MG")
n_wraps.run()
print(n_wraps.results.wrapping_per_frame)
print(n_wraps.results.n_wraps)

However, I'm not sure that addresses your actual aim:

I was wanting to correct that by implementing a way to unwrap my wrapped coordinates for ions on the edge of my box.

Are you computing the ion positions relative to anything? Potentially you could add bonds between each ion and the center so that you can use the AtomGroup.unwrap() function. Alternatively, is your data compatible with GROMACS? GROMACS has an unwrapping utility called "nojump" that unwraps atoms jumping across box edges, e.g.

gmx trjconv -f my_trajectory.xtc -s my_topology.gro -pbc nojump -o my_unwrapped_trajectory.xtc

Answer 2

As Lily mentioned, you could write your own analysis to do this or use GROMACS. However, both Lily's example and the GROMACS implementation of 'nojump' fail to account for box size fluctuations under the NPT ensemble (assuming you've used NPT). von Bulow et al. wrote about this widespread problem a couple of years ago. As far as I'm aware, the only implementation of nojump unwrapping that accounts for box size fluctuations is in LiPyphilic (disclaimer: I am the author of LiPyphilic).

Using LiPyphilic, you can unwrap your trajectory like so:

import MDAnalysis as mda
from lipyphilic.transformations import nojump

u = mda.Universe(pdb, xtc)
ions = u.select_atoms('name NA CLA')

u.trajectory.add_transformations(
    nojump(
        ag=ions,
        nojump_x=True,
        nojump_y=True,
        nojump_z=True)
)

Then, when you do further analysis with your MDAnalysis Universe, the atoms will automatically be unwrapped at each frame.