From 42a601e35509f4ed6fb657512e6095dbdb7962f2 Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Tue, 11 Aug 2026 20:23:00 +0900
Subject: [PATCH] Support custom potentials in Rayleigh-Ritz

---
 Project.toml              |  2 ++
 docs/src/Rayleigh-Ritz.md |  1 +
 src/Rayleigh-Ritz.jl      | 24 ++++++++++++++++++++++++
 test/Rayleigh-Ritz.jl     | 22 ++++++++++++++++++++++
 4 files changed, 49 insertions(+)

diff --git a/Project.toml b/Project.toml
index 74c6852..d8b14f6 100644
--- a/Project.toml
+++ b/Project.toml
@@ -10,6 +10,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
@@ -20,6 +21,7 @@ ArnoldiMethod = "0.4.0"
 FiniteDifferenceMatrices = "0.1.0"
 ForwardDiff = "0.10, 1"
 Optim = "1.9.4"
+QuadGK = "2.11"
 SpecialFunctions = "2.3.1"
 Subscripts = "0.1.3"
 julia = "1.7"
diff --git a/docs/src/Rayleigh-Ritz.md b/docs/src/Rayleigh-Ritz.md
index 8bb2536..1a7f9ed 100644
--- a/docs/src/Rayleigh-Ritz.md
+++ b/docs/src/Rayleigh-Ritz.md
@@ -228,4 +228,5 @@ element(o::Linear, SGB1::SimpleGaussianBasis, SGB2::SimpleGaussianBasis)
 element(o::Coulomb, SGB1::SimpleGaussianBasis, SGB2::SimpleGaussianBasis)
 element(o::PowerLaw, SGB1::SimpleGaussianBasis, SGB2::SimpleGaussianBasis)
 element(o::Gaussian, SGB1::SimpleGaussianBasis, SGB2::SimpleGaussianBasis)
+element(o::Custom, SGB1::SimpleGaussianBasis, SGB2::SimpleGaussianBasis)
 ```
diff --git a/src/Rayleigh-Ritz.jl b/src/Rayleigh-Ritz.jl
index e735cfd..15816e5 100644
--- a/src/Rayleigh-Ritz.jl
+++ b/src/Rayleigh-Ritz.jl
@@ -3,6 +3,7 @@ export solve, optimize
 import LinearAlgebra
 import Optim
 import Printf
+import QuadGK
 import SpecialFunctions
 import Subscripts
 
@@ -412,6 +413,17 @@ function element(o::Gaussian, SGB1::SimpleGaussianBasis, SGB2::SimpleGaussianBas
   return o.coefficient * (π/(o.exponent+SGB1.a+SGB2.a))^(3/2)
 end
 
+function element(o::Custom, SGB1::SimpleGaussianBasis, SGB2::SimpleGaussianBasis)
+  value, _ = QuadGK.quadgk(
+    r -> 4π * r^2 * φ(SGB1, r) * V(o, r) * φ(SGB2, r),
+    0.0,
+    Inf,
+    rtol=1e-10,
+    atol=1e-12,
+  )
+  return value
+end
+
 function element(o::Hamiltonian, B1::Basis, B2::Basis)
   return sum(element(term, B1, B2) for term in o.terms)
 end
@@ -828,6 +840,18 @@ Integral Formula:
 ```
 """ element(o::Gaussian, SGB1::SimpleGaussianBasis, SGB2::SimpleGaussianBasis)
 
+@doc raw"""
+`element(o::Custom, SGB1::SimpleGaussianBasis, SGB2::SimpleGaussianBasis)`
+
+The matrix element for a custom central potential is evaluated numerically with
+adaptive Gauss--Kronrod quadrature:
+```math
+\langle \phi_i | V | \phi_j \rangle
+= 4\pi \int_0^\infty
+r^2 \phi_i(r) V(r) \phi_j(r)\,\mathrm{d}r.
+```
+""" element(o::Custom, SGB1::SimpleGaussianBasis, SGB2::SimpleGaussianBasis)
+
 @doc raw"""
 `element(o::Hamiltonian, SGB1::SimpleGaussianBasis, SGB2::SimpleGaussianBasis)`
 
diff --git a/test/Rayleigh-Ritz.jl b/test/Rayleigh-Ritz.jl
index 18eb4d1..509ed3c 100644
--- a/test/Rayleigh-Ritz.jl
+++ b/test/Rayleigh-Ritz.jl
@@ -121,5 +121,27 @@
       end
     end
 
+    custom = Custom(f=r -> exp(-r^2))
+    gaussian = Gaussian(coefficient=1, exponent=1)
+    for b₁ in BS.basis, b₂ in BS.basis
+      @test TwoBody.element(custom, b₁, b₂) ≈ TwoBody.element(gaussian, b₁, b₂) rtol=1e-9
+    end
+
+  end
+
+  @testset "Custom potential" begin
+    custom_H = Hamiltonian(
+      Kinetic(hbar=1, m=1),
+      Custom(f=r -> -exp(-r^2)),
+    )
+    gaussian_H = Hamiltonian(
+      Kinetic(hbar=1, m=1),
+      Gaussian(coefficient=-1, exponent=1),
+    )
+
+    custom_result = solve(custom_H, BS, info=1)
+    gaussian_result = solve(gaussian_H, BS, info=1)
+    @test custom_result.H ≈ gaussian_result.H rtol=1e-9
+    @test custom_result.E ≈ gaussian_result.E rtol=1e-9
   end
 end