Using Simple atomic distance analysis for atom groups in two different trajectories

[ullaharslan]

Hello Everyone, I am not a very expert user of MDAnalysis, Here I have a question regarding 'simple atomic distance analysis' as explained in the documentation, it says " This module provides a class to efficiently compute distances between two groups of atoms with an equal number of atoms over a trajectory".Since this method is bound to calculate the distances over a trajectory if offer another bound that both gorups of atoms should belong same trajectory/frame. I tried to slice the trajectory following documentation, I accessed two different trajectories and tried to calculate the distance but result in just the zeros. I also tried to extract coordinated info as numpy array and process but it resulted in loss of periodic boundary condition(pbc) thus effecting eventual results. Can anyone suggest a smart way to select atom groups from first trajectory u.trajectory[0] and subtract distances from all over the trajetory ? does Time step class can help in this regard ?
please accept my best regards
Irslan Ullah Ashraf

[orbeckst]

Can you show your code?

Can you explain again what you want to calculate, e.g., by writing down equations? Do you want to calculate the displacement of atom i relative to frame 0 as a function of time

d_i(t) = x_i(t)  - x_i(0)

where x_i is the position of atom i.
Or do you want the distance between atom i at time t and atom j at time 0,

d_ij(t) = x_i(t)  - x_j(0)

You'll probably do something like (for the first case) using calc_bonds()

import numpy as np
import MDAnalysis as mda
from MDAnalysis.lib.distances import calc_bonds

u = mda.Universe(TOP, TRAJ)
ag = u.select_atoms("protein and name CA")  # replace with your own selection 

# frame 0
u.trajectory[0]
x_ref = ag.positions  # makes a copy as a numpy array

d_i = np.zeros((u.trajectory.n_frames, ag.n_atoms), dtype=np.float64)
for iframe, ts in enumerate(u.trajectory):
   # calculate distances with PBC taken into account (use ts.dimensions)
   d_i[iframe] = calc_bonds(x_ref, ag, box=ts.dimensions)

The array d_i now contains the displacements for the atoms in ag for all time steps.

[ullaharslan]

Dear @orbeckst Thank you so much for kind attention and wonderful suggestion, I would definitely add equations but your mentioned equations d_i(t) for atom i at time 0 and t, and d_ij(t) for atoms i and j in the respective time are absolutely correct(with only exception that d_ij(t) refers to atomgroup rather than single atom distance, and this is obvious). Here is my code where I extracted coordinate information for post processing after some failed attempt.

import MDAnalysis as mda
import numpy as np
u=mda.Universe('data_all_1800_1ns','dump_all_1800_1ns', topology_format="data", format="lammpsdump",atom_style='id type x y z',timeunit='ps',dt=0.0001)
#here i selected reference atom positions from trajec[1], I selected a desired region (sphere) to save computation
f1=u.trajectory[1]
sphere1=u.select_atoms('point 4.37956  45.9854  116.058 115',updating = False) 
u.trajectory.dt
pos1=sphere1.positions
#than i selected atomic positions at a later time of trajec[100]
f2=u.trajectory[100]
sphere2=u.select_atoms('point 4.37956  45.9854  116.058 115',updating = False) 
u.trajectory.dt
pos2=sphere1.positions
#at first i tried to calculate one atom distance using "my_dists = ad.AtomicDistances(ag1, ag2).run()" but it just gave me zero as it have to(understood) than I simply used numpy to extract the displacement
pos_diff=pos1-pos2
print(pos_diff)
# manual Calculation for magnitude for each displacement vector
magnitudes = np.sqrt(np.sum(pos_diff ** 2, axis=1))
data = np.column_stack((pos_diff, magnitudes))
np.savetxt('disp.txt',data, fmt='%.18f')

this code didn't solved my problem because

1. It didn't considered PBC and involved specific processing of frames(solved after your suggestion Thanks)
2. I also need atom type and atom coordinates information for atom d_i(t) or d_j(t) along with their displacement magnitude at respective times, is there a built in method to store this information ? I want to build something like "atom_type d_x(t) d_y(t) d_z(t) |d_i(t)|" I wonder If I can use MDAnalysis.analysis.distances.dist(A, B, offset=0, box=None) for this purpose.

[orbeckst]

In your code, sphere1 and sphere2 are probably atomgroups with different numbers of atoms (look at sphere1.n_atoms (or len(sphere1)) and sphere2.n_atoms. If there are different number of atoms in each group then I don't know what it would mean to calculate the distance between the two groups.

You could, however, try to match by atom id and only calculate distances for the same atom. Here I am using np.isin() with the atom indices to find those atoms that are in both groups

from MDAnalysis.lib.distances import calc_bonds, minimize_vectors

# ... use your code where you create sphere1 and sphere2 ...

in1 = np.isin(sphere1.ix, sphere2.ix)   # indices of 1 in 2
in2 = np.isin(sphere2.ix, sphere1.ix)   # indices of 2 in 1

# in1 and in2 are boolean arrays that select the common atoms from sphere1 and sphere2;
# they have different sizes but contain the SAME number of True elements (in1.sum() adds 
# 0 for False and 1 for True and thus counts total number of True): let's check
assert in1.sum() == in2.sum()

# get positions of matching atoms
x1 = pos1[in1]
x2 = pos2[in2]

# calculate distance with PBC
# (note that I am not 100% sure what the correct choice for box is: this may depend on your problem;
# both sphere1.dimensions and sphere2.dimensions currently contain the dimensions of the timestep 
# u.trajectory[100] so you can equally say sphere1.dimensions: 
# atomgroup coordinates and box information always update with the trajectory!)
d = calc_bonds(x1, x2, box=sphere2.dimensions)

# calculate connecting vectors between x1 and x2 (with PBC)
dxyz = minimize_vectors(x1 - x2, box=sphere2.dimensions)

# you can get all information for these common atoms
common = ag1[in1]    

# common is the same as ag2[in2], so we can double check
assert common == ag2[in2]

# You can now get information about the common atoms
print("atom indices", common.indices)
print("atom names", common.names)
print("atom type", common.types)

# combine everything in a table with
# atomtype, dx, dy, dz, |d|
table = list(zip(common.types, dxyz[:, 0], dxyz[:, 1], dxyz[:, 2], d))

The minimize_vectors() function is the equivalent of calc_bonds() but for vectors.

You may actually save some time by only calculating dxyz as above and the calculating the lengths as d = np.linalg.norm(dxyz, axis=1) instead because then you don't have to do PBC fixes twice.

There are probably better ways to organize the data but the above should produce what you asked for.

Note that this code only shows how to do a single calculation between frame index 1 (btw, why did you not go to the first frame in the trajectory u.trajectory[0] at frame index 0?) and frame index 100. In order to generalize it for a trajectory you have to decide which frame is a reference and then you might want to make the second selection an updating AtomGroup. However, you have to decide what scientific questions you want to answer. At least the code above should point you towards tools in MDAnalysis that you can use.

I wonder If I can use MDAnalysis.analysis.distances.dist(A, B, offset=0, box=None) for this purpose.

The code above does pretty much the same and instead of providing residueIDs it gives you access to whatever you want to know about the common atoms. (Just look at the source code for dist() — it's a very naive function).

Also note that you can do interesting set operations on AtomGroups: for instance

common = sphere1 & sphere2

would have also worked (but we still need in1 and in2 to select the appropriate coordinates from the coordinate arrays so I extracted common manually to show how to do it).

[ullaharslan]

Dear @orbeckst
Sorry for late response as I was busy with exam, Thank you so much for in detail and generous response, this astonishingly revolutionised me the powerful build of MDAnalysis.Now I can see number of options to perform my job, but I will think and plan to write the smart possible option. Your reply is so much helpful to me. Have a wonderful day

[orbeckst]

I am glad this was helpful. If it answered your question please mark it as an answer as this helps everybody else who's looking at the forum. Otherwise please let us know what is still not working. Thanks.