yapanuwan
/
BergEC-jl
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Numerical V matrix elements

This commit is contained in:
Nuwan Yapa 2024-04-26 17:02:32 -04:00
parent dc28873bce
commit 98e8cdb9d3
1 changed files with 14 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
p_space.jl
View File

@ -1,9 +1,5 @@
using LinearAlgebra using LinearAlgebra
using FastGaussQuadrature using SpecialFunctions, FastGaussQuadrature, QuadGK
# Gaussian potentials in momentum space
g0(R, p, q) = (exp(-(1/4)*(p + q)^2*R^2)*(-1 + exp(p*q*R^2))*R)/(2*sqrt(π))
g1(R, p, q) = (exp(-(1/4)*(p + q)^2*R^2)*(2 + p*q*R^2 + exp(p*q*R^2)*(-2 + p*q*R^2)))/(2*p*sqrt(π)*q*R)
function gausslegendre_shifted(a, b, n) function gausslegendre_shifted(a, b, n)
scale = (b - a) / 2 scale = (b - a) / 2
@ -50,4 +46,16 @@ function quick_pole_E(V_pq, μ=0.5; cs_angle=0.4, cutoff=8.0, meshpoints=256)
p, w = get_mesh([0, cutoff * exp(-1im * cs_angle)], [meshpoints]) p, w = get_mesh([0, cutoff * exp(-1im * cs_angle)], [meshpoints])
evals = eigvals(get_H_matrix(V_pq, p, w, μ)) evals = eigvals(get_H_matrix(V_pq, p, w, μ))
return evals[identify_pole_i(p, evals, μ)] return evals[identify_pole_i(p, evals, μ)]
end end
# Gaussian potentials in momentum space
g0(R, p, q) = (exp(-(1/4)*(p + q)^2*R^2)*(-1 + exp(p*q*R^2))*R)/(2*sqrt(π))
g1(R, p, q) = (exp(-(1/4)*(p + q)^2*R^2)*(2 + p*q*R^2 + exp(p*q*R^2)*(-2 + p*q*R^2)))/(2*p*sqrt(π)*q*R)
# general potential (numerical integration)
jHat(l, z) = z * sphericalbesselj(l, z)
function Vl_mat_elem(V_of_r, l, p, q; atol=0, maxevals=10^7, R_cutoff=Inf)
integrand(r) = jHat(l, p * r) * V_of_r(r) * jHat(l, q * r)
(integral, _) = quadgk(integrand, 0, R_cutoff; atol=atol, maxevals=maxevals)
return (2 / pi) * integral
end