On computing distances

[horlust]

Hi!

I have a simulation containing water and two solutes and I want to compute the distance between these solutes over time. I have accounted for PBC effects with GROMACS. A screenshot of the simulation box:

Because of PBC effects, the two molecules are actually close, but when I compute the distances with MDAnalysis I get bigger distances than expected. Is there a way to compute the "smallest" distance between these two solutes?

My code is:

u = mda.Universe("prod.gro", "prod.xtc")

group1 = u.select_atoms("resname MCD") 
group2 = u.select_atoms("resname UNL") 

distances = []

for ts in u.trajectory:
    com1 = group1.center_of_mass(pbc=True)  # COM of MCD
    com2 = group2.center_of_mass(pbc=True)  # COM of UNL
    distance = np.linalg.norm(com1 - com2)  # distance between COMs
    distances.append(distance)

Thanks!

[orbeckst]

Use the PBC-enabled distance functions in MDAnalysis.lib.distances (or for very simple calculations, analysis.distances.AtomicDistances) and do not use np.linalg.norm or any other naive distance calculation.

The important part is to calculate distances under minimum image convention by providing the box=ts.dimensions argument.

Below is what I would try — not tested, take it as a starting point for exploring.

COM distance

Calculate the COM (as in your question) and the compute the min-image distance between COM1 and COM1:

u = mda.Universe("prod.gro", "prod.xtc")

group1 = u.select_atoms("resname MCD") 
group2 = u.select_atoms("resname UNL") 

distances = []

for ts in u.trajectory:
    com1 = group1.center_of_mass(pbc=True)  # COM of MCD
    com2 = group2.center_of_mass(pbc=True)  # COM of UNL
    d = MDAnalysis.lib.distances.calc_bonds(com1, com2, box=ts.dimensions)
    distances.append(d)

closest interatomic distance

Calculate all pairwise distances ( $d_{ij} = |\mathbf{r}^\mathrm{MCD}_i - \mathbf{r}^\mathrm{UNL}_j|_\mathrm{min image}$ ) with minimum image convention and then take the minimum $d = \mathrm{min}_{i,j} d_{ij}$):

u = mda.Universe("prod.gro", "prod.xtc")

group1 = u.select_atoms("resname MCD") 
group2 = u.select_atoms("resname UNL") 

distances = []

for ts in u.trajectory:
    all_d = MDAnalysis.lib.distances.distance_array(group1, group2, box=ts.dimensions)
    shortest_d = all_d.min()
    distances.append(shortest_d)