#2345 not fixed in NSGrid

[zemanj]

Expected behavior

The grid-based distance search in NSGrid works for an fcc crystal in a triclinic box.

Actual behavior

When calculating coordination numbers for all atoms in the correctly reproduced fcc-system from #2345, the NSGrid-based result is wrong. The test testsuite/MDAnalysisTests/lib/test_nsgrid.py::test_issue_2670 passes only because it loads the system from a pdb file, which lacks the coordinate precision required to reproduce the issue.

Code to reproduce the behavior

from collections import defaultdict, Counter
from numpy.testing import assert_allclose
import numpy as np
import MDAnalysis as mda
from MDAnalysisTests.datafiles import SURFACE_PDB, SURFACE_TRR
from MDAnalysis.lib.distances import distance_array, apply_PBC
from MDAnalysis.lib.mdamath import triclinic_box, triclinic_vectors
from MDAnalysis.lib.nsgrid import FastNS
from ase.build import fcc111

#search radius:
cutoff = 2.9

print("selected cutoff: {}".format(cutoff))

#generate triclinic fcc slab system with ASE:
ase_system = fcc111('Ag', size=(5, 5, 4), vacuum=10.0)
ase_box = triclinic_box(*ase_system.cell)
ase_positions = ase_system.positions.astype(np.float32)

#load similar system from corresponding PDB file:
u = mda.Universe(SURFACE_PDB)
pdb_box = u.dimensions
pdb_positions = u.atoms.positions


# Assert that systems are almost equal. PDB rounds to 3 decimal places so that
# the difference must not be larger than 0.0005:
assert_allclose(ase_box, pdb_box, atol=0.0005)
assert_allclose(ase_positions, pdb_positions, atol=0.0005)

#print boxes:
print()
print("ASE box: [{}]"
      "".format(" ".join(["{:15.10f}"] * len(ase_box))).format(*ase_box))
print("PDB box: [{}]"
      "".format(" ".join(["{:15.10f}"] * len(pdb_box))).format(*pdb_box))

#ASE bruteforce distance search:
ase_da = distance_array(ase_positions, ase_positions, box=ase_box)
ase_bf_cn = []
for i in range(ase_da.shape[0]):
    ase_bf_cn.append(len([ix for ix in range(ase_da.shape[1])
                 if 0.0 < ase_da[i, ix] < cutoff]))
bf_cns, bf_cn_counts = np.unique(ase_bf_cn, return_counts=True)
ase_bf_cn = dict(list(zip(list(bf_cns), list(bf_cn_counts))))

#ASE grid-based distance search:
ase_fns = FastNS(cutoff=cutoff, coords=ase_positions, box=ase_box, pbc=True)
ase_fns_results = ase_fns.self_search()

ase_nsgrid_pairs = ase_fns_results.get_pairs()
ase_nsgrid_indices = defaultdict(list)
for i, j in ase_nsgrid_pairs:
    ase_nsgrid_indices[i].append(j)
    ase_nsgrid_indices[j].append(i)
ase_nsgrid_cn = dict(Counter(len(v) for v in ase_nsgrid_indices.values()))

#PDB bruteforce distance search:
pdb_da = distance_array(pdb_positions, pdb_positions, box=pdb_box)
pdb_bf_cn = []
for i in range(pdb_da.shape[0]):
    pdb_bf_cn.append(len([ix for ix in range(pdb_da.shape[1])
                 if 0.0 < pdb_da[i, ix] < cutoff]))
bf_cns, bf_cn_counts = np.unique(pdb_bf_cn, return_counts=True)
pdb_bf_cn = dict(list(zip(list(bf_cns), list(bf_cn_counts))))

#PDB grid-based distance search:
pdb_fns = FastNS(cutoff=cutoff, coords=pdb_positions, box=pdb_box, pbc=True)
pdb_fns_results = pdb_fns.self_search()

pdb_nsgrid_pairs = pdb_fns_results.get_pairs()
pdb_nsgrid_indices = defaultdict(list)
for i, j in pdb_nsgrid_pairs:
    pdb_nsgrid_indices[i].append(j)
    pdb_nsgrid_indices[j].append(i)
pdb_nsgrid_cn = dict(Counter(len(v) for v in pdb_nsgrid_indices.values()))

# print results:
print()
print("ASE brute force coordination numbers: {}".format(ase_bf_cn))
print("PDB brute force coordination numbers: {}".format(pdb_bf_cn))
print()
print("ASE NSGrid coordination numbers:      {}".format(ase_nsgrid_cn))
print("PDB NSGrid coordination numbers:      {}".format(pdb_nsgrid_cn))

Output:

selected cutoff: 2.9

ASE box: [  14.4603338242   14.4603338242   27.0840873718   90.0000000000   90.0000000000   60.0000000000]
PDB box: [  14.4600000381   14.4600000381   27.0839996338   90.0000000000   90.0000000000   60.0000000000]

Brute force results:
ASE coordination numbers: {9: 50, 12: 50}
PDB coordination numbers: {9: 50, 12: 50}

NSGrid results:
ASE coordination numbers: {9: 49, 12: 47, 11: 3, 6: 1}     <== Fails for exactly reproduced system!
PDB coordination numbers: {9: 50, 12: 50}

Current version of MDAnalysis

• Which version are you using? 1.0.2-dev0
• Which version of Python? Python 3.8.5
• Which operating system? Ubuntu 20.04.2 LTS