Add support for charge-dipole interaction with PPPM#5059
Conversation
Charge-charge and charge–dipole interactions in
|
| channel | array | p | numerator | denominator |
|---|---|---|---|---|
| charge–charge | greensfn_qq |
1 | (k.k_n)^1 |
` |
| charge–dipole | greensfn_qmu |
2 | (k.k_n)^2 |
` |
| dipole–dipole | greensfn |
3 | (k.k_n)^3 |
` |
Each dipole in the interaction adds one power of (k.k_n) in the numerator and one factor of 1/|k|^2 in the denominator.
Consistency checks
p = 1reproducesPPPM::compute_gf_ikfrompppm.cppexactly (numerator = 4*pi/|k|^2,sum += (k.k_n)/|k_n|^2 * ...).- Continuum limit (
W -> 1, no aliasing,k_n -> k,gf_denom -> 1): all three collapse to the common Ewald kernel4*pi*exp(-k^2/4a^2)/k^2.
All three are filled in the single sweep in compute_gf_dipole().
Field combination (the subtle part)
Total reciprocal charge density of charges + dipoles:
rho~(k) = Q(k) - i*D(k), Q = FFT(charge density), D = k . P, P = FFT(dipole density)
The -i (from rho_bound = -div P) only matters for the cross term; energy and forces are derived from this single convention so they stay mutually consistent regardless of the FFT sign convention.
Per-channel potentials felt by each target species:
phi_on_q = gqq*Q - i*gqm*D (charge target: qq + cross)
phi_on_mu = gqm*Q - i*gdd*D (dipole target: cross + dd)
Imaginary-unit handling matches the source type:
- monopole terms use the re/im swap (multiply by
i), as in basePPPM::poisson_ik(work2[n] = -fkx*work1[n+1],work2[n+1] = +fkx*work1[n]); - dipole terms use a straight multiply, because the dipole potential already carries an implicit
-i.
k-space brick definitions (real/imag parts written out in poisson_ik_dipole()):
charge E-field vd{x,y,z}_brick : u_a = i k_a gqq Q + k_a gqm D
dipole E-field u{x,y,z}_brick_dipole : u_a = k_a gdd D + i k_a gqm Q
dipole gradient vd{ab}_brick_dipole : g_ab = i k_a k_b gdd D - k_a k_b gqm Q
The cross channel uses greensfn_qmu symmetrically for charge->dipole and dipole->charge, preserving Newton's third law.
Energy and virial
Per-k energy density (accumulated with the same s2 / 0.5*volume prefactor the dipole code already used):
gqq*|Q|^2 + gdd*|D|^2 + 2*gqm*Im(Q* D)
Self-energy (added in compute()):
dipole: -musqsum * 2 g^3 / (3 sqrt(pi)) (existing)
charge: -g*qsqsum/sqrt(pi) - (pi/2)*qsum^2/(g^2 V) (new; mirrors base PPPM)
Global virial: the common-kernel vg part applies to all three channels (they share the same continuum kernel). The two channels containing D pick up structure-factor strain corrections:
dipole-dipole: +2 s2 gdd k_b Re(P_a* D)
charge-dipole: +2 s2 gqm k_b Im(Q* P_a)
File-by-file changes
src/KSPACE/pppm_dipole.h
- New members:
double *greensfn_qq, *greensfn_qmu;andFFT_SCALAR *work5;.
src/KSPACE/pppm_dipole.cpp
- Constructor: initialize the new pointers to
nullptr. init(): remove the charge ban;qsum_qsq(0)sets up the charge channel.compute():- gate per-atom energy with charges (per-atom virial already gated);
- update
qsum_qsq()on atom-count change; - early-return only if both
musqsum == 0andqsqsum == 0; - reverse-comm count 3 -> 4, forward-comm count 9 -> 12;
- add charge self-energy term.
allocate()/deallocate(): allocate/freegreensfn_qq,greensfn_qmu,work5, and reuse inheriteddensity_brick,density_fft,vd{x,y,z}_brickfor the charge channel;npergrid9 -> 12.compute_gf_dipole(): fill all three Green's functions in one sweep.make_rho_dipole(): also spreadq[i]ontodensity_brick.brick2fft_dipole(): also remapdensity_fft.poisson_ik_dipole(): full rewrite — 4 forward FFTs (Q + 3 P), combined 3-channel energy/virial, 12 backward FFTs (3 charge-field, 3 dipole-field, 6 dipole-gradient).fieldforce_ik_dipole(): interpolate charge E-field fromvd{x,y,z}_brickfor the force on charges; dipole force/torque now include cross terms via the augmented bricks.pack/unpack_{forward,reverse}_grid(): carry the scalar charge density
(reverse) and the 3 charge-field components (forward).timing_1d/timing_3d: 4 forward + 12 backward FFTs (return 16).memory_usage(): account for the new arrays.
Accuracy / grid-size estimate
The RMS force error estimate combines the charge and dipole channels in quadrature, so grid sizing, the g_ewald solve, and the reported accuracy all account for charges.
k-space (compute_df_kspace_dipole):
df_dipole = sqrt(qopt_dipole/N) * mu2 / (3 Lx Ly Lz_slab) (p=3 qopt)
df_charge = sqrt(qopt_charge/N) * q2 / (Lx Ly Lz_slab) (p=1 qopt, base compute_qopt)
df_kspace = sqrt(df_dipole^2 + df_charge^2)
real-space (newton_raphson_f, final_accuracy_dipole):
df_dipole = mu2/sqrt(V g^4 rc^9 N) * sqrt(13/6 Cc^2 + 2/15 Dc^2 - 13/15 Cc Dc) * exp(-(g rc)^2)
df_charge = 2 q2 exp(-(g rc)^2) / sqrt(N rc V)
df_rspace = sqrt(df_dipole^2 + df_charge^2)
These feed set_grid_global() (grid sizing loop), adjust_gewald() (via the newton_raphson_f override), and the printed estimated accuracy. The initial g_ewald guess (find_gewald_dipole) is still dipole-only but is refined by
adjust_gewald() against the combined error.
Validation
See validation/pppm-charge-dipole/ for runnable inputs. Four internal self-consistency tests pass to k-space accuracy:
| test | quantity | relative agreement |
|---|---|---|
| g_ewald invariance | energy & pressure | 3e-7 / 3e-6 |
charge channel vs production pppm |
energy, forces | 6e-8 / 3e-7 |
| finite-difference force = -dE/dx (full system) | force | 6e-7 |
| finite-difference torque = -dE/dtheta (cross field) | torque | 2e-7 |
The g_ewald-invariance pressure test exercises the global virial (including the cross strain terms); the FD force/torque tests confirm the charge-dipole cross force and field.
Cross-validation vs ewald/disp/dipole (independent reference)
ewald/disp/dipole (src/KSPACE/ewald_disp.cpp, paired with
lj/long/dipole/long) already handles the same three channels and is the
natural external reference. It expresses the cross term differently: rather than
three Green's functions it keeps separate charge and dipole structure factors
S_q = sum q e^{ik.r} and S_mu = sum (mu.k) e^{ik.r}, and forms the
charge-dipole energy as 2 Im(S_q* S_mu) (compute_energy() lines ~1037-1040,
with matching force/torque in compute_force() ~952-962). The two
formulations are equivalent in the continuum limit.
Comparing the two solvers on one identical 216-particle net-neutral
charge+dipole configuration (write_data then read_data into both runs;
short-range LJ via cut, only Coulomb+dipole sent to k-space) -- the van der
Waals energy is bit-identical, confirming the configs match, so only the
electrostatics differ. Each solver chooses its own g_ewald, so the real/k
split differs; only the total ecoul + elong is the meaningful comparison:
| accuracy | pppm/dipole total |
ewald/disp/dipole total |
rel. diff |
|---|---|---|---|
| 1e-6 | -0.06525641545 | -0.06525523691 | 1.8e-5 |
| 1e-7 | -0.06525645223 | -0.06525629796 | 2.4e-6 |
The gap shrinks ~8x when the accuracy target tightens 10x, i.e. it is
discretization error (two different gridding / k-cutoff schemes converging to
the same continuum answer), not a systematic cross-term discrepancy. Because
both solvers include the charge-dipole cross term, this is the strongest
external check short of an analytic Ewald sum.
Per-atom energy and virial
Per-atom eatom/vatom are now supported for the combined charge+dipole system (the previous error->all gates were removed). The decomposition mirrors the global energy/virial split: each grid point's contribution is
attributed half to the charge that feels a potential and half to the dipole that feels a field.
Per-atom energy (u brick interpolated, tallied as q*u + mu.E):
charge: phi_q = gqq Q - i gqm D (new scalar u_brick)
dipole: E_mu = gdd D + i gqm Q (reuses u{x,y,z}_brick_dipole)
compute() adds the charge self-energy per atom
(-g q_i^2/sqrt(pi) - (pi/2) q_i qsum/(g^2 V)) alongside the existing dipole
self term.
Per-atom virial (six components; charge tallied with scalar v{0..5}_brick,
dipole with vector v{0..5}{x,y,z}_brick_dipole):
charge brick: scaleinv [ vg_j (gqq Q - i gqm D) - 2 i gqm k_b P_a ]
dipole brick: scaleinv k_c [ vg_j (gdd D + i gqm Q) + 2 gdd k_b P_a ]
The common-kernel vg part carries each species' energy share; the dipole structure-factor strain (2 gdd k_b Re(P_a* D)) stays on the dipole bricks, and the charge-dipole strain (2 gqm k_b Im(Q* P_a)) is carried on the charge bricks (the only species whose Re(Q* .) tally can express it). Summed, the two reproduce the global virial exactly.
Bookkeeping: poisson_peratom_dipole() runs on evflag_atom (energy needs the charge potential brick too); the per-atom forward-comm count grows 18 -> 25 (18 dipole virial + 1 charge potential + 6 charge virial bricks), and
npergrid for allocate_peratom() matches.
Per-atom self-consistency (200 random charge+dipole atoms, run 0)
| quantity | global | sum of per-atom | agreement |
|---|---|---|---|
| kspace energy | elong |
reduce sum pe/atom kspace |
exact (serial & 4 MPI) |
| kspace stress (all 6) | pressure NULL kspace |
-reduce sum stress/atom / V |
exact (serial & 4 MPI) |
Pure-dipole (charges zeroed) still reproduces both exactly, re-enabling the per-atom virial path that the charge-dipole work had temporarily disabled.
Slab correction (charge + dipole)
slabcorr() was rewritten to the unified Yeh-Berkowitz form for a system that carries both point charges and point dipoles. The total z dipole moment of the cell combines the charge first moment and the intrinsic dipoles:
M_z = sum_i (q_i z_i + mu_iz)
E_slab = (2*pi/V) [ M_z^2 - qsum R2 - qsum^2 Lz^2/12 ], R2 = sum_i q_i z_i^2
This replaces the previous dipole-only (2*pi/V)(sum mu_iz)^2/12 self term, whose 1/12 was inconsistent with the code's own -4*pi/V torque/field (the torque already implied the no-1/12 form). The corrected energy, per-atom
energy, charge force, and dipole torque now all derive from the same E_slab:
f_iz = -4*pi/V q_i (M_z - qsum z_i) (charges only)
tau_ix = -4*pi/V M_z mu_iy, tau_iy = +4*pi/V M_z mu_ix (dipoles)
eatom_i = 2*pi/V [ m_i M_z - q_i(0.5(R2 + qsum z_i^2) + qsum Lz^2/12) ]
with m_i = q_i z_i + mu_iz. The per-atom energy reduces exactly to the base PPPM::slabcorr charge form when mu = 0 and sums to E_slab.
Why the 1/12 moved (the e_slabcorr change in detail)
Two distinct things changed on the energy line. Before:
// dipole_all = sum mu_iz (dipoles only)
e_slabcorr = MY_2PI*(dipole_all*dipole_all/12.0)/volume;
// \_____ M^2 term ____/ /12 (no charge terms at all)After:
// dipole_all = sum (q_i z_i + mu_iz) = M_z (combined moment)
e_slabcorr = MY_2PI*(dipole_all*dipole_all
- qsum*dipole_r2
- qsum*qsum*zprd_slab*zprd_slab/12.0)/volume;
// \__ M^2 __/ \_ q.R2 _/ \____ qsum^2 Lz^2/12 ____/So (a) the M^2 term lost its 1/12, and (b) the 1/12 now sits only on the (qsum*Lz)^2 term. Both are deliberate:
-
The
qsum^2 Lz^2/12term is where1/12actually belongs. It is the second moment of a uniformly charged slab -- integratingz^2across the slab thickness yields theLz^2/12. It is purely a property of the net chargeqsumsmeared over the box (the non-neutral background correction of J. Chem. Phys. 131, 094107) and has nothing to do with the dipole moment. -
The
M_z^2self term must have no1/12. The Yeh-Berkowitz slab energy is2*pi M_z^2/V. This is visible from the point-dipole-as-charge-pair limit: a+/-qpair separated bydinzcontributesM_z = q d = mu_z, and the standard charge slab energy2*pi M_z^2/Vgives2*pi mu_z^2/V-- no1/12. The originalmu^2/12was therefore a factor-of-12 error on this term.
The error was self-evident from an internal inconsistency: the old energy carried 1/12 but the code's torque used ffact = -4*pi/V (no 1/12). Deriving the torque from the old energy, tau = -dE/dtheta, gives a value 12x smaller
than what the code applied -- i.e. energy and torque disagreed by exactly the stray factor. The corrected 2*pi M_z^2/V energy is consistent with that -4*pi/V torque, and the finite-difference test below confirms
tau = -dE/dtheta to ~1e-10. For a neutral system (qsum = 0) the q.R2 and qsum^2 Lz^2/12 terms vanish, leaving E = 2*pi M_z^2/V -- the term whose value changed for existing pure-dipole slab runs.
Validation (slab geometry, kspace_modify slab 3.0):
| test | quantity | agreement |
|---|---|---|
| dipole torque = -dE/dtheta (rotate a dipole) | torque | 2e-10 |
| charge force = -dE/dz (displace a charge) | force | 3e-6 (grid-level) |
sum per-atom energy vs global elong (non-neutral q+mu) |
energy | exact (serial & MPI) |
Note: ewald/dipole (src/KSPACE/ewald_dipole.cpp) still carries the old mu^2/12 slab form and the per-atom/non-neutral block; the same correction would apply there but was left untouched (out of scope for the PPPM work).
Slab correction and the stress tensor
slabcorr() writes into the energy and the forces/torques, but not into the stress: it touches energy (global), eatom[i] (per-atom energy), f[i][2] (charge z-forces), and torque[i][...] (dipole torques), and never writes
virial[] (global) or vatom[][] (per-atom). So the slab term is absent from the stress at both the global and per-atom levels. This matches base PPPM::slabcorr, which is likewise energy-and-force only.
In principle the correction does have a virial: E_slab = (2*pi/V)[M_z^2 - qsum R2 - qsum^2 Lz^2/12] is explicitly
geometry-dependent (M_z carries a length, R2 a length^2, Lz a box dimension, plus the 1/V prefactor), so -dE_slab/d(strain) is a nonzero stress. LAMMPS simply does not evaluate it -- a long-standing choice in the base
charge code, not something introduced here. Two reasons it is acceptable:
- Use case. The slab correction only applies to slab geometry (
boundary p p fplus a vacuum gap), and those runs keep the box fixed in z (no barostat into the vacuum), so the slab contribution toPzzis never
used. - Self-consistency. Because the slab term is excluded from both the global virial and the per-atom stress, the two still agree with each other. That is exactly why the "sum per-atom kspace stress == global kspace stress" check still holds exactly under slab -- both sides omit the same term.
Practical consequence: pressure / stress/atom in a slab run exclude the (usually small) slab term, identically to every charge-only pppm + slab run. A fully stress-complete slab correction would mean deriving -dE_slab/d(strain) and adding it to both virial[] and vatom[][] (in base PPPM and here together) -- a separate enhancement.
Known limitations / TODO
slabcorradds no virial term -- neither global nor per-atom -- matching basePPPM::slabcorr(energy and force only). See "Slab correction and the stress tensor" above for why.- Initial
g_ewaldguess infind_gewald_dipoleis dipole-only (refined afterward); a combined initial guess is a minor future refinement. - No KOKKOS / GPU variant yet.
There was a problem hiding this comment.
Pull request overview
This PR extends the PPPMDipole kspace solver to support mixed systems containing both point charges and point dipoles, adding the missing charge–charge and charge–dipole (cross) channels alongside the existing dipole–dipole channel.
Changes:
- Adds separate Green’s functions for charge–charge and charge–dipole interactions, and expands FFT-grid dataflow to compute the additional fields and per-atom tallies.
- Removes the prior runtime restriction that prevented using
pppm/dipolewhen charges are present, and updates energy/virial/self-term accounting accordingly. - Revises slab correction logic to use a combined cell dipole moment including both
q*zandmu_z, and adds corresponding force/torque corrections.
Reviewed changes
Copilot reviewed 2 out of 2 changed files in this pull request and generated 4 comments.
| File | Description |
|---|---|
| src/KSPACE/pppm_dipole.h | Adds new work buffer and separate Green’s-function arrays for the additional interaction channels. |
| src/KSPACE/pppm_dipole.cpp | Implements mixed charge/dipole density spreading, k-space solve, field interpolation, per-atom tallies, and updated slab correction; updates comm sizes and FFT timing counts. |
| fft1->compute(work1,work1,FFT3d::FFT3d::FORWARD); | ||
| fft1->compute(work1,work1,FFT3d::FFT3d::FORWARD); | ||
| fft1->compute(work1,work1,FFT3d::FFT3d::FORWARD); | ||
| fft1->compute(work1,work1,FFT3d::FFT3d::FORWARD); | ||
| fft2->compute(work1,work1,FFT3d::FFT3d::BACKWARD); |
There was a problem hiding this comment.
fixed these errors, which are present in the original version.
| nxlo_out,nxhi_out,"pppm_dipole:vdz_brick"); | ||
| memory->create(density_fft,nfft_both,"pppm_dipole:density_fft"); | ||
|
|
||
| memory->create(densityx_fft_dipole,nfft_both,"pppm_dipole:densityy_fft_dipole"); |
| dipoleflag = atom->mu?1:0; | ||
| qsum_qsq(0); // q[i] might not be declared ? | ||
|
|
||
| if (dipoleflag && (q2 != 0.0)) | ||
| error->all(FLERR,"Cannot (yet) use charges with Kspace style PPPMDipole"); | ||
| qsum_qsq(0); // computes qsum, qsqsum, q2 for the charge channel | ||
|
|
There was a problem hiding this comment.
updated the doc page
| dipoleflag = atom->mu?1:0; | ||
| qsum_qsq(0); // q[i] might not be declared ? | ||
|
|
||
| if (dipoleflag && (q2 != 0.0)) | ||
| error->all(FLERR,"Cannot (yet) use charges with Kspace style PPPMDipole"); | ||
| qsum_qsq(0); // computes qsum, qsqsum, q2 for the charge channel | ||
|
|
|
@ndtrung81 nice work, I was just thinking about the disabled per-atom virial in |
|
@stanmoore1 thanks! You are correct, |
|
Claude reviewing Claude: Review: pppm/dipole charge–dipole support (branch
|
…3. Finding 1 and C4 will be addressed in a separate PR on ewald/disp slabcorr()
The dipole (spin) self-energy term in slabcorr() carried a 1/12 factor on M_z^2, but the conjugate torque/force term uses -4*pi/V with no such factor. This made the energy and the torque/force disagree by exactly 1/12, breaking torque = -dE/dtheta self-consistency. Remove the erroneous 1/12 so the self term is E = (2pi/V) M_z^2, which is consistent with the -4*pi/V field acting on the dipoles/spins. This is the same fix applied upstream to PPPMDipole (lammps#5059), now carried over to the Ewald dipole variants. The matching per-atom energy prefactor is corrected as well. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01S9FMNhdySkR3iAXXtdjbVN
The dipole (spin) self-energy term in slabcorr() carried a 1/12 factor on M_z^2, but the conjugate torque/force term uses -4*pi/V with no such factor. This made the energy and the torque/force disagree by exactly 1/12, breaking torque = -dE/dtheta self-consistency. Remove the erroneous 1/12 so the self term is E = (2pi/V) M_z^2, which is consistent with the -4*pi/V field acting on the dipoles/spins. This is the same fix applied upstream to PPPMDipole (lammps#5059), now carried over to the Ewald dipole variants. The matching per-atom energy prefactor is corrected as well. Co-Authored-By: Stan Moore <stanmoore1@gmail.com> Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01S9FMNhdySkR3iAXXtdjbVN
The dipole (spin) self-energy term in slabcorr() carried a 1/12 factor on M_z^2, but the conjugate torque/force term uses -4*pi/V with no such factor. This made the energy and the torque/force disagree by exactly 1/12, breaking torque = -dE/dtheta self-consistency. Remove the erroneous 1/12 so the self term is E = (2pi/V) M_z^2, which is consistent with the -4*pi/V field acting on the dipoles/spins. This is the same fix applied upstream to PPPMDipole (lammps#5059), now carried over to the Ewald dipole variants. The matching per-atom energy prefactor is corrected as well. Co-Authored-By: Stan Moore <stanmoore1@gmail.com> Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01S9FMNhdySkR3iAXXtdjbVN
The dipole (spin) self-energy term in slabcorr() carried a 1/12 factor on M_z^2, but the conjugate torque/force term uses -4*pi/V with no such factor. This made the energy and the torque/force disagree by exactly 1/12, breaking torque = -dE/dtheta self-consistency. Remove the erroneous 1/12 so the self term is E = (2pi/V) M_z^2, which is consistent with the -4*pi/V field acting on the dipoles/spins. This is the same fix applied upstream to PPPMDipole (lammps#5059), now carried over to EwaldDipole, EwaldDipoleSpin and PPPMDipoleSpin. The matching per-atom energy prefactor is corrected as well. Co-Authored-By: Stan Moore <stanmoore1@gmail.com> Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01S9FMNhdySkR3iAXXtdjbVN
The dipole (spin) self-energy term in slabcorr() carried a 1/12 factor on M_z^2, but the conjugate torque/force term uses -4*pi/V with no such factor. This made the energy and the torque/force disagree by exactly 1/12, breaking torque = -dE/dtheta self-consistency. Remove the erroneous 1/12 so the self term is E = (2pi/V) M_z^2, which is consistent with the -4*pi/V field acting on the dipoles/spins. This is the same fix applied upstream to PPPMDipole (lammps#5059), now carried over to EwaldDipole, EwaldDipoleSpin and PPPMDipoleSpin. The matching per-atom energy prefactor is corrected as well. Co-Authored-By: Stan Moore <stanmoore1@gmail.com> Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01S9FMNhdySkR3iAXXtdjbVN
|
I fixed F1 in #5053. |
stanmoore1
left a comment
stanmoore1
left a comment
There was a problem hiding this comment.
I approve. One nitpick: "charge_flag" may be more canonical than "has_charges", but either is probably fine.
github-project-automation
Bot
moved this from In Progress
to Done
in LAMMPS Pull Requests
Summary
This PR adds the charge-charge and charge-dipole interactions for the current
PPPMDipoleclass when there are also charges in the dipole system. The changes in this PR is necessary to allowpppm/dipoleto be consistent withlj/cut/dipole/long, which already has charge-charge and charge-dipole interactions in the real space.Add a new example input
examples/dipole/in.charge_dipoleand the reference log files.Related Issue(s)
Issue #4675
Author(s)
Trung Nguyen (U Chicago) and Claude Code Opus 4.8
Licensing
By submitting this pull request, I agree, that my contribution will be included in LAMMPS and redistributed under either the GNU General Public License version 2 (GPL v2) or the GNU Lesser General Public License version 2.1 (LGPL v2.1).
Artificial Intelligence (AI) Tools Usage
This pull request was prepared with the assistance of an AI coding tool (Claude Code, Anthropic, model Claude Opus 4.8) working under the direction, review, and testing of the author.
Backward Compatibility
Should maintain backward compatiblity with pure dipole system. There is a potential bug fix with
slabcorr()that needs confirmation through reivews.Implementation Notes
Adding the Greens functions for charge-charge (similar to PPPM) and charge-dipole parts, and field forces. Summary of the changes and justification in the conversation.
Post Submission Checklist
Further Information, Files, and Links