Atom order rearrangemet for [Water dynamics analysis]

[Ji Woong Yu]

Dear MDAnalysis users,

I am trying to use the  Water dynamics analysis Module as implemented in MDAnalysis package (branch 2.4.0-dev).

In my code, the only input is a LAMMPSDUMP trajectory file with enormous size and the atom id is just randomly given. It means one cannot expect atoms in a molecule don't have consecutive IDs. The topology can be reconstructed by linking the two nearest hydrogens from an oxygen atom.

For hydrogen bond analysis, I can see that it checks the charge and mass of an atom to discern whether the atom is hydrogen or oxygen (or an atom with high electronegativity).

However, when I looked into the source code, I got to know that the way it reads the atoms from the data from MDAnalysis.universe is just assuming that the given atom group is ordered in O(molecule 1), H(molecule 1), H(molecule 1), O(molecule 2), H(molecule 2), H(molecule 2), ... (Please correct me if I am wrong)

This is inferred from the code below (mdanalysis/package/MDAnalysis/analysis/waterdynamics.py: line456~460):
for j in range(len(repInd[i]) // 3):
begj = 3 * j
universe.trajectory[t0]
Ot0 = repInd[i][begj]
H1t0 = repInd[i][begj + 1]
H2t0 = repInd[i][begj + 2]

AFAIK, the index cannot be changed in MDAnalysis.

I think the most straightforward way to fix this problem is to rearrange the trajectory file (i.e. reassigning atom ID). However, the size is too large so I want to avoid it as far as I can.

Is there a convenient way to solve this?

Thanks and Best Regards,

Jiwoong Yu

[Oliver Beckstein]

Hi Jiwoong Yu,

welcome to the MDAnalysis mailing list!

The water dynamics has the limitation that you point out: it unfortunately relies on the ordering of atoms. I might be possible to change it; please open an issue in our issue tracker https://github.com/MDAnalysis/mdanalysis/issues.

I cannot think of a simple way to change the index, there’s certainly no way to just relabel the ID.

You could reorder the information in your timestep using on the fly transformations. That’s more advanced and potentially messy but might work. See https://userguide.mdanalysis.org/stable/trajectories/transformations.html and https://docs.mdanalysis.org/stable/documentation_pages/trajectory_transformations.html 

If you have an array of new indices, let’s say new_indices where

new_indices[0] = index of O atom in water 0 that should come at index 0

new_indices[1] = index of H1 atom in water 0 that should come at index 1

new_indices[2] = index iof H2 atom in water 0 that should come at index 2

new_indices[3] = index of O atom in water 1 that should come at index 3

new_indices[4] = index of H1 atom in water 1 that should come at index 4

…

i.e., a permutation of the original indices so that you have the O H H O H H … ordering, then you could use a transformation

def reorder_waters(ts, new_indices=new_indices):

   ts.positions = ts.positions[new_indices]

   return ts

and then apply the transformation to the universe

u.add_transformations(reorder_water)

This is will ONLY change the coordinates around. Any topology information is still in the original order.  I can’t guarantee that this will work with water dynamics. However, this might give you an idea how to get started if you want to change coordinate ordering on the fly.

Oliver

--
You received this message because you are subscribed to the Google Groups "MDnalysis discussion" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discus...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/1b245265-4962-401c-8b52-ea084f0b9fb7n%40googlegroups.com.

--

Oliver Beckstein (he/his/him)

email: orbe...@mdanalysis.org
twitter: @orbeckst

GitHub: @orbeckst

MDAnalysis – a NumFOCUS fiscally sponsored project

https://www.mdanalysis.org/

[Ji Woong Yu]

Dear Oliver,

Thank you for welcoming me and the answer,

(1) I made the request (https://github.com/MDAnalysis/mdanalysis/issues/3897)

(2)
>>> def reorder_waters(ts, new_indices=new_indices):

>>>    ts.positions = ts.positions[new_indices]

>>>    return ts

>>> 

>>> and then apply the transformation to the universe

>>> 

>>> u.add_transformations(reorder_water)

>>> 

>>> This is will ONLY change the coordinates around. Any topology information is still in the original order.  I can’t guarantee that this will work with water dynamics. However, this might >>> give you an idea how to get started if you want to change coordinate ordering on the fly.

One more quick question.

After a quick skimming through the water dynamics code (line 859~864 of mdanalysis/package/MDAnalysis/analysis/waterdynamics.py), I think the subroutine which actually selects the water is implemented inside the code. 

""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
def _selection_serial(self, universe, selection_str):
    selection = []
    for ts in ProgressBar(universe.trajectory, verbose=True,
                                      total=universe.trajectory.n_frames):

        selection.append(universe.select_atoms(selection_str))

    return selection

""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

The only thing user can do is passing the selection string to the module from the outside. Still, do you expect the trick (on-the-fly the change of coordinates) do the work?

Anyway, I'll give it a try if your answer actually does work.

Thanks and Best Regards,

Jiwoong Yu

2022년 10월 31일 월요일 오전 11시 39분 36초 UTC+9에 orbe...@mdanalysis.org님이 작성:

[Oliver Beckstein]

On Oct 30, 2022, at 11:03 PM, Ji Woong Yu <jiwoon...@gmail.com> wrote:

(1) I made the request (https://github.com/MDAnalysis/mdanalysis/issues/3897)

Thank you.

The only thing user can do is passing the selection string to the module from the outside. Still, do you expect the trick (on-the-fly the change of coordinates) do the work?

To be honest, I would be surprised if it worked out of the box. I think you’d also need to relabel the atoms and residues so that the selections work. You’d need to relabel in such a way that the selections point to the right atoms. You can load the universe, add transformation, relabel resids and types (or whatever is needed for the selection). If you show code then we can comment and help.

By far the better solution is to change the waterdynamics code as per #3897. (We always welcome contributions from the community and we have some material on how to contribute to MDAnalysis https://userguide.mdanalysis.org/stable/contributing.html — just in case you start changing the code and want your changes become part of MDAnalysis. We’re happy to help when you have questions.)

Best,

Oliver