Question on distance analysis

[horlust]

Hey everyone!

Just a quick question, couldn't find the answer in the manual. If I'm to compute the distance between two molecules, for example:

dgluol = u.select_atoms("resname DGLUOL")
    waters = u.select_atoms("resname TIP3")
    sizes = []
    
    for _ in tqdm(u.trajectory, total=len(u.trajectory)):
        dists = distance_array(dgluol.positions, waters.positions, box=u.dimensions)

Is this distance between the centers of mass of the two molecules? Or am I completely misusing the tool? Lol

Thanks!

[orbeckst]

No, it's not center of mass distance.

distance_array(A, B) operates on raw coordinates, namely it just calculates the "all-vs-all" distances of all coordinates A from all coordinates B.

You should first calculate the COMs and then use distance_array on the COMs

water = u.select_atoms("resname SOL and not name MW")  # exclude mass-less particles in TIP4P
protein = u.select_atoms("protein")


for ts in u.trajectory:
   w_com = water.center_of_mass(compound="residues", unwrap=True)
   protein_residues_com = protein.center_of_mass(compound="residues", unwrap=True)
   
   # produces arrays of shape (n_residues, 3) with the COM (correct under PBC) for each residue
   assert w_com.shape == (water.n_residues, 3)
   assert protein_residues_com.shape == (protein.n_residues, 3)

   # calc distances (under PBC with minimum image)
   dists = distance_array(protein_residues_com, w_com, box=ts.dimensions)

Calculating the centers of mass with unwrap=True is a bit slow. If your trajectory is already unwrapped (i.e., the residues of interest are whole and not broken across the periodic boundaries) then you don't need that.