diff --git a/ffprime/electrostatics/multipole.py b/ffprime/electrostatics/multipole.py
new file mode 100644
index 0000000..eb99102
--- /dev/null
+++ b/ffprime/electrostatics/multipole.py
@@ -0,0 +1,228 @@
+"""
+Electrostatic potential and electric field from atom-centered multipoles.
+
+Implements monopole, dipole, and quadrupole contributions following the
+multipole expansion of the electrostatic potential in atomic units.
+
+References
+----------
+Stone, A.J. "The Theory of Intermolecular Forces", Oxford University Press.
+"""
+
+import numpy as np
+
+
+def monopole_potential(charges, coords, points):
+    """
+    Compute electrostatic potential from atomic monopoles (charges).
+
+    Parameters
+    ----------
+    charges : np.ndarray, shape (N,)
+        Atomic charges in atomic units.
+    coords : np.ndarray, shape (N, 3)
+        Atomic coordinates in atomic units.
+    points : np.ndarray, shape (M, 3)
+        Field points at which to evaluate the potential.
+
+    Returns
+    -------
+    potential : np.ndarray, shape (M,)
+        Electrostatic potential at each field point in atomic units.
+    """
+    charges = np.asarray(charges)
+    coords = np.asarray(coords)
+    points = np.asarray(points)
+    r_vecs = points[:, np.newaxis, :] - coords[np.newaxis, :, :]
+    r = np.linalg.norm(r_vecs, axis=-1)
+    safe_r = np.where(r < 1e-12, np.inf, r)
+    return np.sum(charges[np.newaxis, :] / safe_r, axis=1)
+
+
+def monopole_field(charges, coords, points):
+    """
+    Compute electric field from atomic monopoles (charges).
+
+    Parameters
+    ----------
+    charges : np.ndarray, shape (N,)
+    coords : np.ndarray, shape (N, 3)
+    points : np.ndarray, shape (M, 3)
+
+    Returns
+    -------
+    field : np.ndarray, shape (M, 3)
+    """
+    charges = np.asarray(charges)
+    coords = np.asarray(coords)
+    points = np.asarray(points)
+    r_vecs = points[:, np.newaxis, :] - coords[np.newaxis, :, :]
+    r = np.linalg.norm(r_vecs, axis=-1)
+    safe_r = np.where(r < 1e-12, np.inf, r)
+    field = charges[np.newaxis, :, np.newaxis] * r_vecs / safe_r[:, :, np.newaxis] ** 3
+    return np.sum(field, axis=1)
+
+
+def dipole_potential(dipoles, coords, points):
+    """
+    Compute electrostatic potential from atomic dipoles.
+
+    V(r) = sum_i [ p_i . (r - r_i) ] / |r - r_i|^3
+
+    Parameters
+    ----------
+    dipoles : np.ndarray, shape (N, 3)
+        Atomic dipole moment vectors in atomic units.
+    coords : np.ndarray, shape (N, 3)
+    points : np.ndarray, shape (M, 3)
+
+    Returns
+    -------
+    potential : np.ndarray, shape (M,)
+    """
+    dipoles = np.asarray(dipoles)
+    coords = np.asarray(coords)
+    points = np.asarray(points)
+    r_vecs = points[:, np.newaxis, :] - coords[np.newaxis, :, :]
+    r = np.linalg.norm(r_vecs, axis=-1)
+    safe_r = np.where(r < 1e-12, np.inf, r)
+    p_dot_r = np.einsum("ij,mij->mi", dipoles, r_vecs)
+    return np.sum(p_dot_r / safe_r ** 3, axis=1)
+
+
+def dipole_field(dipoles, coords, points):
+    """
+    Compute electric field from atomic dipoles.
+
+    E(r) = sum_i [ 3(p_i . r̂)r̂ - p_i ] / |r - r_i|^3
+
+    Parameters
+    ----------
+    dipoles : np.ndarray, shape (N, 3)
+    coords : np.ndarray, shape (N, 3)
+    points : np.ndarray, shape (M, 3)
+
+    Returns
+    -------
+    field : np.ndarray, shape (M, 3)
+    """
+    dipoles = np.asarray(dipoles)
+    coords = np.asarray(coords)
+    points = np.asarray(points)
+    r_vecs = points[:, np.newaxis, :] - coords[np.newaxis, :, :]
+    r = np.linalg.norm(r_vecs, axis=-1)
+    safe_r = np.where(r < 1e-12, np.inf, r)
+    r_hat = r_vecs / safe_r[:, :, np.newaxis]
+    p_dot_rhat = np.einsum("ij,mij->mi", dipoles, r_hat)
+    term = (3 * p_dot_rhat[:, :, np.newaxis] * r_hat
+            - dipoles[np.newaxis, :, :])
+    return np.sum(term / safe_r[:, :, np.newaxis] ** 3, axis=1)
+
+
+def quadrupole_potential(quadrupoles, coords, points):
+    """
+    Compute electrostatic potential from atomic quadrupoles (traceless).
+
+    V(r) = sum_i [ Q_i:rr ] / (2 * |r - r_i|^5)
+
+    Parameters
+    ----------
+    quadrupoles : np.ndarray, shape (N, 3, 3)
+        Traceless quadrupole moment tensors in atomic units.
+    coords : np.ndarray, shape (N, 3)
+    points : np.ndarray, shape (M, 3)
+
+    Returns
+    -------
+    potential : np.ndarray, shape (M,)
+    """
+    quadrupoles = np.asarray(quadrupoles)
+    coords = np.asarray(coords)
+    points = np.asarray(points)
+    r_vecs = points[:, np.newaxis, :] - coords[np.newaxis, :, :]
+    r = np.linalg.norm(r_vecs, axis=-1)
+    safe_r = np.where(r < 1e-12, np.inf, r)
+    Qrr = np.einsum("nab,mna,mnb->mn", quadrupoles, r_vecs, r_vecs)
+    return np.sum(Qrr / (2 * safe_r ** 5), axis=1)
+
+
+def quadrupole_field(quadrupoles, coords, points):
+    """
+    Compute electric field from atomic quadrupoles (traceless).
+
+    E(r) = sum_i [ 5*(Q_i:rr)*r / (2r^7) - (Q_i.r) / r^5 ]
+
+    Parameters
+    ----------
+    quadrupoles : np.ndarray, shape (N, 3, 3)
+    coords : np.ndarray, shape (N, 3)
+    points : np.ndarray, shape (M, 3)
+
+    Returns
+    -------
+    field : np.ndarray, shape (M, 3)
+    """
+    quadrupoles = np.asarray(quadrupoles)
+    coords = np.asarray(coords)
+    points = np.asarray(points)
+    r_vecs = points[:, np.newaxis, :] - coords[np.newaxis, :, :]
+    r = np.linalg.norm(r_vecs, axis=-1)
+    safe_r = np.where(r < 1e-12, np.inf, r)
+    Qrr = np.einsum("nab,mna,mnb->mn", quadrupoles, r_vecs, r_vecs)
+    Qr = np.einsum("nab,mnb->mna", quadrupoles, r_vecs)
+    term1 = (5 * Qrr[:, :, np.newaxis] * r_vecs
+             / (2 * safe_r[:, :, np.newaxis] ** 7))
+    term2 = Qr / safe_r[:, :, np.newaxis] ** 5
+    return np.sum(term1 - term2, axis=1)
+
+
+def total_potential(coords, points, charges=None, dipoles=None, quadrupoles=None):
+    """
+    Compute total electrostatic potential from all multipole contributions.
+
+    Parameters
+    ----------
+    coords : np.ndarray, shape (N, 3)
+    points : np.ndarray, shape (M, 3)
+    charges : np.ndarray, shape (N,), optional
+    dipoles : np.ndarray, shape (N, 3), optional
+    quadrupoles : np.ndarray, shape (N, 3, 3), optional
+
+    Returns
+    -------
+    potential : np.ndarray, shape (M,)
+    """
+    potential = np.zeros(points.shape[0])
+    if charges is not None:
+        potential += monopole_potential(charges, coords, points)
+    if dipoles is not None:
+        potential += dipole_potential(dipoles, coords, points)
+    if quadrupoles is not None:
+        potential += quadrupole_potential(quadrupoles, coords, points)
+    return potential
+
+
+def total_field(coords, points, charges=None, dipoles=None, quadrupoles=None):
+    """
+    Compute total electric field from all multipole contributions.
+
+    Parameters
+    ----------
+    coords : np.ndarray, shape (N, 3)
+    points : np.ndarray, shape (M, 3)
+    charges : np.ndarray, shape (N,), optional
+    dipoles : np.ndarray, shape (N, 3), optional
+    quadrupoles : np.ndarray, shape (N, 3, 3), optional
+
+    Returns
+    -------
+    field : np.ndarray, shape (M, 3)
+    """
+    field = np.zeros((points.shape[0], 3))
+    if charges is not None:
+        field += monopole_field(charges, coords, points)
+    if dipoles is not None:
+        field += dipole_field(dipoles, coords, points)
+    if quadrupoles is not None:
+        field += quadrupole_field(quadrupoles, coords, points)
+    return field
\ No newline at end of file
diff --git a/tests/test_multipole.py b/tests/test_multipole.py
new file mode 100644
index 0000000..2d5e79a
--- /dev/null
+++ b/tests/test_multipole.py
@@ -0,0 +1,333 @@
+"""
+Tests for electrostatic multipole potential and field functions.
+All tests validate against analytical solutions derived from first principles.
+"""
+
+import numpy as np
+import pytest
+from numpy.testing import assert_allclose
+from ffprime.electrostatics.multipole import (
+    monopole_potential,
+    monopole_field,
+    dipole_potential,
+    dipole_field,
+    quadrupole_potential,
+    quadrupole_field,
+    total_potential,
+    total_field,
+)
+
+
+# ─────────────────────────────────────────────
+# MONOPOLE TESTS
+# ─────────────────────────────────────────────
+
+class TestMonopolePotential:
+
+    def test_single_unit_charge(self):
+        """V = q/r → at distance 1.0, V = 1.0"""
+        charges = np.array([1.0])
+        coords = np.array([[0.0, 0.0, 0.0]])
+        points = np.array([[1.0, 0.0, 0.0]])
+        result = monopole_potential(charges, coords, points)
+        assert_allclose(result, [1.0], rtol=1e-10)
+
+    def test_inverse_distance_scaling(self):
+        """V scales as 1/r"""
+        charges = np.array([1.0])
+        coords = np.array([[0.0, 0.0, 0.0]])
+        points = np.array([[2.0, 0.0, 0.0],
+                           [4.0, 0.0, 0.0]])
+        result = monopole_potential(charges, coords, points)
+        assert_allclose(result[0] / result[1], 2.0, rtol=1e-10)
+
+    def test_negative_charge(self):
+        """Negative charge gives negative potential"""
+        charges = np.array([-1.0])
+        coords = np.array([[0.0, 0.0, 0.0]])
+        points = np.array([[1.0, 0.0, 0.0]])
+        result = monopole_potential(charges, coords, points)
+        assert_allclose(result, [-1.0], rtol=1e-10)
+
+    def test_superposition(self):
+        """Two equal charges: potential is sum of individual contributions"""
+        charges = np.array([1.0, 1.0])
+        coords = np.array([[1.0, 0.0, 0.0],
+                           [-1.0, 0.0, 0.0]])
+        points = np.array([[0.0, 2.0, 0.0]])
+        result = monopole_potential(charges, coords, points)
+        expected = 2.0 / np.sqrt(5.0)
+        assert_allclose(result, [expected], rtol=1e-10)
+
+    def test_singularity_safe(self):
+        """Potential at exact atom location returns finite value"""
+        charges = np.array([1.0])
+        coords = np.array([[0.0, 0.0, 0.0]])
+        points = np.array([[0.0, 0.0, 0.0]])
+        result = monopole_potential(charges, coords, points)
+        assert np.isfinite(result).all()
+
+    def test_multiple_field_points(self):
+        """Vectorized evaluation over multiple points"""
+        charges = np.array([1.0])
+        coords = np.array([[0.0, 0.0, 0.0]])
+        points = np.array([[1.0, 0.0, 0.0],
+                           [2.0, 0.0, 0.0],
+                           [3.0, 0.0, 0.0]])
+        result = monopole_potential(charges, coords, points)
+        assert_allclose(result, [1.0, 0.5, 1/3], rtol=1e-10)
+
+
+class TestMonopoleField:
+
+    def test_unit_charge_along_x(self):
+        """E = q*r/r^3 → at (1,0,0), E = (1,0,0)"""
+        charges = np.array([1.0])
+        coords = np.array([[0.0, 0.0, 0.0]])
+        points = np.array([[1.0, 0.0, 0.0]])
+        result = monopole_field(charges, coords, points)
+        assert_allclose(result, [[1.0, 0.0, 0.0]], rtol=1e-10)
+
+    def test_inverse_square_scaling(self):
+        """E magnitude scales as 1/r^2"""
+        charges = np.array([1.0])
+        coords = np.array([[0.0, 0.0, 0.0]])
+        points = np.array([[2.0, 0.0, 0.0],
+                           [4.0, 0.0, 0.0]])
+        result = monopole_field(charges, coords, points)
+        ratio = np.linalg.norm(result[0]) / np.linalg.norm(result[1])
+        assert_allclose(ratio, 4.0, rtol=1e-10)
+
+    def test_field_is_radial(self):
+        """Field points radially away from positive charge"""
+        charges = np.array([1.0])
+        coords = np.array([[0.0, 0.0, 0.0]])
+        points = np.array([[1.0, 1.0, 1.0]])
+        result = monopole_field(charges, coords, points)
+        r = np.array([1.0, 1.0, 1.0])
+        r_hat = r / np.linalg.norm(r)
+        result_hat = result[0] / np.linalg.norm(result[0])
+        assert_allclose(result_hat, r_hat, rtol=1e-10)
+
+    def test_singularity_safe(self):
+        """Field at exact atom location returns finite value"""
+        charges = np.array([1.0])
+        coords = np.array([[0.0, 0.0, 0.0]])
+        points = np.array([[0.0, 0.0, 0.0]])
+        result = monopole_field(charges, coords, points)
+        assert np.isfinite(result).all()
+
+
+# ─────────────────────────────────────────────
+# DIPOLE TESTS
+# ─────────────────────────────────────────────
+
+class TestDipolePotential:
+
+    def test_dipole_along_axis(self):
+        """
+        p=(1,0,0) at origin, point at (2,0,0):
+        V = (p.r)/r^3 = 2/8 = 0.25
+        """
+        dipoles = np.array([[1.0, 0.0, 0.0]])
+        coords = np.array([[0.0, 0.0, 0.0]])
+        points = np.array([[2.0, 0.0, 0.0]])
+        result = dipole_potential(dipoles, coords, points)
+        assert_allclose(result, [0.25], rtol=1e-10)
+
+    def test_dipole_perpendicular_zero(self):
+        """
+        p=(1,0,0), point at (0,2,0):
+        V = (p.r)/r^3 = 0 (perpendicular)
+        """
+        dipoles = np.array([[1.0, 0.0, 0.0]])
+        coords = np.array([[0.0, 0.0, 0.0]])
+        points = np.array([[0.0, 2.0, 0.0]])
+        result = dipole_potential(dipoles, coords, points)
+        assert_allclose(result, [0.0], atol=1e-12)
+
+    def test_dipole_antisymmetry(self):
+        """V at +r and -r are equal and opposite"""
+        dipoles = np.array([[1.0, 0.0, 0.0]])
+        coords = np.array([[0.0, 0.0, 0.0]])
+        points = np.array([[2.0, 0.0, 0.0],
+                           [-2.0, 0.0, 0.0]])
+        result = dipole_potential(dipoles, coords, points)
+        assert_allclose(result[0], -result[1], rtol=1e-10)
+
+    def test_dipole_scales_as_inverse_r_squared(self):
+        """V_dipole scales as 1/r^2"""
+        dipoles = np.array([[1.0, 0.0, 0.0]])
+        coords = np.array([[0.0, 0.0, 0.0]])
+        points = np.array([[1.0, 0.0, 0.0],
+                           [2.0, 0.0, 0.0]])
+        result = dipole_potential(dipoles, coords, points)
+        assert_allclose(result[0] / result[1], 4.0, rtol=1e-10)
+
+
+class TestDipoleField:
+
+    def test_dipole_field_along_axis(self):
+        """
+        p=(0,0,1) at origin, point at (0,0,2):
+        E = [3(p.r̂)r̂ - p] / r^3
+          = [3*(0,0,1) - (0,0,1)] / 8
+          = (0, 0, 0.25)
+        """
+        dipoles = np.array([[0.0, 0.0, 1.0]])
+        coords = np.array([[0.0, 0.0, 0.0]])
+        points = np.array([[0.0, 0.0, 2.0]])
+        result = dipole_field(dipoles, coords, points)
+        assert_allclose(result, [[0.0, 0.0, 0.25]], rtol=1e-10)
+
+    def test_dipole_field_perpendicular(self):
+        """
+        p=(0,0,1) at origin, point at (2,0,0):
+        r̂=(1,0,0), p.r̂=0
+        E = [0 - (0,0,1)] / 8 = (0, 0, -0.125)
+        """
+        dipoles = np.array([[0.0, 0.0, 1.0]])
+        coords = np.array([[0.0, 0.0, 0.0]])
+        points = np.array([[2.0, 0.0, 0.0]])
+        result = dipole_field(dipoles, coords, points)
+        assert_allclose(result, [[0.0, 0.0, -0.125]], rtol=1e-10)
+
+    def test_singularity_safe(self):
+        """Field at exact atom location returns finite value"""
+        dipoles = np.array([[1.0, 0.0, 0.0]])
+        coords = np.array([[0.0, 0.0, 0.0]])
+        points = np.array([[0.0, 0.0, 0.0]])
+        result = dipole_field(dipoles, coords, points)
+        assert np.isfinite(result).all()
+
+
+# ─────────────────────────────────────────────
+# QUADRUPOLE TESTS
+# ─────────────────────────────────────────────
+
+class TestQuadrupolePotential:
+
+    def test_quadrupole_along_z(self):
+        """
+        Q = diag(-1,-1,2) (traceless), point at (0,0,2):
+        Qrr = Q_ab r_a r_b = 2*4 = 8
+        V = Qrr / (2*r^5) = 8 / (2*32) = 0.125
+        """
+        Q = np.array([[[-1.0, 0.0, 0.0],
+                       [0.0, -1.0, 0.0],
+                       [0.0,  0.0, 2.0]]])
+        coords = np.array([[0.0, 0.0, 0.0]])
+        points = np.array([[0.0, 0.0, 2.0]])
+        result = quadrupole_potential(Q, coords, points)
+        assert_allclose(result, [0.125], rtol=1e-10)
+
+    def test_traceless_symmetry(self):
+        """Trace of Q should not contribute — pure traceless result"""
+        Q = np.array([[[2.0, 0.0, 0.0],
+                       [0.0, -1.0, 0.0],
+                       [0.0,  0.0, -1.0]]])
+        coords = np.array([[0.0, 0.0, 0.0]])
+        points = np.array([[2.0, 0.0, 0.0]])
+        result = quadrupole_potential(Q, coords, points)
+        assert np.isfinite(result).all()
+
+    def test_singularity_safe(self):
+        """Potential at exact atom location returns finite value"""
+        Q = np.array([[[1.0, 0.0, 0.0],
+                       [0.0, -0.5, 0.0],
+                       [0.0, 0.0, -0.5]]])
+        coords = np.array([[0.0, 0.0, 0.0]])
+        points = np.array([[0.0, 0.0, 0.0]])
+        result = quadrupole_potential(Q, coords, points)
+        assert np.isfinite(result).all()
+
+
+class TestQuadrupoleField:
+
+    def test_quadrupole_field_finite(self):
+        """Field from quadrupole is finite at valid points"""
+        Q = np.array([[[-1.0, 0.0, 0.0],
+                       [0.0, -1.0, 0.0],
+                       [0.0,  0.0, 2.0]]])
+        coords = np.array([[0.0, 0.0, 0.0]])
+        points = np.array([[0.0, 0.0, 2.0],
+                           [1.0, 0.0, 0.0],
+                           [0.0, 1.0, 1.0]])
+        result = quadrupole_field(Q, coords, points)
+        assert np.isfinite(result).all()
+
+    def test_singularity_safe(self):
+        """Field at exact atom location returns finite value"""
+        Q = np.array([[[1.0, 0.0, 0.0],
+                       [0.0, -0.5, 0.0],
+                       [0.0, 0.0, -0.5]]])
+        coords = np.array([[0.0, 0.0, 0.0]])
+        points = np.array([[0.0, 0.0, 0.0]])
+        result = quadrupole_field(Q, coords, points)
+        assert np.isfinite(result).all()
+
+
+# ─────────────────────────────────────────────
+# TOTAL POTENTIAL AND FIELD TESTS
+# ─────────────────────────────────────────────
+
+class TestTotalPotential:
+
+    def test_monopole_only(self):
+        """total_potential with only charges matches monopole_potential"""
+        charges = np.array([1.0])
+        coords = np.array([[0.0, 0.0, 0.0]])
+        points = np.array([[1.0, 0.0, 0.0]])
+        result = total_potential(coords, points, charges=charges)
+        expected = monopole_potential(charges, coords, points)
+        assert_allclose(result, expected, rtol=1e-10)
+
+    def test_superposition_all_terms(self):
+        """Total potential is sum of all individual contributions"""
+        charges = np.array([1.0])
+        dipoles = np.array([[1.0, 0.0, 0.0]])
+        Q = np.array([[[-1.0, 0.0, 0.0],
+                       [0.0, -1.0, 0.0],
+                       [0.0,  0.0, 2.0]]])
+        coords = np.array([[0.0, 0.0, 0.0]])
+        points = np.array([[2.0, 1.0, 0.0]])
+
+        total = total_potential(coords, points,
+                                charges=charges,
+                                dipoles=dipoles,
+                                quadrupoles=Q)
+        expected = (monopole_potential(charges, coords, points) +
+                    dipole_potential(dipoles, coords, points) +
+                    quadrupole_potential(Q, coords, points))
+        assert_allclose(total, expected, rtol=1e-10)
+
+
+class TestTotalField:
+
+    def test_monopole_only(self):
+        """total_field with only charges matches monopole_field"""
+        charges = np.array([1.0])
+        coords = np.array([[0.0, 0.0, 0.0]])
+        points = np.array([[1.0, 0.0, 0.0]])
+        result = total_field(coords, points, charges=charges)
+        expected = monopole_field(charges, coords, points)
+        assert_allclose(result, expected, rtol=1e-10)
+
+    def test_superposition_all_terms(self):
+        """Total field is sum of all individual contributions"""
+        charges = np.array([1.0])
+        dipoles = np.array([[1.0, 0.0, 0.0]])
+        Q = np.array([[[-1.0, 0.0, 0.0],
+                       [0.0, -1.0, 0.0],
+                       [0.0,  0.0, 2.0]]])
+        coords = np.array([[0.0, 0.0, 0.0]])
+        points = np.array([[2.0, 1.0, 0.0]])
+
+        total = total_field(coords, points,
+                            charges=charges,
+                            dipoles=dipoles,
+                            quadrupoles=Q)
+        expected = (monopole_field(charges, coords, points) +
+                    dipole_field(dipoles, coords, points) +
+                    quadrupole_field(Q, coords, points))
+        assert_allclose(total, expected, rtol=1e-10)
\ No newline at end of file