Allow cell vectors to have length zero if pbc is False

[jameskermode]

Consider a dimer:

from ase import Atoms
from ase.units import eV, Ang, GPa
from matscipy.neighbours import neighbour_list

d = 2.5*Ang
a = Atoms('2Cu', positions=[(0., 0., 0.), (0., 0., d)])
ij = neighbour_list('ij', a, 3.0)

The ase.Atoms constructor correctly sets pbc=False for the zero-length vectors.
The current implementation fails with a LinAlgError when trying to invert the cell.

Traceback (most recent call last):
  File "neigh_pbc_test.py", line 7, in <module>
    ij = neighbour_list('ij', a, 2.0)
  File "/Users/jameskermode/.julia/v0.5/Conda/deps/usr/lib/python2.7/site-packages/matscipy/neighbours.py", line 127, in neighbour_list
    np.linalg.inv(a.cell.T), a.pbc,
  File "/Users/jameskermode/.julia/v0.5/Conda/deps/usr/lib/python2.7/site-packages/numpy/linalg/linalg.py", line 513, in inv
    ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
  File "/Users/jameskermode/.julia/v0.5/Conda/deps/usr/lib/python2.7/site-packages/numpy/linalg/linalg.py", line 90, in _raise_linalgerror_singular
    raise LinAlgError("Singular matrix")
numpy.linalg.linalg.LinAlgError: Singular matrix

Since the inverse shouldn't be needed when pbc=False, I tried the following patch:

index 30096ba..c8567de 100644
--- a/matscipy/neighbours.py
+++ b/matscipy/neighbours.py
@@ -123,7 +123,12 @@ def neighbour_list(quantities, a, cutoff, *args):
                    np.bincount(i, weights=pair_forces[:, 2], minlength=len(a))
     """

+    if np.any(a.pbc):
+        invcell = np.linalg.inv(a.cell.T)
+    else:
+        invcell = np.eye(3)
+
     return _matscipy.neighbour_list(quantities, a.cell,
-                                    np.linalg.inv(a.cell.T), a.pbc,
+                                    invcell, a.pbc,
                                     a.positions, cutoff, *args)

However, we need to think about this a bit more, as this leads to zero cell volume, and would screw up the division of box into cells within C neighbour list code.

$ python neigh_pbc_test.py
Traceback (most recent call last):
  File "neigh_pbc_test.py", line 7, in <module>
    ij = neighbour_list('ij', a, 2.0)
  File "/Users/jameskermode/.julia/v0.5/Conda/deps/usr/lib/python2.7/site-packages/matscipy/neighbours.py", line 133, in neighbour_list
    a.positions, cutoff, *args)
RuntimeError: Zero cell volume.

[pastewka]

I think we should use some sort of shrink-wrapped cell in that case. The neighbor list needs a cell because it uses cell subdivision internally to look for neighbor. The shrink-wrapped cell should be fastest, because it minimizes the number of cells with no atoms in them. (The neighbor list will not necessarily loop over those, but it needs to initialize them and they may destroy cache coherency.)

[jameskermode]

Sounds good, but if the shrink-wrapped lengths are zero (as they would be in x and y here) we will still get zero volume and the subdivision will fail. Round up to at least cutoff?

[pastewka]

Yes, set to at least cutoff. Any small values should do but cutoff is probably reasonable.

[pastewka]

Part of this could involve changing inv to pinv. ASE uses pinv throughout to compute reciprocal cell.

[stenczelt]

+1 - I have been seeing the same issue with LJ clusters. I can manually set the cell to some shrink-wrapped version though.

[stenczelt]

may relate te usability in ACEsuit/mace#52 as well