Logical error in PotentialTersoffGPU causes mismatch with CPU (revcross-specific) #2222
Description
A graduate student working with me was doing some simulations that used the RevCross potential, and was finding that hoomd-blue's GPU and CPU implementations give very different results. It seems that this is a known issue in the community that uses this potential (see, for instance, Section 3.3 of this thesis).
It's been a long time (since v2.0!) since I compiled hoomd from scratch, and I am hesitant to submit a PR without being able to fully validate the build pipeline (or be 100% confident in the fix). However, looking through the source code, I think it's clear where the problem is.
In particular: looking at the PotentialTersoffGPU.cuh code, I believe the issue is inside the if (evaluator::flag_for_RevCross) code path, around line 337 of the current trunk. The neighbor iteration loop contains this condition:
// I continue only if k is not the same as j
if ((cur_k > cur_j) && (cur_j > idx))In this check idx is the central particle, and the cur_j > idx condition filters out any triplet interaction where j has a lower tag/index than the central particle. But since the three-body potentials are calculated per central atom, excluding neighbors based on their index relative to i is physically incorrect. This check does not exist in the CPU implementation, which correctly iterates over all neighbors. Simply removing the (cur_j > idx) conditional should resolve the discrepancy.
Also, it's not entirely clear to me that the cur_k / cur_j logic for doing the neighbor searching is correct. I'm not exactly sure how this interacts with the strided tag loading in j loop, etc. At the very least, it seems to me that replacing the tag-based check with a more straightforward neighbor-list index check (i.e., make the conditional around line 337 neigh_idy > neigh_idx) to ensure each triplet is only counted once might be robust?