Unexpected neighborlist behaviors when (n_atoms - 1) % 32 == 0

[maxentile]

Recently noticed two small oddities when interacting with the neighborlist from Python. Cause currently unclear.

Possibly related to size-dependent test failures @mcwitt has recently observed.

---

1. On some inputs, the neighborlist uses block bounds where one bounding box has extent == (0,0,0).

Randomly sampling (conf, box) input pairs with random # atoms suffices to find a few examples of this behavior:

from timemachine.lib import custom_ops
import numpy as np

def num_zero_extent_blocks(conf, box):
    nblist = custom_ops.Neighborlist_f64(len(conf))
    bb_centers, bb_extents = nblist.compute_block_bounds(conf, box, 32)
    return sum((bb_extents == 0).all(1))

np.random.seed(0)

inputs = []
results = []

for _ in range(1000):
    n_atoms = np.random.randint(100, 3000)
    scale = np.random.rand() + 3
    
    conf = scale * np.random.randn(n_atoms, 3)
    box = scale * np.eye(3)
    
    n = num_zero_extent_blocks(conf, box)
    if n > 0:
        print(f'found an instance that produces zero-extent blocks!\n(# atoms = {n_atoms}, # zero-extent blocks = {n})')
        inputs.append((conf, box))
        results.append(n)

print(len(results))  # 42

Looking at the number of atoms in each of these, I noticed that they all have in common (n_atoms - 1) % 32 == 0.

---

2. Segfault

Modifying the loop to only try inputs where (n_atoms - 1) % 32 == 0 reliably triggers segfaults after a random number of iterations.

(Note: Segfaults are not triggered by this loop if n_atoms = 32 * x + 1 is changed to n_atoms = 32 * x + 2, for example.)

from timemachine.lib import custom_ops
import numpy as np

def num_zero_extent_blocks(conf, box):
    nblist = custom_ops.Neighborlist_f64(len(conf))
    bb_centers, bb_extents = nblist.compute_block_bounds(conf, box, 32)
    return sum((bb_extents == 0).all(1))

np.random.seed(0)

inputs = []
results = []

for t in range(1000):
    x = np.random.randint(10, 50)
    n_atoms = 32 * x + 1
    scale = np.random.rand() + 3
    
    conf = scale * np.random.randn(n_atoms, 3)
    box = scale * np.eye(3)

    n = num_zero_extent_blocks(conf, box)
    if n > 0:
        print(f'iteration {t}: found an instance that produces zero-extent blocks!\n(# atoms = {n_atoms}, # zero-extent blocks = {n})')
        inputs.append((conf, box))
        results.append(n)

print(len(results))

[maxentile]

(Note: the 2nd issue appears completely unrelated to the 1st -- see @mcwitt 's minimized reproduction + fix for (2) in #664 )

[maxentile]

Closing since behavior (1) is actually expected!