Fix non-periodic cell in get_neighborhood: extent-based padding + physical cell on partial PBC#1533
Conversation
get_neighborhood extends the cell along non-periodic directions so that matscipy can bin atoms, but the extended cell leaked out of the function and was stored in AtomicData. Stress is normalized by det(cell), so for slabs the analytic stress was silently scaled down by the ratio of the artificial extent to the true vacuum extent (~18x in the reported case), while energies and forces were unaffected. Keep the extended cell internal to the neighbour-list call and return the physical cell whenever at least one direction is periodic. Shifts are unchanged because unit_shifts are zero along non-periodic directions. Fully non-periodic systems retain the previous behaviour. Add a finite-difference regression test comparing analytic stress with the strain derivative of the energy on a slab.
The combined ACEsuit#1473 + ACEsuit#1509 fix makes get_neighborhood return the physical cell for partially periodic systems. That cell also becomes AtomicData's volume/rcell, which drive PolarMACE's k-space electrostatics, so the change is not stress-only for electrostatic models. - test_polar_slab_partial_pbc_cell_contract: deterministic guard that a slab gets the physical cell (not the matscipy blow-up), that volume/rcell derive from it, that the short-range graph is invariant to the vacuum, and that the forward is finite on a partial-PBC electrostatic model (a path neither PR exercised). - test_polar_slab_electrostatics_converge_with_vacuum: physics check for the reviewer's concern that the vacuum affects the electrostatic part. On a dipole-free slab the 3D k-space sum must converge as the vacuum grows; non-convergence would signal the need for a 2D / slab-corrected sum.
… cell The extent-based non-periodic cell (from the neighborhood fix) has ~2 orders smaller volume than the old max(|pos|)*5*cutoff cell, so stress = virial/det(cell) scales up proportionally. The e3nn<->cueq parity is unchanged, but the absolute stress bounds (1e-8 float32, 1e-11 float64) were calibrated against the old inflated cell and now trip at ~2.3e-8. Loosen to 1e-6 / 1e-9; energy and force bounds are untouched (they don't divide by det(cell)).
Mutation testing showed two gaps in the added tests: - Nothing guarded the extent-based non-periodic cell (the ACEsuit#1473 OOM fix): reverting to max(|positions|)*5*cutoff kept every test green. Add test_nonperiodic_cell_is_extent_based, which asserts the fictitious cell is origin-independent (same molecule shifted far away gets the same cell) and tracks the atom extent + cutoff padding. Verified to fail on the old formula. - test_polar_slab_electrostatics_converge_with_vacuum was documented as if it exercised the Yeh-Berkowitz slab correction, but its slab is dipole-free and converges with the correction on or off (verified by toggling include_pbc_corrections). Reword the docstring to scope it honestly as a finite/non-diverging smoke check, not a guard of the correction.
Has anyone carefully thought about whether the implemented solution works when the cell vectors are not orthogonal? It's easy for naive heuristics to fail under these circumstances |
get_neighborhood replaces the non-periodic (vacuum) row of the cell with an axis-aligned vector for the matscipy binning. This verifies the resulting neighbour list is still exact on non-orthogonal slabs, where the periodic vectors keep their true skewed directions. Parametrized over hexagonal, monoclinic and triclinic in-plane lattices x two cutoffs x three random configurations (18 cases), compared edge-by-edge (indices + shift vectors) against a brute-force reference whose image range grows until the edge set converges (complete regardless of skew).
Good point. I added tests that check the neighbour list on skewed (non-orthogonal) cells: hexagonal, monoclinic and triclinic slabs, with different cutoffs and random atom positions. Each one is compared atom by atom against a brute-force reference and they all match, so the extent-based cell works fine here. |
| for dim in range(3): | ||
| if not pbc[dim]: | ||
| extent = positions[:, dim].max() - positions[:, dim].min() | ||
| extended_cell[dim, :] = (extent + 2 * cutoff + 1) * identity[dim, :] |
There was a problem hiding this comment.
I want to think about it more, but I don't understand how this can possibly work in general. It replaces periodic vector a_i with a vector along Cartesian direction i, with 0-x, 1-y, and 2-z, which doesn't necessarily have anything to do with a_i.
There was a problem hiding this comment.
For example, will this work with
[[0, 5, 0],
[5, 0, 0],
[0, 0, -5]]
and pbc=[False, True, True]? Won't it become
[[L, 0, 0],
[5, 0, 0],
[0, 0, -5]]
where L is the related to the x extent of the system?
If a system is non-periodic along a dimension, I think you need to make that vector along np.cross(a[i+1 % 3], a[i + 2 % 3]) with some sufficient magnitude, which is not entirely trivial to calculate. And the possibility of magnitude 0 vectors complicates this a bit more. You probably have to construct 1 or 2 (or 3) vectors that are independent, to use as dummy directions, before doing the cross product.
|
Can we propose that neighbourlist packages, such as matscipy take on this functionality, of determining the minimum padding? |
|
I agree that should be the case. |
6 participants
Summary
Combines and reconciles two open fixes to
get_neighborhood's handling of the fictitious cell used along non-periodic directions. They touch the same lines and conflict textually, so they are merged into one coherent change. Supersedes #1473 and #1523.max(abs(positions)) * 5 * cutoffcell depended on the absolute coordinate origin and produced enormous cells, OOM-ing PolarMACE's k-space electrostatics for ordinary molecules far from the origin.AtomicData, sostress = virial / det(cell)was scaled by a large factor (~18x) for partially periodic (slab) systems.What changed
get_neighborhoodnow:extent + 2*cutoff + 1) instead ofmax(abs(positions)) * 5 * cutoff— origin-independent, far smaller, identical neighbour lists (fixes the OOM).Effect on electrostatic models (PolarMACE)
The returned cell also becomes
AtomicData'svolume/rcell, which drive PolarMACE's electrostatics.graph_longrangealready applies the correct boundary treatment but they divide bydet(cell). With the old inflated cell that volume was fake, so the corrections were silently scaled away; returning the physical cell gives them theright volume. So this change is not only a stress fix, it also makes PolarMACE's slab/molecule electrostatic corrections numerically correct. (
graph_longrangeonly supportsTTT/TTFand raises otherwise, so no other partial-PBC case is affected.) Adequate vacuum along the non-periodic axis remains the caller's responsibility.Note
Medium Risk
Touches core graph construction used by all models; slab stress and PolarMACE electrostatics numerics change by design, though behaviour is heavily regression-tested.
Overview
get_neighborhoodno longer mutates the caller’s cell and fixes how fictitious non-periodic cell rows are built and what gets returned downstream.For matscipy neighbour binning only, non-periodic axes use an extent + cutoff padded cell (
extent + 2*cutoff + 1) instead ofmax(|positions|) * 5 * cutoff, so padding is origin-independent and much smaller—avoiding huge k-space boxes that could OOM PolarMACE electrostatics while keeping the same neighbour lists.When any axis is periodic, the function returns the physical cell (with a fallback extended row if a non-periodic axis is all zeros), so
AtomicDatavolume/rcell, slab stress (virial / det(cell)), and long-range dipole corrections use the real box volume. Fully aperiodic systems still return the extended cell.Tests cover partial-PBC slabs (PolarMACE cell contract and vacuum convergence), brute-force neighbour parity on skewed slabs, extent-based aperiodic cells, slab stress vs finite differences, and looser cueq stress parity bounds after the smaller cell volume.
Reviewed by Cursor Bugbot for commit 6cca0b2. Bugbot is set up for automated code reviews on this repo. Configure here.