Runs the Smooth Particle Ewald Mesh (SPME) algorithm for n-body simulations with periodic boundary conditions

Original paper describing the SPME method

This library is for Python and Rust. Supports GPU, or thread-pooled CPU.

This has applications primarily in structural biology. For example, molecular dynamics. Compared to other n-body approximations for long-range forces, this has utility when periodic bounday conditions are used. If not using these, for example in cosmology simulations, consider Barnes Hut, or Fast Multipole Methods (FMM) instead.

Uses Rayon to parallelize as thread pools. Support for SIMD (256-bit and 512-bit), is planned. To use on an nVidia GPU, enable the cuda feature in Cargo.toml.

Used by the Daedalus protein viewer and molecular dynamics program, and the Dynamics library.

Uses f32 for Lennard Jones and Coulomb interactions. Energy sums are computed as f64.

Here's an example of use. The Python API is equivalent.

use rayon::prelude::*;
use ewald::{force_coulomb_ewald_real, force_coulomb_ewald_real};

const LONG_RANGE_CUTOFF: f32 = 10.0;

// A bigger α means more damping, and a smaller real-space contribution. (Cheaper real), but larger
// reciprocal load.
const EWALD_ALPHA: f32 = 0.35; // Å^-1. 0.35 is good for cutoff = 10.

impl System {
    // Primary application:
    fn apply_forces(&self) {
        pairs
            .par_iter()
            .map(|(i_0, i_1)| {
                let atom_0 = &self.atoms[i_0];
                let atom_1 = &self.atoms[i_1];
                let diff = atom_1.pos - atom_0.pos;
                let r = diff.magnitude();
                let dir = diff / r;
    
                let mut f = Vec3::zero();
    
                let (f, energy) = force_coulomb_short_range(
                    dir,
                    r,
                    // We include 1/r as it's likely shared between this and Lennard Jones;
                    // improves efficiency.
                    1./r,
                    atom_0.charge,
                    atom_1.charge,
                    // e.g. (8Å, 10Å)
                    LONG_RANGE_CUTOFF,
                    ALPHA,
                );
    
                atom_0.force += f;
                atom_1.force -= f;
            });

        let (recip_forces_per_atom, energy_recip) = self.pme_recip.forces(&atom_posits, &[atom_charges]);
    }

    /// Run this at init, and whenever you update the sim box.
    pub fn regen_pme(&mut self) {
        self.pme_recip = PmeRecip::new((SPME_N, SPME_N, SPME_N), self.cell.extent, EWALD_ALPHA);
    }
}