feat: Allow custom callable potentials in exchange kernel functions

Author: skilledwolf

README.md

@@ -37,7 +37,7 @@ The package performs calculations in dimensionless units where lengths are scale
 - **Coulomb interaction**: The code assumes a potential of the form $V(q) = \kappa \frac{2\pi e^2}{q \ell_B}$ (in effective dimensionless form).
   - If you set `kappa = 1.0`, the resulting exchange kernels are in units of the Coulomb energy scale $E_C = e^2 / (\epsilon \ell_B)$.
   - To express results in units of the cyclotron energy $\hbar \omega_c$, set $\kappa = E_C / (\hbar \omega_c) = (e^2/\epsilon \ell_B) / (\hbar \omega_c)$.
-- **General potential**: For a general $V(q)$, `V_of_q` should return values in your desired energy units. The integration measure $d^2q/(2\pi)^2$ introduces a factor of $1/\ell_B^2$, so ensure your potential scaling is consistent.
+- **Custom potential**: Provide a callable `potential(q)` that returns values in your desired energy units. The integration measure $d^2q/(2\pi)^2$ introduces a factor of $1/\ell_B^2$, so ensure your potential scaling is consistent.
 
 ## Installation
 
@@ -74,7 +74,7 @@ print("F shape:", F.shape)
 print("X shape:", X.shape)
 ```
 
-To use a user-provided Coulomb interaction, pass a callable `V_of_q` and select the `"general"` potential:
+To use a user-provided interaction, pass a callable directly as `potential`:
 
 ```python
 def V_coulomb(q, kappa=1.0):
@@ -86,8 +86,7 @@ X_coulomb = get_exchange_kernels(
     thetas,
     nmax,
     method="gausslegendre",
-    potential="general",
-    V_of_q=lambda q: V_coulomb(q, kappa=1.0),
+    potential=lambda q: V_coulomb(q, kappa=1.0),
 )
 ```
 

src/quantumhall_matrixelements/exchange_gausslag.py

@@ -61,9 +61,8 @@ def get_exchange_kernels_GaussLag(
     G_angles,
     nmax: int,
     *,
-    potential: str = "coulomb",
+    potential: str | callable = "coulomb",
     kappa: float = 1.0,
-    V_of_q=None,
     nquad: int = 200,
     ell: float = 1.0,
 ) -> "ComplexArray":
@@ -76,13 +75,11 @@ def get_exchange_kernels_GaussLag(
     nmax :
         Number of Landau levels.
     potential :
-        Either ``'coulomb'`` (default) or ``'general'``. In the latter case
-        a callable ``V_of_q(q)`` must be provided.
+        ``'coulomb'`` (default), ``'constant'``, or a callable ``V(q)`` giving
+        the interaction in 1/ℓ units.
     kappa :
         Interaction strength prefactor. For Coulomb this corresponds to
         :math:`\\kappa = e^2/(\\varepsilon\\ell_B)/\\hbar\\omega_c`.
-    V_of_q :
-        Callable ``V_of_q(q) -> V(q)`` used when ``potential='general'``.
     nquad :
         Number of Gauss–Laguerre quadrature points.
     ell :
@@ -99,6 +96,21 @@ def get_exchange_kernels_GaussLag(
         raise ValueError("G_magnitudes and G_angles must have the same shape.")
     nG = G_magnitudes.size
 
+    # Resolve potential
+    if callable(potential):
+        pot_kind = "callable"
+        pot_fn = potential
+    else:
+        pot_kind = str(potential).strip().lower()
+        pot_fn = None
+
+    if pot_kind in {"coulomb", "constant"}:
+        pass
+    elif pot_kind == "callable":
+        pass
+    else:
+        raise ValueError("potential must be 'coulomb', 'constant', or a callable V(q).")
+
     Gscaled = G_magnitudes * float(ell)
     Xs = np.zeros((nG, nmax, nmax, nmax, nmax), dtype=np.complex128)
 
@@ -131,16 +143,12 @@ def get_exchange_kernels_GaussLag(
                         phase_factor = (1j) ** (d1 - d2)
                         pref = (kappa * C / np.sqrt(2.0)) * phase_factor
                     else:
-                        if not callable(V_of_q):
-                            raise ValueError(
-                                "For potential='general', provide V_of_q: callable(q)->V(q)."
-                            )
                         alpha = 0.5 * (d1 + d2)
                         z, w = _lag_nodes_weights(nquad, alpha)
                         L1 = _laguerre_on_grid(p, d1, alpha, nquad, z)
                         L2 = _laguerre_on_grid(q, d2, alpha, nquad, z)
                         qvals = np.sqrt(2.0 * z) / float(ell)
-                        Veff = V_of_q(qvals) / (2.0 * np.pi * float(ell) ** 2)
+                        Veff = pot_fn(qvals) / (2.0 * np.pi * float(ell) ** 2)
                         W = w * L1 * L2 * Veff
                         key = (absN, float(alpha))
                         J_abs = J_cache.get(key)
@@ -160,4 +168,3 @@ def get_exchange_kernels_GaussLag(
 
 
 __all__ = ["get_exchange_kernels_GaussLag"]
-

src/quantumhall_matrixelements/exchange_hankel.py

@@ -75,14 +75,28 @@ def get_exchange_kernels_hankel(
     G_magnitudes: "RealArray",
     G_angles: "RealArray",
     nmax: int,
-    **kwargs,
+    *,
+    potential: str | callable = "coulomb",
+    kappa: float = 1.0,
 ) -> "ComplexArray":
     """Compute X_{n1,m1,n2,m2}(G) via Hankel transforms (κ=1 convention).
 
     This backend parametrizes the radial integral via Hankel transforms with
     robust Laguerre-based normalization and explicit control over the Bessel
     order. It is numerically more intensive than the Gauss–Laguerre backend
     but can be useful for cross-checks or alternative potentials.
+
+    Parameters
+    ----------
+    G_magnitudes, G_angles :
+        Arrays describing |G| and polar angle θ_G (same shape, no broadcasting).
+    nmax :
+        Number of Landau levels.
+    potential :
+        ``'coulomb'`` (default), ``'constant'``, or a callable ``V(q)`` giving
+        the interaction in 1/ℓ units.
+    kappa :
+        Prefactor for Coulomb/constant cases.
     """
     G_magnitudes = np.asarray(G_magnitudes, dtype=float)
     G_angles = np.asarray(G_angles, dtype=float)
@@ -94,6 +108,18 @@ def get_exchange_kernels_hankel(
     k_unique, inv_idx = np.unique(G_magnitudes, return_inverse=True)
     nG = G_magnitudes.size
 
+    # Resolve potential
+    if callable(potential):
+        pot_callable = potential
+        pot_key = ("callable", id(potential))
+    else:
+        pot_name = str(potential).strip().lower()
+        if pot_name in {"coulomb", "constant"}:
+            pot_callable = pot_name
+            pot_key = (pot_name, float(kappa))
+        else:
+            raise ValueError("potential must be 'coulomb', 'constant', or callable V(q)")
+
     # Precompute angular phase vectors for all possible N values
     # N = (n1 - m1) - (m2 - n2) ∈ [Nmin, Nmax]
     N_min = -2 * (nmax - 1)
@@ -115,8 +141,8 @@ def get_exchange_kernels_hankel(
         np.arange(nmax)[None, :].repeat(nmax, axis=0),
     )
 
-    # Cache for radial Hankel transforms keyed by (d1,N1,d2,N2,absN)
-    radial_cache: dict[tuple[int, int, int, int, int], np.ndarray] = {}
+    # Cache for radial Hankel transforms keyed by (d1,N1,d2,N2,absN,pot_key)
+    radial_cache: dict[tuple[int, int, int, int, int, tuple], np.ndarray] = {}
 
     # Output
     Xs = np.zeros((nG, nmax, nmax, nmax, nmax), dtype=np.complex128)
@@ -127,13 +153,15 @@ def get_radial_block(n1: int, m1: int, n2: int, m2: int, absN: int) -> np.ndarra
         d2 = int(d_mat[n2, m2])
         N1 = int(Nmin_mat[n1, m1])
         N2 = int(Nmin_mat[n2, m2])
-        key = (d1, N1, d2, N2, int(absN))
+        key = (d1, N1, d2, N2, int(absN), pot_key)
         arr = radial_cache.get(key)
         if arr is not None:
             return arr
 
         def integrand(q):
-            return _radial_exchange_integrand_rgamma(q, n1, m1, n2, m2)
+            return _radial_exchange_integrand_rgamma(
+                q, n1, m1, n2, m2, potential=pot_callable, kappa=kappa
+            )
 
         ht = _get_hankel_transformer(absN)
         # Compute only on unique radii, then cache
@@ -162,4 +190,3 @@ def integrand(q):
 
 
 __all__ = ["get_exchange_kernels_hankel"]
-

src/quantumhall_matrixelements/exchange_legendre.py

@@ -57,9 +57,8 @@ def get_exchange_kernels_GaussLegendre(
     G_angles,
     nmax: int,
     *,
-    potential: str = "coulomb",
+    potential: str | callable = "coulomb",
     kappa: float = 1.0,
-    V_of_q=None,
     nquad: int = 1000,
     scale: float = 0.5,
     ell: float = 1.0,
@@ -78,11 +77,9 @@ def get_exchange_kernels_GaussLegendre(
     nmax :
         Number of Landau levels.
     potential :
-        Either ``'coulomb'`` (default) or ``'general'``.
+        ``'coulomb'`` (default) or a callable ``V(q)`` returning the interaction.
     kappa :
         Interaction strength prefactor.
-    V_of_q :
-        Callable ``V_of_q(q) -> V(q)`` used when ``potential='general'``.
     nquad :
         Number of quadrature points (default 1000).
     scale :
@@ -104,6 +101,21 @@ def get_exchange_kernels_GaussLegendre(
     Gscaled = G_magnitudes * float(ell)
     Xs = np.zeros((nG, nmax, nmax, nmax, nmax), dtype=np.complex128)
 
+    # Resolve potential
+    if callable(potential):
+        pot_kind = "callable"
+        pot_fn = potential
+    else:
+        pot_kind = str(potential).strip().lower()
+        pot_fn = None
+
+    if pot_kind == "coulomb":
+        pass
+    elif pot_kind == "callable":
+        pass
+    else:
+        raise ValueError("potential must be 'coulomb' or a callable V(q).")
+
     # Get mapped grid
     z, w = _legendre_nodes_weights_mapped(nquad, scale)
 
@@ -154,16 +166,13 @@ def get_exchange_kernels_GaussLegendre(
                         pref = (kappa * C / np.sqrt(2.0)) * phase_factor
 
                     else:
-                        # General potential
-                        if not callable(V_of_q):
-                            raise ValueError("Provide V_of_q for potential='general'")
-                        
+                        # General/callable potential
                         alpha = 0.5 * (d1 + d2)
                         L1 = sps.eval_genlaguerre(p, d1, z)
                         L2 = sps.eval_genlaguerre(q, d2, z)
 
                         qvals = np.sqrt(2.0 * z) / float(ell)
-                        Veff = V_of_q(qvals) / (2.0 * np.pi * float(ell) ** 2)
+                        Veff = pot_fn(qvals) / (2.0 * np.pi * float(ell) ** 2)
 
                         if absN not in J_cache:
                             arg = np.sqrt(2.0 * z)[None, :] * Gscaled[:, None]

tests/test_exchange_legendre.py

@@ -30,3 +30,24 @@ def test_legendre_large_n_stability():
     # This should not raise an error
     X = get_exchange_kernels_GaussLegendre(Gs_dimless, thetas, nmax, nquad=500)
     assert np.isfinite(X).all()
+
+
+def test_legendre_callable_potential_matches_coulomb():
+    """Callable potential should reproduce Coulomb when given V(q)=2πκ/q."""
+    nmax = 3
+    Gs_dimless = np.array([0.3, 1.1])
+    thetas = np.array([0.0, 0.7])
+    kappa = 1.2
+    nquad = 300
+
+    def V_coulomb(q):
+        return kappa * 2.0 * np.pi / q
+
+    X_coulomb = get_exchange_kernels_GaussLegendre(
+        Gs_dimless, thetas, nmax, potential="coulomb", kappa=kappa, nquad=nquad
+    )
+    X_callable = get_exchange_kernels_GaussLegendre(
+        Gs_dimless, thetas, nmax, potential=V_coulomb, nquad=nquad
+    )
+
+    assert np.allclose(X_callable, X_coulomb, rtol=1e-4, atol=1e-4)

tests/test_validation.py

@@ -93,3 +93,23 @@ def test_gausslag_convergence():
 
     assert np.allclose(X_low, X_high, rtol=1e-6, atol=1e-4), \
         "Gauss-Laguerre quadrature not converged between nquad=50 and nquad=200"
+
+
+def test_hankel_callable_potential_matches_coulomb():
+    """Hankel backend should respect a user-supplied callable identical to Coulomb."""
+    nmax = 2
+    Gs_dimless = np.array([0.8, 1.6])
+    thetas = np.array([0.0, 0.4])
+    kappa = 0.9
+
+    def V_coulomb(q):
+        return kappa * 2.0 * np.pi / q
+
+    X_coulomb = get_exchange_kernels(
+        Gs_dimless, thetas, nmax, method="hankel", potential="coulomb", kappa=kappa
+    )
+    X_callable = get_exchange_kernels(
+        Gs_dimless, thetas, nmax, method="hankel", potential=V_coulomb
+    )
+
+    assert np.allclose(X_callable, X_coulomb, rtol=1e-4, atol=1e-4)