Features

Coulomb is a library for working with electrolyte solutions and calculating electrostatic interactions in and between molecules and particles. The main purpose is to offer support for molecular simulation software.

• Temperature dependent dielectric permittivity models for common solvents.
• Handling of ionic strength and the Debye screening length.
• Automatic stoichiometry deduction for arbitrary salts.
• Extensive library of truncated electrostatic interaction schemes such as Wolf, Reaction field, Real-space Ewald, generalized through a short-range function trait.
• Reciprocal-space Ewald summation for Coulomb and Yukawa potentials (PBC and IPBC), with automatic parameter optimization, mixed f32/f64 precision, and SIMD acceleration on x86_64.
• Ewald summation with and without implicit salt.
• Multipole expansion for energies, forces, fields between ions, dipoles, and quadrupoles.
• Extensively unit tested and may serve as reference for other implementations or approximations.
• Partial support for static unit of measure analysis via uom. To enable, use the uom feature flag.
• Vector types use mint for interoperability with glam, cgmath, nalgebra, and other math libraries.

This is largely a Rust rewrite and extension of the CoulombGalore C++ library.

Examples

Dielectric Media and Electrolytes

Simple polynomial models are provided to obtain the relative permittivity or a Medium as a function of temperature. For working with the ionic strength, Salt of arbitrary valency can be given and the required stoichiometry is automatically worked out.

use coulomb::{Medium, Salt};
let molarity = 0.1;
let medium = Medium::salt_water(298.15, Salt::CalciumChloride, molarity);
assert_eq!(medium.permittivity()?, 78.35565171480539);
assert_eq!(medium.ionic_strength()?, 0.3);             // mol/l
assert_eq!(medium.debye_length()?, 5.548902662386284); // angstrom

Multipolar Interactions

All pairwise schemes support calculation of potential, energy, field, force from or between multipolar particles, up to second order (ion-ion, ion-dipole, dipole-dipole; ion-quadrupole). Most scheme can be evaluated with or without a Debye-Hückel screening length.

Vector arguments accept any type implementing Into<mint::Vector3<f64>>, including simple arrays and vectors from nalgebra, glam, cgmath, and other libraries with mint support.

use coulomb::pairwise::*;
let scheme = Plain::default();                     // Vanilla Coulomb scheme, 𝒮(𝑞)=1
let z = 1.0;                                       // point charge, 𝑧 
let r = [3.0, 5.0, 0.0];                           // distance vector, 𝒓
let r_len = (r[0]*r[0] + r[1]*r[1] + r[2]*r[2]).sqrt();
let ion_pot = scheme.ion_potential(z, r_len);      // potential |𝒓| away from charge 
assert_eq!(ion_pot, z / r_len);

let mu = [0.2, 3.0, -1.0];                         // point dipole, 𝝁
let dipole_pot = scheme.dipole_potential(mu, r);   // potential 𝒓 away from dipole
let energy = scheme.ion_dipole_energy(z, mu, r);   // interaction energy assuming 𝒓 = 𝒓(𝜇) - 𝒓(𝑧)

The image below is generated by examples/potential.rs and shows the calculated electric potential on a plane containing a monopole and a dipole.

Reciprocal-Space Ewald Summation

The reciprocal module provides k-space Ewald summation for Coulomb (κ=0) and Yukawa (κ>0) potentials with automatic parameter selection, O(N_k) Monte Carlo trial moves, and mixed-precision SIMD acceleration.

use coulomb::reciprocal::*;

// Automatic parameter selection for cutoff 5 Å and target accuracy 1e-5
let mut ewald = EwaldReciprocal::new([10.0; 3], 5.0, 1e-5, None);

let (x, y, z) = (vec![-0.5, 0.5], vec![0.0, 0.0], vec![0.0, 0.0]);
let charges = vec![1.0, -1.0];

// Build structure factors from all particle positions
ewald.update_all(&x, &y, &z, &charges);

// Reciprocal-space energy + self-energy correction
let energy = ewald.energy() + ewald.self_energy(&charges);

// MC trial move: compute ΔE without mutating state
let de = ewald.energy_change(1.0, [-0.5, 0.0, 0.0], [-0.3, 0.0, 0.0]);

// Forces for MD
let (mut fx, mut fy, mut fz) = (vec![0.0; 2], vec![0.0; 2], vec![0.0; 2]);
ewald.add_forces(&x, &y, &z, &charges, &mut fx, &mut fy, &mut fz);

Both f32 and f64 input slices are accepted (generic over Into<f64> + Copy). On x86_64 with the simd feature, a polynomial sin/cos with SIMD vectorization accelerates the inner loops.

Unit analysis

Experimental support for static unit analysis can be activated with the uom feature.

use coulomb::{units::*, pairwise::{Plain, MultipoleEnergySI}};
let z1 = ElectricCharge::new::<elementary_charge>(1.0);
let z2 = ElectricCharge::new::<elementary_charge>(2.0);
let r = Length::new::<nanometer>(2.3);
let energy = Plain::without_cutoff().ion_ion_energy(z1, z2, r);
assert_eq!(energy.get::<kilojoule_per_mole>(), 362.4403242896922);