chemx-props 0.3.0

//! Electric dipole moment from the one-particle density.

use chemx_core::Molecule;
use chemx_integrals::IntegralProvider;

/// Electric dipole moment in atomic units (electron·bohr).
///
/// μ_k = Σ_A Z_A (R_{A,k} − O_k)  (nuclear)
///      − Σ_μν D_μν ⟨μ|(r−O)_k|ν⟩  (electronic, charge = −1)
///
/// `density` is the **total** AO density (Dα + Dβ), row-major n×n.
/// `origin` is the gauge origin; use the center of mass for translational
/// invariance. The return value is translation-invariant for neutral molecules
/// regardless of origin.
///
/// Multiply each component by [`chemx_core::units::AU_DIPOLE_TO_DEBYE`] to
/// convert to Debye.
pub fn dipole_moment<P: IntegralProvider>(
    provider: &P,
    molecule: &Molecule,
    density: &[f64],
    origin: [f64; 3],
) -> [f64; 3] {
    let n = provider.n_basis();
    let [dx, dy, dz] = provider.dipole_integrals(origin);

    let mut mu = [0.0f64; 3];

    // Electronic contribution: −Tr[D · r_k]
    for ao_mu in 0..n {
        for ao_nu in 0..n {
            let d = density[ao_mu * n + ao_nu];
            mu[0] -= d * dx[ao_mu * n + ao_nu];
            mu[1] -= d * dy[ao_mu * n + ao_nu];
            mu[2] -= d * dz[ao_mu * n + ao_nu];
        }
    }

    // Nuclear contribution: +Σ_A Z_A (R_A − O)
    for atom in &molecule.atoms {
        let z = atom.element.z() as f64;
        let r = atom.position;
        mu[0] += z * (r[0] - origin[0]);
        mu[1] += z * (r[1] - origin[1]);
        mu[2] += z * (r[2] - origin[2]);
    }

    mu
}

/// Center of mass of the molecule in bohr (for use as the dipole origin).
pub fn center_of_mass(molecule: &Molecule) -> [f64; 3] {
    let mut com = [0.0f64; 3];
    let mut total_mass = 0.0f64;
    for atom in &molecule.atoms {
        let m = atom.element.mass();
        total_mass += m;
        for (k, c) in com.iter_mut().enumerate() {
            *c += m * atom.position[k];
        }
    }
    for c in &mut com {
        *c /= total_mass;
    }
    com
}