Fix stress normalization for partially periodic systems (#1509)#1523
Fix stress normalization for partially periodic systems (#1509)#1523darthjaja6 wants to merge 1 commit into
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.
9f27457 to
b1c210e
Compare
|
@cursor review |
Bugbot couldn't runBugbot is not enabled for your user on this team. Ask your team administrator to increase your team's hard limit for Bugbot seats or add you to the allowlist in the Cursor dashboard. |
|
@cursor review |
Bugbot couldn't runBugbot is not enabled for your user on this team. Ask your team administrator to increase your team's hard limit for Bugbot seats or add you to the allowlist in the Cursor dashboard. |
|
@cursor review |
Bugbot couldn't runBugbot is not enabled for your user on this team. Ask your team administrator to increase your team's hard limit for Bugbot seats or add you to the allowlist in the Cursor dashboard. |
|
@cursor review |
Bugbot couldn't runBugbot is not enabled for your user on this team. Ask your team administrator to increase your team's hard limit for Bugbot seats or add you to the allowlist in the Cursor dashboard. |
|
@cursor review |
Bugbot couldn't runBugbot is not enabled for your user on this team. Ask your team administrator to increase your team's hard limit for Bugbot seats or add you to the allowlist in the Cursor dashboard. |
|
@cursor review |
Bugbot couldn't runBugbot is not enabled for your user on this team. Ask your team administrator to increase your team's hard limit for Bugbot seats or add you to the allowlist in the Cursor dashboard. |
|
@cursor review |
Bugbot couldn't runBugbot is not enabled for your user on this team. Ask your team administrator to increase your team's hard limit for Bugbot seats or add you to the allowlist in the Cursor dashboard. |
|
@cursor review |
Bugbot couldn't runBugbot is not enabled for your user on this team. Ask your team administrator to increase your team's hard limit for Bugbot seats or add you to the allowlist in the Cursor dashboard. |
|
@cursor review |
![cursor[bot]](https://avatars.githubusercontent.com/in/1210556?s=60&v=4){kind=small}
There was a problem hiding this comment.
✅ Bugbot reviewed your changes and found no new issues!
Comment @cursor review or bugbot run to trigger another review on this PR
Reviewed by Cursor Bugbot for commit b1c210e. Configure here.
aacostadiaz
left a comment
aacostadiaz
left a comment
There was a problem hiding this comment.
Thanks for fixing this. I left a few comments but the fix itself looks correct to me.
| # stress divides by det(cell) later, so the fake volume shrinks it (#1509). | ||
| # keep it for fully aperiodic systems though (no stress there anyway) | ||
| if any(pbc): | ||
| cell = np.array(cell, dtype=float, copy=True) |
There was a problem hiding this comment.
Edge case worth considering: a zero row along a non-periodic direction now gives det(cell)=0, so stress comes out as zeros and training with stress in the loss silently NaNs all gradients. Hard to debug when it happens.
| cell = np.array(cell, dtype=float, copy=True) | |
| cell = np.array(cell, dtype=float, copy=True) | |
| # a zero row along a non-periodic direction would make det(cell)=0 | |
| # downstream (stress division, rcell) — keep the extended row there | |
| for dim in range(3): | |
| if not pbc[dim] and not cell[dim].any(): | |
| cell[dim] = extended_cell[dim] |
| edge_index, shifts, _, cell = get_neighborhood( | ||
| config.positions, cutoff=2.9, pbc=(True, True, False), cell=config.cell | ||
| ) | ||
| # should give back the real cell, not the blown-up one (#1509) | ||
| assert np.allclose(cell, config.cell) |
There was a problem hiding this comment.
I don't think this assert can catch the bug it references: the old code mutated the passed-in cell in place and returned that same array, so on the old code this compares the inflated cell with itself and passes. Snapshotting the cell before the call makes it a real regression guard.
| edge_index, shifts, _, cell = get_neighborhood( | |
| config.positions, cutoff=2.9, pbc=(True, True, False), cell=config.cell | |
| ) | |
| # should give back the real cell, not the blown-up one (#1509) | |
| assert np.allclose(cell, config.cell) | |
| cell_before = config.cell.copy() | |
| edge_index, shifts, _, cell = get_neighborhood( | |
| config.positions, cutoff=2.9, pbc=(True, True, False), cell=config.cell | |
| ) | |
| # should give back the real cell, not the blown-up one (#1509) | |
| assert np.allclose(cell, cell_before) |
|
Thanks for incorporating this fix and the regression test into #1533. Since the combined PR supersedes this one and also addresses the review comments here, I'll close this and follow the discussion there. |
4 participants
Summary
Fixes #1509. The analytic stress returned for partially periodic (slab) systems was silently scaled down by a large factor (~18x in the reported case), while energies and forces were unaffected.
Root cause
get_neighborhoodextends the cell along non-periodic directions so thatmatscipycan bin the atoms:This inflated cell was then returned from the function and stored in
AtomicData. Stress is normalized bydet(cell)incompute_forces_virials, so the artificial volume along the vacuum direction shrinks the stress by exactly the ratio of the inflated extent to the true extent. For the phosphorene example in the issue this ratio ismax_positions * 5 * cutoff / z_true ≈ 17.8, matching the reported17.81. Energies and forces don't touchdet(cell), which is why only stress was affected and the bug went unnoticed.Fix
Use the extended cell only for the neighbour-list call, and return the physical cell whenever at least one direction is periodic.
unit_shiftsare zero along non-periodic directions, so the shifts (and therefore the energy/forces) are identical. Fully non-periodic systems keep the extended cell, since stress is meaningless there and long-range models rely on a non-degenerate cell — so their behaviour is unchanged.Verification
Finite-difference check (analytic stress vs. strain derivative of the energy), model-independent:
Added a regression test (
test_stress_partial_pbc) that fails on the current code and passes with the fix, plus an assertion intest_half_periodicthat the physical cell is returned. Full existing test suite passes.Note
Medium Risk
Touches core neighbourhood/cell plumbing used by all models; behaviour change is scoped to partial PBC (fully periodic and fully aperiodic paths are preserved) but stress and any volume-derived quantities now differ for slabs.
Overview
Fixes incorrect analytic stress for slabs and other partially periodic systems, where stress was scaled down by the ratio of the artificially inflated vacuum cell to the true cell.
get_neighborhoodstill builds an extended cell along non-periodic axes formatscipyneighbour binning, but passes that cell only intoneighbour_list. When any direction is periodic, it now returns a copy of the physical cell (not the blown-up one) so downstream stress normalization viadet(cell)uses the real volume. Fully aperiodic systems still return the extended cell, unchanged.Tests add a physical-cell assertion in
test_half_periodicand atest_stress_partial_pbcregression that compares finite-difference strain energy to analyticsigma_yyon a Si TTF slab.Reviewed by Cursor Bugbot for commit b1c210e. Bugbot is set up for automated code reviews on this repo. Configure here.