PeriodicBox fuzzing

[aster31]

I have noticed that you have a PeriodicBox that supports triclinic boxes.

For context, this is something I have tried to do as well for one of my rust projects (currently not available publicly), and when I tried to fuzz it against mdtraj or mdanalysis, I have noticed that my implementation was behaving incorrectly. Thus, when I noticed that you also supported triclinic boxes, I have decided to fuzz your implementation as well, hoping that I could maybe copy it if it worked. I have noticed a potential problem, but I'm not entirely sure yet if I am just using your library wrong, or if I was fuzzing it wrong, or if this is indeed a bug. While the periodic box in my fuzzing attempt is far from orthogonal, I have made sure to obey the triclinic box conditions as described here: https://manual.gromacs.org/2021/reference-manual/algorithms/periodic-boundary-conditions.html

#!/usr/bin/env python3

import numpy as np
import pymolar
import mdtraj
import MDAnalysis as mda

pbc = np.array([[
    [10., 0., 0.],
    [4., 10., 0.],
    [-4., 0., 10.]
]])

xyz = np.array([[
    [38.9214, 40.0078, -34.0795],
    [-26.6187, 40.8926, 30.9709]
]])

traj = mdtraj.Trajectory(xyz, None)
traj.unitcell_vectors = pbc
dist_mdtraj = mdtraj.compute_distances(traj, [(0, 1)])[0][0]

p = pymolar.PeriodicBox(pbc[0])
dist_pymolar = p.distance(xyz[0, 0], xyz[0, 1], [True, True, True])

mda_box = [l for l in traj.unitcell_lengths[0]] + [a for a in traj.unitcell_angles[0]]

dist_mda = mda.lib.distances.distance_array(xyz[0, 0], xyz[0, 1], box=mda_box)[0][0]

print("points:")
print(xyz[0, 0], xyz[0, 1])

print("dist (mdtraj, mda, pymolar):")
print(dist_mdtraj, dist_mda, dist_pymolar)


closest = p.closest_image(xyz[0][0], xyz[0][1])
print("supposed closest image:", closest)

lowest = (0, 0, 0)
lowest_dist = np.linalg.norm(xyz[0, 1] - xyz[0, 0])
closest_image = xyz[0, 0]
# try to improve on closest image manually
for i in range(-10, 11):
    for j in range(-10, 11):
        for k in range(-10, 11):
            new_image = xyz[0, 0] + i * pbc[0, 0] + j * pbc[0, 1] + k * pbc[0, 2]
            dist = np.linalg.norm(xyz[0, 1] - new_image)
            if dist < lowest_dist:
                lowest_dist = dist
                lowest = (i, j, k)
                closest_image = new_image

print("manual lowest dist:", lowest_dist)
print("manual lowest offset:", lowest)
print("closest image:", closest_image)

I get the following output:

points:
[ 38.9214  40.0078 -34.0795] [-26.6187  40.8926  30.9709]
dist (mdtraj, mda, pymolar):
5.3536267 5.353625200162365 5.462098121643066
supposed closest image: [-31.078602  42.0078    33.920498]
manual lowest dist: 5.353626734280232
manual lowest offset: (-4, 0, 6)
closest image: [-25.0786  40.0078  25.9205]

I was using the latest pymolar on PyPI, 1.3.1 I think

[yesint]

Hi, @aster31. Thank you for reporting!

There are two unrelated issues here:

1. MDTraj and MDAnalysis expect the box matrix to be raw-major, while molar expect is to be column-major. If you transpose your matrix before feeding it to molar you'd get an expected answer. I've added comments to the docs to make this distinction clear.
2. There was a bug in treatment of triclinic boxes, but it's not captured by your test :) I've fixed it.

Everything should be working correctly now in pymolar >=1.3.3. Please give it a try and let me know if it works.

Btw, your code contains unnecessary extra dimension for arrays, which makes it a bit hard to comprehend.