Long range¶

Long range electrostatics split the Coulomb energy into a short range part computed in real space and a long range part computed in reciprocal space:

with \(S(\mathbf{k}) = \sum_i q_i e^{i\mathbf{k}\cdot\mathbf{r}_i}\), \(g\) the Ewald splitting parameter (g_ewald), \(V\) the cell volume and \(Q=\sum_i q_i\) the total charge. The last two terms are the self energy and the neutralizing background energy for a non-neutral system.

exaStamp provides two reciprocal space solvers:

Both methods have the following properties:

• Orthogonal and triclinic periodic cells are supported.
• Charges are read per particle (charge field) by default. To use the species charges instead, set per_atom_charge: false on coulombic_ewald_short_range.
• Per-atom energies and the full virial tensor, including the reciprocal part, are computed.
• Both methods run on CPU (OpenMP) and GPU (CUDA or HIP). PPPM FFTs use cuFFT on NVIDIA GPUs, hipFFT on AMD GPUs and the bundled pocketfft library on CPU. In a HIP build, hipFFT is found in the ROCm installation (ROCM_INSTALL_ROOT, default /opt/rocm); if it is missing, the configure step says so and PPPM runs on CPU.
• When the cell changes (NPT, deformation), the reciprocal space data is updated automatically at the next call of the initialization operator. Place that operator in compute_force when the cell changes during the run.

All accuracies are relative: accuracy_relative is the target RMS force error divided by the force between two unit charges 1 Å apart (14.399645 eV/Å). The Coulomb constant is \(1/(4\pi\varepsilon_0) = 14.399645\) eV·Å/e².

Ewald Summation¶

coulombic_ewald_init builds the list of k vectors. coulombic_ewald_long_range computes the structure factors \(S(\mathbf{k})\) and the reciprocal space forces, per-atom energies and virial. Only half of k space is stored: \(\mathbf{k}\) and \(-\mathbf{k}\) give the same contribution. For a triclinic cell with matrix \(H\) (columns = cell vectors), the k vectors are \(\mathbf{k} = 2\pi H^{-T}\mathbf{n}\).

coulombic_ewald_init:
  accuracy_relative: 1.0e-5
  g_ewald: 0.0            # automatic
  radius: 10.0 ang
  kmax: 0                 # automatic

init_parameters:
  - coulombic_ewald_init  # sets rcut / rcut_max before the system exists

compute_force:
  - coulombic_ewald_short_range
  - coulombic_ewald_long_range

setup_system:
  # ... domain, particles ...
  - copy_charge_species_to_particle
  - sum_charges
  - coulombic_ewald_init

The chosen parameters (g_ewald, kmax per direction, number of k vectors in half k space) are printed at initialization in an Ewald configuration block. For the 12000-atom UO2 validation system with \(r_c\) = 10 Å and accuracy_relative \(10^{-5}\), the automatic choice is g_ewald = 0.32644139 with 14335 half-space k vectors.

Particle-Particle Particle-Mesh¶

coulombic_pppm spreads the charges onto a regular mesh with a B-spline of order order, then solves Poisson's equation by FFT with the optimal influence function (Hockney–Eastwood). The field is interpolated back onto the particles. coulombic_pppm_init chooses g_ewald and the mesh size from accuracy_relative (mesh sizes with factors 2, 3 and 5 only). coulombic_pppm_init also fills the parameters used by coulombic_ewald_short_range, so the real space part is shared with the Ewald method.

coulombic_pppm_init:
  accuracy_relative: 1.0e-5
  g_ewald: 0.0            # automatic
  radius: 10.0 ang
  order: 5
  mesh: [ 0, 0, 0 ]       # automatic

init_parameters:
  - coulombic_pppm_init

compute_force:
  - coulombic_ewald_short_range
  - coulombic_pppm

setup_system:
  # ... domain, particles ...
  - copy_charge_species_to_particle
  - sum_charges
  - coulombic_pppm_init

The chosen parameters (g_ewald, mesh, order, differentiation, decomposition, slab factor, estimated accuracy) are printed at initialization in a PPPM configuration block. For the 12000-atom UO2 validation system, the automatic choice is g_ewald = 0.33415064 with a 72³ mesh (72×75×75 for the triclinic cell).

Differentiation: ik and ad¶

• ik (default) computes the field \(\mathbf{E}(\mathbf{k}) = -i\mathbf{k}\,\phi(\mathbf{k})\) in k space. It works on orthogonal and triclinic cells.
• ad (analytic differentiation) differentiates the assignment function in real space. It needs fewer inverse FFTs and a larger mesh for the same accuracy, and it works on orthogonal cells only. The self-force correction (sf_coeff) is applied.

Slab correction¶

For systems that are periodic in x and y only (surfaces, interfaces), slab / slab_auto add the EW3DC correction of Yeh and Berkowitz:

• The z boundary of the domain must be non-periodic, and x and y must be periodic.
• The mesh is built on a cell extended by the factor slab in z (typically 3.0). The dipole correction removes the interaction with the periodic images in z.
• With slab_auto: true, the extension factor is computed from accuracy_relative and g_ewald.
• Triclinic cells are supported with an xy tilt only: the third cell vector must be along z (\(H_{13}=H_{23}=H_{31}=H_{32}=0\)).
• The slab correction adds forces and energies, but no virial contribution.

coulombic_pppm_init:
  accuracy_relative: 1.0e-5
  radius: 10.0 ang
  slab: 3.0               # or slab_auto: true

setup_system:
  - domain:
      periodic: [ true, true, false ]
      # ...

Mesh decomposition across MPI ranks¶

Both decompositions give the same results. The cost of replicated grows with the number of ranks, since every rank runs the full FFTs and the whole mesh is summed over all ranks. Time per step of coulombic_pppm on 12000 atoms, 72³ mesh, 1 thread per rank, energy and virial computed every step:

On 2 ranks the two are close; from 4 ranks on, distributed is faster. On 1 rank the distributed layout only adds data transposes (30.0 vs 20.9 ms above), so the default distributed runs the replicated path on a single rank. The PPPM configuration block printed at initialization shows the decomposition actually used.

Real space part¶

coulombic_ewald_short_range computes the \(\operatorname{erfc}(g r)/r\) pair term. It is used by both Ewald and PPPM. It reads the parameters set by coulombic_ewald_init or coulombic_pppm_init, so it has no physics parameter of its own.

With use_symmetry: true alone, the pair loop also visits ghost cells. With ghost_fold_back: true, it visits owned cells only. On 12000 atoms, one thread, \(r_c\) = 10 Å, this brings the real space time from 140 to 89 ms/step. The ghost forces, energies and virials must then be zeroed before the force computation and folded back after it:

coulombic_ewald_short_range:
  use_symmetry: true
  ghost_fold_back: true

compute_force_prolog:
  - zero_force_energy: { ghost: true }

compute_force_epilog:
  - update_virial_force_energy_from_ghost
  - force_to_accel

# optional with use_symmetry: half neighbor lists, halves the neighbor list memory
chunk_neighbors:
  config:
    half_symmetric: true

Setting ghost_fold_back: true without use_symmetry: true is a fatal error.

Tests and timings¶

The reference test system is coulomb-only UO2 with 12000 atoms (disturbed fluorite, \(r_c\) = 10 Å, accuracy \(10^{-5}\)), orthogonal and triclinic (tilts xy=5, xz=3, yz=4 Å). All Ewald and PPPM options were checked on it (Ewald automatic and user kmax; PPPM ik, ad, slab, slab_auto, triclinic slab; per-atom and species charges, use_symmetry, ghost_fold_back), on 1 rank, 2 MPI ranks × 2 OpenMP threads, GPU with 1 and 2 ranks, and with the distributed PPPM mesh on up to 8 MPI ranks. The distributed mesh gives identical results on 3, 4 and 8 MPI ranks.

The ctest regression cases in data/regression_new/potentials/coulombic/ (ewald_*, pppm*, *_fold) run the same operators on a smaller generated system (768 atoms).

Timings¶

12000 atoms UO2, energy and virial computed every step, time per step of the reciprocal space operator. CPU is 1 thread.