Code for KC local normal optimized DispFn doesn't look right?

[ExpHP]

It seems to me that, if atoms A1 and A2 are bonded in the same layer, and atom B is in another layer, then displacing A1 should modify at least all of the following potential terms:

• A1 <-> B
• A1 <-> A2
• A2 <-> B (because the normal changed)

It might even impact A2 <-> A3 and B <-> B2 (where A3 is bonded to A2, B is bonded to B2), as those terms originate from the computation of A2 <-> B.

Looking over the code I wrote, it seems to me that the first two of these terms are produced, but not the others? Then again, I could've sworn that I did write something to make the code compute interlayer terms for all neighbors (I just can't seem to find the part that takes care of this right now, and it has me worried...).

I should compare the optimized DispFn to the naive, O(n) helper DispFn (the one created by _default_initialize_disp_fn) on a large structure to see whether any terms really are missing. I'll be very disappointed in past me if they don't match...

[ExpHP]

Oops, false alarm. I totally forgot about this:

https://github.com/ExpHP/rsp2/blob/56249d1bc5d6ca93a709bc431ec2d4bd9546b54d/src/tasks/potential/homestyle.rs#L380-L383

That is to say, local normals don't actually use any optimization for the DispFn, instead doing things in the definitely-correct (but slow!) way. Oftentimes the cost of diagonalization totally swamps force constant calculation anyways since I switched to dense diagonalization, so I guess it hasn't been missed too sorely.

IIRC it was the intralayer terms that made me chicken out. Because those terms are affected by multiple interlayer pairs, this complicated any attempt to compute the difference between the equilibrium forces and the perturbed forces (because I have a full sum in the equilibrium forces, but I might only have computed a partial sum in the perturbed forces).