diff --git a/ffprime/electrostatics/spherical.py b/ffprime/electrostatics/spherical.py new file mode 100644 index 0000000..fd7ee60 --- /dev/null +++ b/ffprime/electrostatics/spherical.py @@ -0,0 +1,112 @@ +import numpy as np + + +def dipole_cartesian_to_spherical(p: np.ndarray) -> np.ndarray: + """ + Convert a Cartesian dipole moment to real spherical (solid harmonic) form. + + Parameters + ---------- + p : np.ndarray, shape (3,) + Cartesian dipole moment [px, py, pz] in atomic units (e * bohr). + + Returns + ------- + np.ndarray, shape (3,) + Spherical components [Q_10, Q_11c, Q_11s]. + + Notes + ----- + Convention follows Stone, A.J. "Theory of Intermolecular Forces", 2nd ed. + Q_10 = pz, Q_11c = px, Q_11s = py + """ + px, py, pz = p + return np.array([pz, px, py]) + + +def dipole_spherical_to_cartesian(q: np.ndarray) -> np.ndarray: + """ + Convert real spherical dipole components back to Cartesian form. + + Parameters + ---------- + q : np.ndarray, shape (3,) + Spherical components [Q_10, Q_11c, Q_11s]. + + Returns + ------- + np.ndarray, shape (3,) + Cartesian dipole moment [px, py, pz] in atomic units (e * bohr). + """ + Q_10, Q_11c, Q_11s = q + return np.array([Q_11c, Q_11s, Q_10]) + + +def quadrupole_cartesian_to_spherical(theta: np.ndarray) -> np.ndarray: + """ + Convert a traceless Cartesian quadrupole tensor to real spherical form. + + Parameters + ---------- + theta : np.ndarray, shape (3, 3) + Traceless symmetric quadrupole tensor in atomic units (e * bohr^2). + Must satisfy np.trace(theta) == 0. + + Returns + ------- + np.ndarray, shape (5,) + Spherical components [Q_20, Q_21c, Q_21s, Q_22c, Q_22s]. + + Raises + ------ + ValueError + If the input tensor is not traceless within numerical tolerance. + + Notes + ----- + Convention follows Stone, A.J. "Theory of Intermolecular Forces", 2nd ed. + Q_20 = Θzz + Q_21c = Θxz + Q_21s = Θyz + Q_22c = (Θxx - Θyy) / 2 + Q_22s = Θxy + """ + if not np.isclose(np.trace(theta), 0.0, atol=1e-10): + raise ValueError( + f"Quadrupole tensor must be traceless, got trace = {np.trace(theta):.6e}" + ) + + Q_20 = theta[2, 2] + Q_21c = theta[0, 2] + Q_21s = theta[1, 2] + Q_22c = (theta[0, 0] - theta[1, 1]) / 2.0 + Q_22s = theta[0, 1] + + return np.array([Q_20, Q_21c, Q_21s, Q_22c, Q_22s]) + + +def quadrupole_spherical_to_cartesian(q: np.ndarray) -> np.ndarray: + """ + Convert real spherical quadrupole components back to Cartesian tensor form. + + Parameters + ---------- + q : np.ndarray, shape (5,) + Spherical components [Q_20, Q_21c, Q_21s, Q_22c, Q_22s]. + + Returns + ------- + np.ndarray, shape (3, 3) + Traceless symmetric quadrupole tensor in atomic units (e * bohr^2). + """ + Q_20, Q_21c, Q_21s, Q_22c, Q_22s = q + + theta = np.zeros((3, 3)) + theta[2, 2] = Q_20 + theta[0, 0] = Q_22c - Q_20 / 2.0 + theta[1, 1] = -Q_22c - Q_20 / 2.0 + theta[0, 2] = theta[2, 0] = Q_21c + theta[1, 2] = theta[2, 1] = Q_21s + theta[0, 1] = theta[1, 0] = Q_22s + + return theta \ No newline at end of file diff --git a/tests/test_spherical.py b/tests/test_spherical.py new file mode 100644 index 0000000..18cb98e --- /dev/null +++ b/tests/test_spherical.py @@ -0,0 +1,57 @@ +import numpy as np +import pytest +from ffprime.electrostatics.spherical import ( + dipole_cartesian_to_spherical, + dipole_spherical_to_cartesian, + quadrupole_cartesian_to_spherical, + quadrupole_spherical_to_cartesian, +) + + +def test_dipole_z_aligned(): + """Pure z-dipole maps entirely to Q_10, others zero.""" + p = np.array([0.0, 0.0, 1.0]) + sph = dipole_cartesian_to_spherical(p) + assert np.isclose(sph[0], 1.0) + assert np.isclose(sph[1], 0.0) + assert np.isclose(sph[2], 0.0) + + +def test_dipole_roundtrip(): + """Cartesian -> spherical -> Cartesian recovers original.""" + p = np.array([1.5, -2.0, 3.0]) + assert np.allclose( + dipole_spherical_to_cartesian(dipole_cartesian_to_spherical(p)), p + ) + + +def test_quadrupole_traceless_check(): + """Non-traceless tensor raises ValueError.""" + bad_tensor = np.eye(3) + with pytest.raises(ValueError, match="traceless"): + quadrupole_cartesian_to_spherical(bad_tensor) + + +def test_quadrupole_z_axial(): + """Axially symmetric quadrupole along z has only Q_20 nonzero.""" + theta = np.array([ + [-0.5, 0.0, 0.0], + [ 0.0, -0.5, 0.0], + [ 0.0, 0.0, 1.0] + ]) + sph = quadrupole_cartesian_to_spherical(theta) + assert np.isclose(sph[0], 1.0) + assert np.allclose(sph[1:], 0.0) + + +def test_quadrupole_roundtrip(): + """Cartesian -> spherical -> Cartesian recovers original tensor.""" + theta = np.array([ + [-0.5, 0.3, 0.1], + [ 0.3, -0.5, 0.2], + [ 0.1, 0.2, 1.0] + ]) + recovered = quadrupole_spherical_to_cartesian( + quadrupole_cartesian_to_spherical(theta) + ) + assert np.allclose(recovered, theta) \ No newline at end of file