Skip to content

Fix: Incorrect capped_distance in some triclinic box situations - #5481

Merged
BradyAJohnston merged 13 commits into
MDAnalysis:developfrom
BradyAJohnston:fix-4906-triclinic-pbc
Oct 19, 2026
Merged

BradyAJohnston merged 13 commits into
MDAnalysis:developfrom
BradyAJohnston:fix-4906-triclinic-pbc

Conversation

@BradyAJohnston

Copy link
Copy Markdown
Member

Fixes #4906

In triclinic boxes where both b_x and c_y are non-zero (e.g. truncated octahedra), the nsgrid and pkdtree methods of capped_distance / self_capped_distance miss pairs within the cutoff. This affects bond guessing, selections, etc. Reproduced with the issue's struc.txt and the truncated-octahedron script; both now match bruteforce.

Changes made in this Pull Request:

  • _triclinic_pbc: fix the a-axis wrapping bound, which could leave coordinates just outside the primary image (affects apply_PBC, pkdtree and nsgrid)
  • FastNS.coord2cellxyz: fix the same missing term in the x cell index (as suspected here)
  • FastNS._prepare_box: size grid cells from the box's perpendicular widths, so cells are never thinner than the cutoff
  • Add regression tests for apply_PBC, capped_distance, self_capped_distance and FastNS in skewed boxes

LLM / AI generated code disclosure

LLMs or other AI-powered tools (beyond simple IDE use cases) were used in this contribution: no

I did use Claude to scan through and track down the issue, but I wrote the code and the tests.

PR Checklist

  • Issue raised/referenced?
  • Tests updated/added?
  • Documentation updated/added?
  • package/CHANGELOG file updated?
  • Is your name in package/AUTHORS? (If it is not, add it!)
  • LLM/AI disclosure was updated.

Developers Certificate of Origin

I certify that I can submit this code contribution as described in the Developer Certificate of Origin, under the MDAnalysis LICENSE.

@BradyAJohnston BradyAJohnston changed the title Fix: #4096: Incorrect capped_distance in some triclinic box situations Fix: Incorrect capped_distance in some triclinic box situations Oct 9, 2026
@read-the-docs-community

read-the-docs-community Bot commented Oct 9, 2026 •

Copy link
Copy Markdown

@BradyAJohnston

Copy link
Copy Markdown
Member Author

Seems one of the old snapshots is failing but I'm pretty sure they are just incorrectly snapshotted values, looking into it

@codecov

codecov Bot commented Oct 9, 2026 •

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 93.91%. Comparing base (4ba8fe7) to head (dfb5613).

Additional details and impacted files
@@           Coverage Diff            @@
##           develop    #5481   +/-   ##
========================================
  Coverage    93.91%   93.91%           
========================================
  Files          184      184           
  Lines        22644    22644           
  Branches      3219     3219           
========================================
+ Hits         21266    21267    +1     
  Misses         917      917           
+ Partials       461      460    -1     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@IAlibay

IAlibay commented Oct 9, 2026

Copy link
Copy Markdown
Member

This looks great, unfortunately I won't have time to review. I'm going to un-request my review here and pick on someone else ;)

@IAlibay
IAlibay requested review from tylerjereddy and removed request for IAlibay October 9, 2026 11:42
@orbeckst

Copy link
Copy Markdown
Member

@IAlibay can we get this fix into the next release? It looks important enough as it can affect scientific correctness across many different quantities.

@orbeckst

Copy link
Copy Markdown
Member

(Yes, I get that it needs to be reviewed...)

@IAlibay

IAlibay commented Oct 13, 2026

Copy link
Copy Markdown
Member

@IAlibay can we get this fix into the next release? It looks important enough as it can affect scientific correctness across many different quantities.

Yes, as long as it can get merged by Sunday.

@BradyAJohnston BradyAJohnston added this to the 2.11 milestone Oct 13, 2026
@orbeckst

Copy link
Copy Markdown
Member

At first glance this looks very good, but I do need a bit more time to dig into it.

@IAlibay

IAlibay commented Oct 14, 2026

Copy link
Copy Markdown
Member

At first glance this looks very good, but I do need a bit more time to dig into it.

If "by Sunday" isn't possible - could you propose an alternative deadline? I do agree that this one probably falls under "critical fixes", so we could delay the release a tiny bit.

@orbeckst

Copy link
Copy Markdown
Member

I'll try to keep the "by Sunday" deadline for this one. Primarily I have to convince myself that what the LLMs are telling me makes sense to me, too. It's been a while since I dug into PBC handling. I even looked up the original algorithm in Frenkel & Smit again – the classics, in FORTRAN!

@BradyAJohnston

Copy link
Copy Markdown
Member Author

I've just made some tweaks, there was a 2x slowdown in FastNS construction from the additional numpy calls, but that's now been reduced to the same speeds. Tweaked the new tests a bit to clean them up.

- extensive CHANGELOG entry
- versionchanged (because the fix is important)
- code docs/comment updates
@orbeckst
orbeckst force-pushed the fix-4906-triclinic-pbc branch from 8a50ba1 to a7b88e1 Compare October 19, 2026 08:29

@orbeckst orbeckst left a comment

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It took me a while to understand the fixes. Both make sense to me.

I added more docs (and additional test).

Overall I am happy with it and think it definitely needs to go into 2.11.0.

Comment thread package/MDAnalysis/lib/nsgrid.pyx Outdated
Comment thread package/MDAnalysis/lib/nsgrid.pyx
Comment thread package/CHANGELOG
* 2.11.0

Fixes
* Fixed PBC handling in triclinic boxes in `_triclinic_pbc()` (used by

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I added an extensive entry because the impact of the fix is potentially big as it can directly affect computed properties. I know that we normally keep this short but in this case I'd like to be pretty clear.

Comment thread package/CHANGELOG
- `capped_distance()` and `self_capped_distance()` with the `nsgrid` or
`pkdtree` methods (and hence bond guessing, distance-based selections)
could miss pairs within the cutoff; `nsgrid` was affected in triclinic
boxes with two or more non-90º angles (e.g. rhombic dodecahedron),

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note that the patch contains multiple fixes and in fact affected more box shapes than originally noted.

if (!msr) {
// translate remainder of crd to central cell along a-axis
lbound = crd[1] * a_ax_yfactor + crd[2] * a_ax_zfactor;
lbound = -(crd[1] * bi3 + crd[2] * bi6) * box[0];

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes, this fix makes sense — I looked at it for a long time with diagrams, algebra, and LLM help, and I concur that this is necessary

Comment on lines +346 to +363
# Cells must be at least `cutoff` thick perpendicular to each pair of
# box faces, so size the grid by the box's perpendicular widths.
# The box matrix is lower triangular, so width = volume / face area
# has a closed form.
ax = self.triclinic_dimensions[XX]
bx = self.triclinic_dimensions[XY]
by = self.triclinic_dimensions[YY]
cx = self.triclinic_dimensions[XZ]
cy = self.triclinic_dimensions[YZ]
cz = self.triclinic_dimensions[ZZ]
volume = ax * by * cz
widths[0] = widths[1] = widths[2] = 0.0
# invalid boxes give an all-zero matrix, leave widths at zero
if volume > 0:
widths[0] = volume / math.sqrt((by * cz) ** 2 + (bx * cz) ** 2
+ (bx * cy - by * cx) ** 2)
widths[1] = by * cz / math.sqrt(cy * cy + cz * cz)
widths[2] = cz

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This fixes issues with boxes that have more than one non-90º angle

xyz : numpy.ndarray int of shape (3,)
Cell coordinate for the given coordinate

Background

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it's useful to write out how the calculation is done, given that the second key bug fix is the missing s_b b_x term


eps = DielectricConstant(ag, make_whole=False).run()
assert_allclose(eps.results["eps_mean"], 721.711, rtol=1e-03)
# reference values depend on correct PBC handling

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I wanted to make clear that this test is actually sensitive to PBC (which is good); the previous comment in this PR only said that the ref value changed but this historical info gets confusing over time without context

Comment thread testsuite/MDAnalysisTests/lib/test_distances.py
Comment thread testsuite/MDAnalysisTests/lib/test_nsgrid.py Outdated
@orbeckst

Copy link
Copy Markdown
Member

I quickly ran the new tests against the develop branch (without the fix).

  • The nsgrid test_issue_4906 test does not fail for the rhombic dodecahedron (box2).
  • All distance tests fail for all three example boxes.
$ pytest -v test_nsgrid_issue4906.py test_distances_issue4906.py
...
==================================================================== short test summary info =====================================================================
FAILED test_nsgrid_issue4906.py::test_issue_4906[box0] - AssertionError: assert {(np.int64(0)...64(325)), ...} == {(np.int64(0)...64(325)), ...}
FAILED test_nsgrid_issue4906.py::test_issue_4906[box1] - AssertionError: assert {(np.int64(0)...64(325)), ...} == {(np.int64(0)...64(325)), ...}
FAILED test_distances_issue4906.py::TestCappedDistances::test_capped_distance_skewed_triclinic[pkdtree-box0] - AssertionError: assert {(np.int64(0)...64(837)), ...} == {(np.int64(0)...64(837)), ...}
FAILED test_distances_issue4906.py::TestCappedDistances::test_capped_distance_skewed_triclinic[pkdtree-box1] - AssertionError: assert {(np.int64(0)...64(720)), ...} == {(np.int64(0)...64(720)), ...}
FAILED test_distances_issue4906.py::TestCappedDistances::test_capped_distance_skewed_triclinic[pkdtree-box2] - AssertionError: assert {(np.int64(0)...64(733)), ...} == {(np.int64(0)...64(733)), ...}
FAILED test_distances_issue4906.py::TestCappedDistances::test_capped_distance_skewed_triclinic[nsgrid-box0] - AssertionError: assert {(np.int64(0)...64(837)), ...} == {(np.int64(0)...64(837)), ...}
FAILED test_distances_issue4906.py::TestCappedDistances::test_capped_distance_skewed_triclinic[nsgrid-box1] - AssertionError: assert {(np.int64(0)...64(720)), ...} == {(np.int64(0)...64(720)), ...}
FAILED test_distances_issue4906.py::TestCappedDistances::test_capped_distance_skewed_triclinic[nsgrid-box2] - AssertionError: assert {(np.int64(0)...64(733)), ...} == {(np.int64(0)...64(733)), ...}
FAILED test_distances_issue4906.py::TestCappedDistances::test_self_capped_distance_skewed_triclinic[nsgrid-box0] - AssertionError: assert {(np.int64(0)...64(204)), ...} == {(np.int64(0)...64(204)), ...}
FAILED test_distances_issue4906.py::TestCappedDistances::test_self_capped_distance_skewed_triclinic[nsgrid-box1] - AssertionError: assert {(np.int64(0)...64(204)), ...} == {(np.int64(0)...64(204)), ...}
FAILED test_distances_issue4906.py::TestCappedDistances::test_self_capped_distance_skewed_triclinic[nsgrid-box2] - AssertionError: assert {(np.int64(0)...64(449)), ...} == {(np.int64(0)...64(449)), ...}
FAILED test_distances_issue4906.py::Test_apply_PBC::test_coords_in_central_image_skewed_tric[box0-serial] - assert np.False_
FAILED test_distances_issue4906.py::Test_apply_PBC::test_coords_in_central_image_skewed_tric[box0-openmp] - assert np.False_
FAILED test_distances_issue4906.py::Test_apply_PBC::test_coords_in_central_image_skewed_tric[box1-serial] - assert np.False_
FAILED test_distances_issue4906.py::Test_apply_PBC::test_coords_in_central_image_skewed_tric[box1-openmp] - assert np.False_
FAILED test_distances_issue4906.py::Test_apply_PBC::test_coords_in_central_image_skewed_tric[box2-serial] - assert np.False_
FAILED test_distances_issue4906.py::Test_apply_PBC::test_coords_in_central_image_skewed_tric[box2-openmp] - assert np.False_
=================================================== 17 failed, 5359 passed, 1 skipped, 1248 warnings in 14.34s ===================================================

@BradyAJohnston

Copy link
Copy Markdown
Member Author

Thanks for the thorough review and for updating the docs @orbeckt . I've just updated that test to definitely fail on develop and now pass. The bug was only showing up for certain cutoff values so I had managed to pick a combo which just slipped under on develop initially. Should have double checked it failed first though.

Should be all good to go.

@orbeckst

Copy link
Copy Markdown
Member

You beat me to it, I was just about to push a similarly updated test.

I have one that is guaranteed to fail with specially engineered points so I'll add it.

@BradyAJohnston

Copy link
Copy Markdown
Member Author

Always good to have more tests!

@BradyAJohnston
BradyAJohnston merged commit 6cd1f7e into MDAnalysis:develop Oct 19, 2026
25 checks passed
@BradyAJohnston
BradyAJohnston deleted the fix-4906-triclinic-pbc branch October 19, 2026 11:04
@orbeckst

Copy link
Copy Markdown
Member

🎉

@IAlibay I hope this is still acceptable.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Bugs for distance calculation functions

3 participants