I’m seeing a clean, exactly-reproducible factor-of-2 discrepancy between the energy pair_style lj/cut/dipole/cut reports and the textbook dipole-dipole interaction formula, for the simplest possible test case. I’ve ruled out every explanation I can check from outside the source, and I’m hoping someone familiar with the DIPOLE package’s internals can point out what I’m missing.
Setup: 2 particles, LJ epsilon set to 0 (isolated dipole-dipole test), aligned head-to-tail along their separation vector.
- Particle 1 at (-1.5, 0, 0), Particle 2 at (1.5, 0, 0) – separation r = 3.0
- Both dipoles set to (1.0, 0, 0) via the set command
- pair_style lj/cut/dipole/cut 5.0, pair_coeff * * 0.0 1.0 5.0
- units lj, boundary f f f (no periodicity, no kspace involved)
- run 0 (single-point energy evaluation)
Expected energy, from the textbook dipole-dipole formula (which matches LAMMPS’s own documented formula on the pair_style lj/cut/dipole/cut page exactly):
E = (mu_i . mu_j)/r^3 - 3*(mu_i . r)(mu_j . r)/r^5 = (1 - 3)/27 = -0.0740740741
Actual reported EPair: -0.037037037 – exactly half.
What I’ve confirmed rules out several possible explanations:
- Dumped the actual mux/muy/muz values LAMMPS stores after the set command – confirmed exactly (1.0, 0.0, 0.0) for both particles, no internal rescaling.
- Checked LAMMPS’s source (force->qqr2e in update.cpp) – confirmed qqr2e = 1.0 for lj units, and we never touched the dielectric command (default 1.0), so the C/epsilon prefactor mentioned in the pair_dipole docs is exactly 1.0.
- Checked the GPU-accelerated variant’s source (pair_lj_cut_dipole_cut_gpu.cpp) for the actual C++ energy expression – it computes ecoul += r3invpdotp - 3.0r5invpidotrpjdotr, then ecoul = factor_coulqqrd2e – this matches the textbook formula exactly with no missing factor, and factor_coul is 1.0 for two unrelated, unbonded atoms.
- Reproduced this exact discrepancy on two completely different machines/builds (a 29Aug2024 stable build linked against KISS FFT, and a fresh 2025.7.22 build linked against FFTW3) – identical -0.037037037 both times, so it’s not a build-specific fluke.
The one difference I haven’t been able to check directly: the GPU variant I found requests a full neighbor list and tallies via ev_tally_xyz_full, while my actual runs use the plain CPU pair_lj_cut_dipole_cut.cpp with a half neighbor list (confirmed via the printed “attributes: half, newton on”). I wasn’t able to get read access to the exact half-list tally call in that specific file to check whether something in that bookkeeping halves the contribution.
Given a half list is supposed to visit each unique pair exactly once (not twice like a full list), I wouldn’t expect any halving to be needed or applied – which is exactly why the result confuses me. Is there something specific to how lj/cut/dipole/cut’s compute() calls ev_tally for a half list that would produce this? Or is there a normalization convention for the dipole moment or the mu array’s 4th component (I noticed mu[i][3] gates the dipole-dipole branch in the GPU source) that differs from what I’m assuming?
Happy to share the full input script if useful. Using the 29Aug2024 stable release, DIPOLE and KSPACE packages built in, no MPI.