Support custom potentials in Rayleigh-Ritz

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"

docs/src/Rayleigh-Ritz.md

@@ -266,4 +266,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)
 ```

src/Rayleigh-Ritz.jl

@@ -3,6 +3,7 @@ export solve, optimize
 import LinearAlgebra
 import Optim
 import Printf
+import QuadGK
 import SpecialFunctions
 import Subscripts
 
@@ -428,6 +429,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
+
 @inline function element(B1::ContractedBasis, B2::Basis)
   return _contracted_bra_element(
     getfield(B1, :coefficients),
@@ -910,6 +922,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)`
 

test/Rayleigh-Ritz.jl

@@ -218,5 +218,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