From 5c512e73469e0235a0f97f41ae869937eaa1a753 Mon Sep 17 00:00:00 2001
From: RAHUL BALAJI <rahul.balaji13289@gmail.com>
Date: Sun, 8 Mar 2026 17:52:48 +0530
Subject: [PATCH] Add dipole electrostatic potential and field functions

Implements dipole_potential() and dipole_field() following the same vectorized NumPy pattern as the monopole functions. Includes proper handling of edge cases (r->0) and comprehensive docstrings with physics formulas.

Builds on PR #2 to extend electrostatics module toward full multipole support.
---
 ffprime/electrostatics/multipole.py | 55 +++++++++++++++++++++++++++++
 1 file changed, 55 insertions(+)
 create mode 100644 ffprime/electrostatics/multipole.py

diff --git a/ffprime/electrostatics/multipole.py b/ffprime/electrostatics/multipole.py
new file mode 100644
index 0000000..6786fac
--- /dev/null
+++ b/ffprime/electrostatics/multipole.py
@@ -0,0 +1,55 @@
+import numpy as np
+
+
+def dipole_potential(p: np.ndarray, r_vec: np.ndarray) -> float:
+    """Compute electrostatic potential from dipole moment.
+    
+    Physics: V = (p · r̂) / r² in atomic units
+    
+    Parameters
+    ----------
+    p : np.ndarray
+        Dipole moment vector (3,)
+    r_vec : np.ndarray
+        Position vector(s) (N, 3) or (3,)
+        
+    Returns
+    -------
+    float or np.ndarray
+        Potential at given position(s)
+    """
+    r = np.linalg.norm(r_vec, axis=-1 if r_vec.ndim > 1 else 0)
+    
+    if np.any(r < 1e-12):
+        return 0.0 if np.isscalar(r) else np.where(r < 1e-12, 0.0, 
+                                                    np.dot(r_vec, p) / r**3)
+    
+    return np.dot(r_vec, p) / r**3
+
+
+def dipole_field(p: np.ndarray, r_vec: np.ndarray) -> np.ndarray:
+    """Compute electric field from dipole moment.
+    
+    Physics: E = (1/r³)[3(p·r̂)r̂ - p] in atomic units
+    
+    Parameters
+    ----------
+    p : np.ndarray
+        Dipole moment vector (3,)
+    r_vec : np.ndarray
+        Position vector(s) (N, 3) or (3,)
+        
+    Returns
+    -------
+    np.ndarray
+        Electric field at given position(s), shape (N, 3) or (3,)
+    """
+    r = np.linalg.norm(r_vec, axis=-1, keepdims=True if r_vec.ndim > 1 else False)
+    
+    if np.any(r < 1e-12):
+        return np.zeros_like(r_vec)
+    
+    r_hat = r_vec / r
+    p_dot_r = np.dot(r_vec, p) if r_vec.ndim == 1 else np.sum(r_vec * p, axis=1, keepdims=True)
+    
+    return (3 * p_dot_r * r_hat - p) / r**3