coulomb 0.5.0

// Copyright 2023 Björn Stenqvist and Mikael Lund
//
// Converted to Rust with modification from the C++ library "CoulombGalore":
// https://zenodo.org/doi/10.5281/zenodo.3522058
//
// Licensed under the Apache license, version 2.0 (the "license");
// you may not use this file except in compliance with the license.
// You may obtain a copy of the license at
//
//     http://www.apache.org/licenses/license-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the license is distributed on an "as is" basis,
// without warranties or conditions of any kind, either express or implied.
// See the license for the specific language governing permissions and
// limitations under the license.

use crate::debye_length::DebyeLength;
use crate::pairwise::{SelfEnergyPrefactors, ShortRangeFunction};
#[cfg(feature = "serde")]
use serde::{Deserialize, Deserializer, Serialize};

/// Fixed capacity of the compile-time coefficient table ([`Poisson::COEFFS`]).
///
/// The table is sized to a constant rather than to `C` because stable Rust
/// cannot use a const-generic expression as an array length: `[_; C as usize]`
/// requires the unstable `generic_const_exprs`. Every published scheme has
/// `C <= 4`, so 8 is ample headroom; a larger `C` trips the assert in `COEFFS`.
const MAX_C: usize = 8;

/// Binomial coefficient `n` over `k`, usable in `const` context.
///
/// Returns 0 when `k < 0`, `n < 0`, or `k > n` — matching the only cases that
/// arise when building the Poisson coefficient table. (`num_integer::binomial`
/// is not `const`, hence this small local version.)
const fn binomial(mut n: i64, k: i64) -> i64 {
    if k < 0 || n < 0 {
        return 0;
    }
    let mut numerator: i64 = 1;
    let mut denominator: i64 = 1;
    let mut i: i64 = 0;
    while i < k {
        numerator *= n;
        denominator *= i + 1;
        n -= 1;
        i += 1;
    }
    numerator / denominator
}

/// # Scheme for the Poisson short-range function
///
/// This is a general scheme for the short-ranged part of the electrostatic interaction
/// which can be used to arbitrarily cancel derivatives at the origin and at the cut-off.
/// From the abstract of <https://doi.org/c5fr>:
///
/// _"Electrostatic pair-potentials within molecular simulations are often based on empirical data,
/// cancellation of derivatives or moments, or statistical distributions of image-particles.
/// In this work we start with the fundamental Poisson equation and show that no truncated Coulomb
/// pair-potential, unsurprisingly, can solve the Poisson equation. For any such pair-potential
/// the Poisson equation gives two incompatible constraints, yet we find a single unique expression
/// which, pending two physically connected smoothness parameters, can obey either one of these.
/// This expression has a general form which covers several recently published pair-potentials.
/// For sufficiently large degree of smoothness we find that the solution implies a Gaussian
/// distribution of the charge, a feature which is frequently assumed in pair-potential theory.
/// We end up by recommending a single pair-potential based both on theoretical arguments and
/// empirical evaluations of non-thermal lattice- and thermal water-systems.
/// The same derivations have also been made for the screened Poisson equation,
/// i.e. for Yukawa potentials, with a similar solution."_
///
/// The general short-range function is:
/// $$
/// S(q) = (1 - q)^{D + 1} \sum_{c = 0}^{C - 1} \frac{C - c}{C} \binom{D - 1 + c}{c} q^c
/// $$
///
/// where $C$ is the number of cancelled derivatives at origin -2 (starting from the second derivative),
/// and $D$ is the number of cancelled derivatives at the cut-off (starting from the zeroth derivative).
///
/// For infinite Debye-length, $\kappa=0$, the [`Poisson`] scheme captures several
/// other truncation schemes by setting $C$ and $D$ according to this table:
///
/// | Type          | $C$ | $D$ | Reference / Comment
/// |---------------|-----|-----|---------------------
/// | `plain`       | 1   | -1  | Scheme for a vanilla coulomb interaction using the Poisson framework. Same as `Coulomb`.
/// | `wolf`        | 1   | 0   | Scheme for [Undamped Wolf](https://doi.org/10.1063/1.478738)
/// | `fennell`     | 1   | 1   | Scheme for [Levitt/undamped Fennell](https://doi.org/10/fp959p). See also doi:10/bqgmv2.
/// | `kale`        | 1   | 2   | Scheme for [Kale](https://doi.org/10/csh8bg)
/// | `mccann`      | 1   | 3   | Scheme for [McCann](https://doi.org/10.1021/ct300961)
/// | `fukuda`      | 2   | 1   | Scheme for [Undamped Fukuda](https://doi.org/10.1063/1.3582791)
/// | `markland`    | 2   | 2   | Scheme for [Markland](https://doi.org/10.1016/j.cplett.2008.09.019)
/// | `stenqvist`   | 3   | 3   | Scheme for [Stenqvist](https://doi.org/10/c5fr)
/// | `fanourgakis` | 4   | 3   | Scheme for [Fanourgakis](https://doi.org/10.1063/1.3216520),
///
/// Helper struct to store salt screening parameters
#[derive(Debug, Clone, PartialEq)]
#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
struct Screening {
    /// Inverse Debye screening length
    #[cfg_attr(feature = "serde", serde(alias = "κ"))]
    pub kappa: f64,
    /// Reduced kappa = cutoff * kappa
    pub reduced_kappa: f64,
    pub reduced_kappa_squared: f64,
    pub yukawa_denom: f64,
}

impl Screening {
    fn new(debye_length: f64, cutoff: f64) -> Self {
        let reduced_kappa = cutoff / debye_length;
        Screening {
            kappa: 1.0 / debye_length,
            reduced_kappa,
            reduced_kappa_squared: reduced_kappa.powi(2),
            yukawa_denom: 1.0 / (1.0 - (2.0 * reduced_kappa).exp()),
        }
    }
}

#[derive(Debug, Clone, PartialEq)]
#[cfg_attr(feature = "serde", derive(Serialize))]
/// Generalized Poisson scheme for short-range electrostatic interactions.
pub struct Poisson<const C: i32, const D: i32> {
    /// Cutoff radius
    cutoff: f64,
    /// Debye length
    #[cfg_attr(feature = "serde", serde(alias = "debyelength"))]
    debye_length: Option<f64>,
    /// Currently not in use
    #[cfg_attr(feature = "serde", serde(skip))]
    _has_dipolar_selfenergy: bool,
    #[cfg_attr(feature = "serde", serde(skip))]
    screening: Option<Screening>,
}

#[cfg(feature = "serde")]
impl<'de, const C: i32, const D: i32> Deserialize<'de> for Poisson<C, D> {
    fn deserialize<DES>(deserializer: DES) -> Result<Self, DES::Error>
    where
        DES: Deserializer<'de>,
    {
        #[derive(Deserialize)]
        #[serde(deny_unknown_fields)]
        struct PoissonData {
            cutoff: f64,
            debye_length: Option<f64>,
        }

        let PoissonData {
            cutoff,
            debye_length,
        } = PoissonData::deserialize(deserializer)?;
        Ok(Poisson::new(cutoff, debye_length))
    }
}

/// Scheme for a vanilla coulomb interaction using the Poisson framework. Same as `Coulomb`.
pub type _Plain = Poisson<1, -1>;

/// Energy and force shifted Yukawa potential [Levitt/undamped Fennell](https://doi.org/10/fp959p).
///
/// See also doi:10/bqgmv2.
pub type Yukawa = Poisson<1, 1>;

/// Scheme for [Undamped Wolf](https://doi.org/10.1063/1.478738)
pub type UndampedWolf = Poisson<1, 0>;

/// Scheme for [Kale](https://doi.org/10/csh8bg)
pub type Kale = Poisson<1, 2>;

/// Scheme for [McCann](https://doi.org/10.1021/ct300961)
pub type McCann = Poisson<1, 3>;

/// Scheme for [Undamped Fukuda](https://doi.org/10.1063/1.3582791)
pub type UndampedFukuda = Poisson<2, 1>;

/// Scheme for [Markland](https://doi.org/10.1016/j.cplett.2008.09.019)
pub type Markland = Poisson<2, 2>;

/// Scheme for [Stenqvist](https://doi.org/10/c5fr)
pub type Stenqvist = Poisson<3, 3>;

/// Scheme for [Fanourgakis](https://doi.org/10.1063/1.3216520)
pub type Fanourgakis = Poisson<4, 3>;

impl<const C: i32, const D: i32> Poisson<C, D> {
    /// Create a new Poisson scheme with a given cutoff and optional Debye length.
    pub fn new(cutoff: f64, debye_length: Option<f64>) -> Self {
        if C < 1 {
            panic!("`C` must be larger than zero");
        }
        if D < -1 && D != -C {
            panic!("If `D` is less than negative one, then it has to equal negative `C`");
        }
        if D == 0 && C != 1 {
            panic!("If `D` is zero, then `C` has to equal one ");
        }

        let _has_dipolar_selfenergy = C >= 2;

        let screening = debye_length.map(|d| Screening::new(d, cutoff));

        Self {
            cutoff,
            debye_length,
            _has_dipolar_selfenergy,
            screening,
        }
    }

    /// Integer numerators `bᶜ = (D−1+c choose c)·(C−c)` of the short-range
    /// polynomial, evaluated at compile time and shared by every instance of a
    /// given `(C, D)`. The `÷C` is deferred to evaluation time, keeping the
    /// table pure-integer (and thus `const`-evaluable on stable Rust).
    const COEFFS: [i64; MAX_C] = {
        assert!(C >= 1, "Poisson: C must be >= 1");
        assert!(
            C as usize <= MAX_C,
            "Poisson: C exceeds the coefficient-table cap MAX_C; increase MAX_C"
        );
        // Bound `D` so the falling-factorial numerators in `binomial` cannot
        // overflow `i64` during const-eval (which would otherwise surface as a
        // cryptic E0080). 12 sits far above every published scheme (D <= 3) and
        // is safe for all C <= MAX_C; raising it requires wider intermediates.
        assert!(
            D <= 12,
            "Poisson: D exceeds the i64 coefficient-table limit (12)"
        );
        let mut coeffs = [0i64; MAX_C];
        let mut c = 0i32;
        while c < C {
            coeffs[c as usize] = binomial((D - 1 + c) as i64, c as i64) * (C - c) as i64;
            c += 1;
        }
        coeffs
    };

    /// `(C+D choose C)·D`, evaluated at compile time (used by the 2nd and 3rd
    /// derivatives). When `D == -C` this is zero and the derivatives short-circuit
    /// before it is read.
    const BINOM_CDC: i64 = binomial(C as i64 + D as i64, C as i64) * D as i64;

    /// Short-range polynomial `S̃(qp) = Σ_{c=0}^{C-1} (bᶜ/C)·qpᶜ` and its
    /// derivative `dS̃/dqp`, via a single (simultaneous) Horner pass over the
    /// compile-time [`Self::COEFFS`].
    #[inline]
    fn poly_and_deriv(qp: f64) -> (f64, f64) {
        let divc = C as f64;
        let mut poly = Self::COEFFS[(C - 1) as usize] as f64 / divc;
        let mut deriv = 0.0;
        for c in (0..C - 1).rev() {
            deriv = deriv * qp + poly;
            poly = poly * qp + Self::COEFFS[c as usize] as f64 / divc;
        }
        (poly, deriv)
    }

    /// First derivative `dS/dqp` of the full short-range function `S = (1-qp)^{D+1}·S̃`.
    #[inline]
    fn dsdqp(qp: f64) -> f64 {
        let (poly, dpoly) = Self::poly_and_deriv(qp);
        -f64::from(D + 1) * (1.0 - qp).powi(D) * poly + (1.0 - qp).powi(D + 1) * dpoly
    }

    /// Second derivative `d²S/dqp²`, which collapses to a single monomial.
    #[inline]
    fn d2sdqp2(qp: f64) -> f64 {
        Self::BINOM_CDC as f64 * (1.0 - qp).powi(D - 1) * qp.powi(C - 1)
    }
}

impl<const C: i32, const D: i32> crate::Cutoff for Poisson<C, D> {
    fn cutoff(&self) -> f64 {
        self.cutoff
    }
}

impl<const C: i32, const D: i32> crate::DebyeLength for Poisson<C, D> {
    #[inline]
    fn kappa(&self) -> Option<f64> {
        self.screening.as_ref().map(|s| s.kappa)
    }
    fn set_debye_length(&mut self, debye_length: Option<f64>) -> crate::Result<()> {
        self.screening = debye_length.map(|d| Screening::new(d, self.cutoff));
        Ok(())
    }
}

impl<const C: i32, const D: i32> ShortRangeFunction for Poisson<C, D> {
    fn url() -> &'static str {
        match (C, D) {
            (1, -1) => "https://doi.org/msxd",             // plain
            (1, 0) => "https://doi.org/10.1063/1.478738",  // wolf
            (1, 1) => "https://doi.org/10/fp959p",         // fennell
            (1, 2) => "https://doi.org/10/csh8bg",         // kale
            (1, 3) => "https://doi.org/10.1021/ct300961",  // mccann
            (2, 1) => "https://doi.org/10.1063/1.3582791", // fukuda
            (2, 2) => "https://doi.org/dbpbts",            // markland
            (3, 3) => "https://doi.org/10/c5fr",           // stenqvist
            (4, 3) => "https://doi.org/10.1063/1.3216520", // fanourgakis
            _ => "https://doi.org/c5fr",                   // generic poisson
        }
    }

    fn short_range_f0(&self, q: f64) -> f64 {
        // checks on C and D are evaluated at compile time
        if D == -C {
            return 1.0;
        }
        let qp: f64 = self.screening.as_ref().map_or(q, |s| {
            (1.0 - (2.0 * s.reduced_kappa * q).exp()) * s.yukawa_denom
        });

        if D == 0 && C == 1 {
            return 1.0 - qp;
        }

        let (poly, _) = Self::poly_and_deriv(qp);
        (1.0 - qp).powi(D + 1) * poly
    }
    fn short_range_f1(&self, q: f64) -> f64 {
        if D == -C {
            return 0.0;
        }
        if D == 0 && C == 1 {
            return 0.0;
        }
        let (qp, dqpdq) = if let Some(s) = &self.screening {
            let exp2kq = (2.0 * s.reduced_kappa * q).exp();
            let qp = (1.0 - exp2kq) * s.yukawa_denom;
            let dqpdq = -2.0 * s.reduced_kappa * exp2kq * s.yukawa_denom;
            (qp, dqpdq)
        } else {
            (q, 1.0)
        };
        Self::dsdqp(qp) * dqpdq
    }

    fn short_range_f2(&self, q: f64) -> f64 {
        if D == -C {
            return 0.0;
        }
        if D == 0 && C == 1 {
            return 0.0;
        }

        let (qp, dqpdq, d2qpdq2, dsdqp) = if let Some(s) = &self.screening {
            let exp2kq = (2.0 * s.reduced_kappa * q).exp();
            let qp = (1.0 - exp2kq) * s.yukawa_denom;
            let dqpdq = -2.0 * s.reduced_kappa * exp2kq * s.yukawa_denom;
            let d2qpdq2 = -4.0 * s.reduced_kappa_squared * exp2kq * s.yukawa_denom;
            (qp, dqpdq, d2qpdq2, Self::dsdqp(qp))
        } else {
            (q, 1.0, 0.0, 0.0)
        };
        Self::d2sdqp2(qp) * dqpdq * dqpdq + dsdqp * d2qpdq2
    }

    fn short_range_f3(&self, q: f64) -> f64 {
        if D == -C {
            return 0.0;
        }
        if D == 0 && C == 1 {
            return 0.0;
        }

        let binom_cdc = Self::BINOM_CDC as f64;
        let (qp, dqpdq, d2qpdq2, d3qpdq3, d2sdqp2, dsdqp) = if let Some(s) = &self.screening {
            let exp2kq = (2.0 * s.reduced_kappa * q).exp();
            let qp = (1.0 - exp2kq) * s.yukawa_denom;
            let dqpdq = -2.0 * s.reduced_kappa * exp2kq * s.yukawa_denom;
            let d2qpdq2 = -4.0 * s.reduced_kappa_squared * exp2kq * s.yukawa_denom;
            let d3qpdq3 =
                -8.0 * s.reduced_kappa_squared * s.reduced_kappa * exp2kq * s.yukawa_denom;
            (qp, dqpdq, d2qpdq2, d3qpdq3, Self::d2sdqp2(qp), Self::dsdqp(qp))
        } else {
            (q, 1.0, 0.0, 0.0, 0.0, 0.0)
        };
        let d3sdqp3 = binom_cdc
            * (1.0 - qp).powi(D - 2)
            * qp.powi(C - 2)
            * ((2.0 - C as f64 - D as f64) * qp + C as f64 - 1.0);
        d3sdqp3 * dqpdq * dqpdq * dqpdq + 3.0 * d2sdqp2 * dqpdq * d2qpdq2 + dsdqp * d3qpdq3
    }

    fn self_energy_prefactors(&self) -> SelfEnergyPrefactors {
        let mut c1: f64 = -0.5 * (C + D) as f64 / C as f64;
        if let Some(s) = &self.screening {
            c1 = c1 * -2.0 * s.reduced_kappa * s.yukawa_denom;
        }
        SelfEnergyPrefactors {
            monopole: Some(c1),
            dipole: None,
        }
    }
}

impl<const C: i32, const D: i32> core::fmt::Display for Poisson<C, D> {
    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
        write!(
            f,
            "Poisson: 𝐶 = {}, 𝐷 = {}, 𝑟✂ = {:.1} Å",
            C, D, self.cutoff
        )?;
        if let Some(debye_length) = self.kappa().map(f64::recip) {
            write!(f, ", λᴰ = {:.1} Å", debye_length)?;
        }
        write!(f, " <{}>", Self::url())?;
        Ok(())
    }
}

#[test]
fn test_poisson() {
    use super::test_utils::{assert_vec3_eq, assert_vec_x_equals_norm, assert_vec_zero};
    use crate::{
        pairwise::{MultipoleEnergy, MultipoleField, MultipoleForce, MultipolePotential},
        NalgebraMatrix3 as Matrix3, NalgebraVector3 as Vector3,
    };
    use approx::assert_relative_eq;
    let cutoff = 29.0;
    let pot = Stenqvist::new(cutoff, None);
    let eps = 1e-9; // Set epsilon for approximate equality

    // Test Stenqvist short-range function
    assert_relative_eq!(pot.short_range_f0(0.5), 0.15625, epsilon = eps);
    assert_relative_eq!(pot.short_range_f1(0.5), -1.0, epsilon = eps);
    assert_relative_eq!(pot.short_range_f2(0.5), 3.75, epsilon = eps);
    assert_relative_eq!(pot.short_range_f3(0.5), 0.0, epsilon = eps);
    assert_relative_eq!(pot.short_range_f3(0.6), -5.76, epsilon = eps);
    assert_relative_eq!(pot.short_range_f0(1.0), 0.0, epsilon = eps);
    assert_relative_eq!(pot.short_range_f1(1.0), 0.0, epsilon = eps);
    assert_relative_eq!(pot.short_range_f2(1.0), 0.0, epsilon = eps);
    assert_relative_eq!(pot.short_range_f3(1.0), 0.0, epsilon = eps);
    assert_relative_eq!(pot.short_range_f0(0.0), 1.0, epsilon = eps);
    assert_relative_eq!(pot.short_range_f1(0.0), -2.0, epsilon = eps);
    assert_relative_eq!(pot.short_range_f2(0.0), 0.0, epsilon = eps);
    assert_relative_eq!(pot.short_range_f3(0.0), 0.0, epsilon = eps);

    let pot = Stenqvist::new(cutoff, Some(23.0));
    approx::assert_relative_eq!(
        pot.self_energy(&[2.0], &[0.0]),
        -0.03037721287,
        epsilon = eps
    );

    assert_eq!(
        pot.to_string(),
        "Poisson: 𝐶 = 3, 𝐷 = 3, 𝑟✂ = 29.0 Å, λᴰ = 23.0 Å <https://doi.org/10/c5fr>"
    );

    // Test Fanougarkis short-range function
    let pot = Fanourgakis::new(cutoff, None);
    assert_relative_eq!(pot.short_range_f0(0.5), 0.19921875, epsilon = eps);
    assert_relative_eq!(pot.short_range_f1(0.5), -1.1484375, epsilon = eps);
    assert_relative_eq!(pot.short_range_f2(0.5), 3.28125, epsilon = eps);
    assert_relative_eq!(pot.short_range_f3(0.5), 6.5625, epsilon = eps);

    assert_eq!(
        pot.to_string(),
        "Poisson: 𝐶 = 4, 𝐷 = 3, 𝑟✂ = 29.0 Å <https://doi.org/10.1063/1.3216520>"
    );

    let pot = Poisson::<4, 3>::new(cutoff, None);
    let z1 = 2.0;
    let z2 = 3.0;
    let r = Vector3::new(23.0, 0.0, 0.0); // distance vector
    let rq = Vector3::new(
        5.75 * 6.0_f64.sqrt(),
        5.75 * 2.0_f64.sqrt(),
        11.5 * 2.0_f64.sqrt(),
    );
    let rh = r.normalize();
    let mu1 = Vector3::new(19.0, 7.0, 11.0);
    let mu2 = Vector3::new(13.0, 17.0, 5.0);
    //      [3 7 8]
    //  Q = [5 9 6]
    //      [2 1 4]
    let quad1 = Matrix3::from_row_slice(&[3.0, 7.0, 8.0, 5.0, 9.0, 6.0, 2.0, 1.0, 4.0]);
    let quad2 = Matrix3::zeros();
    assert_relative_eq!(quad1.trace(), 16.0, epsilon = eps);

    // Test potentials
    assert_relative_eq!(pot.ion_potential(z1, cutoff), 0.0, epsilon = eps);
    assert_relative_eq!(
        pot.ion_potential(z1, r.norm()),
        0.0009430652121,
        epsilon = eps
    );
    assert_relative_eq!(
        pot.dipole_potential(mu1, rh.scale(cutoff)),
        0.0,
        epsilon = eps
    );
    assert_relative_eq!(pot.dipole_potential(mu1, r), 0.005750206554, epsilon = eps);
    assert_relative_eq!(
        pot.quadrupole_potential(quad1, rq),
        0.000899228165,
        epsilon = eps
    );
    assert_relative_eq!(
        pot.quadrupole_potential(quad1, rq.scale(cutoff / 23.0)),
        0.0,
        epsilon = eps
    );

    // Test fields with no salt
    assert_vec_zero!(pot.ion_field(z1, rh.scale(cutoff)), eps);
    assert_vec_x_equals_norm!(pot.ion_field(z1, r), 0.0006052849004, eps);
    assert_vec_zero!(pot.dipole_field(mu1, rh.scale(cutoff)), eps);

    let ion_field_scalar = pot.ion_field_scalar(z1, r.norm());
    let e_ion: Vector3 = pot.ion_field(z1, r).into();
    assert_relative_eq!(ion_field_scalar, e_ion.norm(), epsilon = eps);

    assert_vec3_eq!(
        pot.dipole_field(mu1, r),
        [0.002702513754, -0.00009210857180, -0.0001447420414],
        eps
    );
    assert_vec3_eq!(
        pot.quadrupole_field(quad1, r),
        [0.00001919309993, -0.00004053806958, -0.00003378172465],
        eps
    );

    // Test energies
    assert_relative_eq!(pot.ion_ion_energy(z1, z2, cutoff), 0.0, epsilon = eps);
    assert_relative_eq!(
        pot.ion_ion_energy(z1, z2, r.norm()),
        0.002829195636,
        epsilon = eps
    );
    assert_relative_eq!(
        pot.ion_dipole_energy(z1, mu2, rh.scale(cutoff)),
        0.0,
        epsilon = eps
    );
    assert_relative_eq!(
        pot.ion_dipole_energy(z1, mu2, r),
        -0.007868703705,
        epsilon = eps
    );
    assert_relative_eq!(
        pot.ion_dipole_energy(z2, mu1, -r),
        0.01725061966,
        epsilon = eps
    );
    assert_relative_eq!(
        pot.dipole_dipole_energy(mu1, mu2, rh.scale(cutoff)),
        0.0,
        epsilon = eps
    );
    assert_relative_eq!(
        pot.dipole_dipole_energy(mu1, mu2, r),
        -0.03284312288,
        epsilon = eps
    );
    assert_relative_eq!(
        pot.ion_quadrupole_energy(z2, quad1, rq.scale(cutoff / 23.0 * 29.0)),
        0.0,
        epsilon = eps
    );
    assert_relative_eq!(
        pot.ion_quadrupole_energy(z2, quad1, rq),
        0.002697684495,
        epsilon = eps
    );
    assert_relative_eq!(
        pot.ion_quadrupole_energy(z1, quad2, -rq.scale(cutoff / 23.0 * 29.0)),
        0.0,
        epsilon = eps
    );
    assert_relative_eq!(
        pot.ion_quadrupole_energy(z1, quad2, -rq),
        0.0,
        epsilon = eps
    );

    // Test forces
    assert_vec_zero!(
        pot.ion_ion_force(z1, z2, Vector3::new(cutoff, 0.0, 0.0)),
        eps
    );
    assert_vec_x_equals_norm!(pot.ion_ion_force(z1, z2, r), 0.001815854701, eps);
    assert_vec_zero!(pot.ion_dipole_force(z2, mu1, rh.scale(cutoff)), eps);
    assert_vec3_eq!(
        pot.ion_dipole_force(z2, mu1, r),
        [0.008107541263, -0.0002763257154, -0.0004342261242],
        eps
    );
    assert_vec3_eq!(
        pot.ion_dipole_force(z1, mu2, -r),
        [0.003698176716, -0.0004473844916, -0.0001315836740],
        eps
    );
    assert_vec_zero!(pot.dipole_dipole_force(mu1, mu2, rh.scale(cutoff)), eps);
    assert_vec3_eq!(
        pot.dipole_dipole_force(mu1, mu2, r),
        [0.009216400961, -0.002797126801, -0.001608010094],
        eps
    );

    // Test with screening
    let debye_length = 23.0;
    let pot = Poisson::<3, 3>::new(cutoff, Some(debye_length));

    // Test short-ranged function with screening
    assert_relative_eq!(pot.short_range_f0(0.5), 0.5673222086324718, epsilon = eps);
    assert_relative_eq!(pot.short_range_f1(0.5), -1.4373727619264975, epsilon = eps);
    assert_relative_eq!(pot.short_range_f2(0.5), -2.552012309527445, epsilon = eps);
    assert_relative_eq!(pot.short_range_f3(0.5), 4.384434366606605, epsilon = eps);

    // Test potentials with screening
    assert_relative_eq!(pot.ion_potential(z1, cutoff), 0.0, epsilon = eps);
    assert_relative_eq!(
        pot.ion_potential(z1, r.norm()),
        0.003344219306,
        epsilon = eps
    );
    assert_relative_eq!(
        pot.dipole_potential(mu1, rh.scale(cutoff)),
        0.0,
        epsilon = eps
    );
    assert_relative_eq!(pot.dipole_potential(mu1, r), 0.01614089171, epsilon = eps);
    assert_relative_eq!(
        pot.quadrupole_potential(quad1, rq),
        0.0016294707475,
        epsilon = eps
    );
    assert_relative_eq!(
        pot.quadrupole_potential(quad1, rq.scale(cutoff / 23.0 * 29.0)),
        0.0,
        epsilon = eps
    );

    // Test fields with screening
    assert_vec_zero!(pot.ion_field(z1, r.scale(cutoff)), eps);
    assert_vec_x_equals_norm!(pot.ion_field(z1, r), 0.001699041230, eps);
    assert_vec_zero!(pot.dipole_field(mu1, r.scale(cutoff)), eps);

    let ion_field_scalar = pot.ion_field_scalar(z1, r.norm());
    let ion_field: Vector3 = pot.ion_field(z1, r).into();
    assert_relative_eq!(ion_field_scalar, ion_field.norm(), epsilon = eps);

    assert_vec3_eq!(
        pot.dipole_field(mu1, r),
        [0.004956265485, -0.0002585497523, -0.0004062924688],
        eps
    );
    assert_vec3_eq!(
        pot.quadrupole_field(quad1, r),
        [-0.00005233355205, -0.00007768480608, -0.00006473733856],
        eps
    );

    // Test energies with screening
    assert_relative_eq!(pot.ion_ion_energy(z1, z2, cutoff), 0.0, epsilon = eps);
    assert_relative_eq!(
        pot.ion_ion_energy(z1, z2, r.norm()),
        0.01003265793,
        epsilon = eps
    );
    assert_relative_eq!(
        pot.ion_dipole_energy(z1, mu2, r.scale(cutoff)),
        0.0,
        epsilon = eps
    );
    assert_relative_eq!(
        pot.ion_dipole_energy(z1, mu2, r),
        -0.02208753604,
        epsilon = eps
    );
    assert_relative_eq!(
        pot.ion_dipole_energy(z2, mu1, -r),
        0.04842267505,
        epsilon = eps
    );
    assert_relative_eq!(
        pot.dipole_dipole_energy(mu1, mu2, r.scale(cutoff)),
        0.0,
        epsilon = eps
    );
    assert_relative_eq!(
        pot.dipole_dipole_energy(mu1, mu2, r),
        -0.05800464321,
        epsilon = eps
    );
    assert_relative_eq!(
        pot.ion_quadrupole_energy(z2, quad1, rq.scale(cutoff / 23.0 * 29.0)),
        0.0,
        epsilon = eps
    );
    assert_relative_eq!(
        pot.ion_quadrupole_energy(z2, quad1, rq),
        0.004888412229,
        epsilon = eps
    );

    // Test forces with screening
    assert_vec_zero!(pot.ion_ion_force(z1, z2, r.scale(cutoff)), eps);
    assert_vec_x_equals_norm!(pot.ion_ion_force(z1, z2, r), 0.005097123689, eps);
    assert_vec_zero!(pot.ion_dipole_force(z2, mu1, r.scale(cutoff)), eps);
    assert_vec3_eq!(
        pot.ion_dipole_force(z2, mu1, r),
        [0.01486879646, -0.0007756492577, -0.001218877402],
        eps
    );
    assert_vec3_eq!(
        pot.ion_dipole_force(z1, mu2, -r),
        [0.006782258035, -0.001255813082, -0.0003693567885],
        eps
    );
    assert_vec_zero!(pot.dipole_dipole_force(mu1, mu2, r.scale(cutoff)), eps);
    assert_vec3_eq!(
        pot.dipole_dipole_force(mu1, mu2, r),
        [0.002987655323, -0.005360251624, -0.003081497314],
        eps
    );
}