Quit[]

B-splines

setup and spline pictures

(*for mathematica's numerical quadrature*)

Needs["NumericalDifferentialEquationAnalysis`"]

(*function definition to call the ith B-spline of order ord defined on the knot points "knots" *)

bsp[i_,x_,ord_,knots_]:=BSplineBasis[{ord-1,knots},i-1,x];

(*dummy function for integration later*)

unit[x_]:=1;

(*function to symmetrize a real matrix*)

sym[x_]:=(x+Transpose[x])-DiagonalMatrix[Diagonal[x]];

(*define a set of b-splines, evaluated at the points needed for Gaussian quadrature*)
(*ord = order of the b-splines. numk = number of knot points. rr = right boundary. rl = left boundary. in = integration order, usually 2 or 3 is fine. knotop = linear, quadratic, or expontential determines the knot poisitons. Potential is the potential function. *)

setupBsplines[ord_,numk_,rl_,rr_,in_,knotop_,potential_]:=Module[
{d,n,intx,gqw,bs,bsd,bb,knots,pot0,pot1},
(*note that since BSplineBasis takes in degree rather than order, we subtract 1...*)
d=ord-1;
(*knots are defined on a radial, quadratic, or exponential grid here*)
If[knotop=="Linear",
knots=Table[i,{i,rl,rr,(rr-rl)/numk//N}]];
If[knotop=="Quadratic",
knots=Table[i^2,{i,Sqrt[rl],Sqrt[rr],(Sqrt[rr]-Sqrt[rl])/numk//N}]];
If[knotop=="Exponential",
knots=Table[Exp[i],{i,Log[rl],Log[rr],(Log[rr]-Log[rl])/numk}]//N];
(*add the correct number of buffers on the boundaries to allow open B.C.s*);
knots=Flatten[Join[{Table[knots[[1]],{i,1,d}],Table[knots[[-1]],{i,1,d}],knots}]]//Sort;
n=ord+numk-1;(*total number of splines in this basis with numk knots and order ord*)
intx=100;(*for interpolating splines later*)
(*weights for gaussian quadrature*);
gqw=Table[GaussianQuadratureWeights[in*ord,knots[[i]],knots[[i+1]]],{i,1,Length[knots]-1}];
gqw =Append[gqw,gqw[[-1]]];
(*pre-evaluate the splines and their derivatives at the points required for gaussian quadrature; this will make integrals fast later*);
bs=ParallelTable[bb[x_]:=bsp[i,x,ord,knots];Table[bb/@gqw[[q,All,1]],{q,1,Length[gqw]}],{i,1,n}];
bsd=ParallelTable[bb[x_]:=D[bsp[i,y,ord,knots],y]/.y->x;Table[bb/@gqw[[q,All,1]],{q,1,Length[gqw]}],{i,1,n}];
pot1=Table[potential/@gqw[[q,All,1]]*gqw[[q,All,2]],{q,1,Length[gqw]}];
pot0=Table[unit/@gqw[[q,All,1]]*gqw[[q,All,2]],{q,1,Length[gqw]}];
Return[{bs,bsd,pot1,pot0,gqw,knots,n}]
]

integrateNew[i_,j_,bs_,bsd_,fun_,order_,d1_]:=Module[
{p},
p=If[i>j||Abs[i-j]>=order,
0,
If[d1==0,
Total[Table[Total[(bs[[i,q]])*(bs[[j,q]])*(fun[[q]])],{q,i,j+order}]],
Total[Table[Total[(bsd[[i,q]])*(bsd[[j,q]])*(fun[[q]])],{q,i,j+order}]]
]];
p
]

(*just a test example: 10 closed B-splines + 2 open splines, all of order 3, on a linear mesh*)

{bsa,bsda,p1,p0,gqw,knots,ntot}=setupBsplines[3,10,10,1*^-10,2,"Linear",unit];

(*plot the splines. note that and the first degree spline and the nth spline are one at the boundary. All others are finite over some range but zero at the boundary. This is because we have added the extra buffer points at the knot edges*)
Plot[Table[bsp[i,x,3,knots],{i,1,ntot}],{x,knots[[1]],knots[[-1]]},PlotStyle->{Red},PlotRange->All]

(*notice that if we want to impose boundary conditions such that the wave function vanishes on the boundary (a typical scenario), we need to exclude the degree one and degree numK states from our basis*)Plot[Table[bsp[i,x,3,knots],{i,2,ntot-1}],{x,knots[[1]],knots[[-1]]},PlotStyle->{Red},PlotRange->All]

(*all of our splines togehter should sum to 1 over the whole range*)
Plot[Evaluate[Sum[bsp[i,x,3,knots],{i,1,ntot}]],{x,knots[[1]],knots[[-1]]},PlotRange->{0.99999,1.00001}]

(*some quadratic and exponential grid examples*)
{bsa,bsda,p1,p0,gqw,knots,ntot}=setupBsplines[3,10,10,1*^-10,2,"Quadratic",unit];
Plot[Table[bsp[i,x,3,knots],{i,1,ntot}],{x,knots[[1]],knots[[-1]]},PlotStyle->{Red},PlotRange->All]

(*some quadratic and exponential grid examples*)
{bsa,bsda,p1,p0,gqw,knots,ntot}=setupBsplines[3,10,10,1*^-10,2,"Exponential",unit];
Plot[Table[bsp[i,x,3,knots],{i,1,ntot}],{x,knots[[1]],knots[[-1]]},PlotStyle->{Red},PlotRange->All]

Spline demonstrations

Harmonic potential.

pot[x_]:=1/2*x^2

μred=1;

ord=6;

{bsa,bsda,p1,p0,gqw,knots,ntot}=setupBsplines[ord,100,-10,10,3,"Linear",pot];

Hdense=sym[ParallelTable[1/(2*μred)*integrateNew[i,j,bsa,bsda,p0,ord,1]+integrateNew[i,j,bsa,bsda,p1,ord,0],{i,2,ntot-1},{j,2,ntot-1}]];
Hsparse=SparseArray[Flatten[Join[{Table[Band[{i,1}]->Diagonal[Hdense,i-1],{i,1,ord+1,1}],Table[Band[{1,i}]->Diagonal[Hdense,i-1],{i,2,ord+1,1}]}],1]];

OVdense = sym[Table[integrateNew[i,j,bsa,bsda,p0,ord,0],{i,2,ntot-1},{j,2,ntot-1}]];
OVsparse=SparseArray[Flatten[Join[{Table[Band[{i,1}]->Diagonal[OVdense,i-1],{i,1,ord+1,1}],Table[Band[{1,i}]->Diagonal[OVdense,i-1],{i,2,ord+1,1}]}],1]];

{v,v2}=Eigensystem[{Hsparse,OVsparse}];

v//Sort

ho[x_,n_]:=1/Sqrt[2^n*n!]/Pi^0.25*Exp[-x^2/2]*HermiteH[n,x]

sol=v2[[-1]];
norm=sol.(OVsparse.sol)//Sqrt;
Plot[{sol.Table[bsp[i,x,ord,knots],{i,2,ntot-1}]/norm,ho[x,0]},{x,-10,10},PlotRange->All,PlotStyle->{{Red,Thick},{Blue,Dashed}}]

Coulomb potential, l = 0

pot[x_]:=-1/x+0*(0+1)/2/x^2

Plot[pot[x],{x,1*^-8,30}]

μred=1;

ord=6;

{bsa,bsda,p1,p0,gqw,knots,ntot}=setupBsplines[ord,100,1*^-8,100,3,"Quadratic",pot];

Hdense=sym[ParallelTable[1/(2*μred)*integrateNew[i,j,bsa,bsda,p0,ord,1]+integrateNew[i,j,bsa,bsda,p1,ord,0],{i,2,ntot-1},{j,2,ntot-1}]];
Hsparse=SparseArray[Flatten[Join[{Table[Band[{i,1}]->Diagonal[Hdense,i-1],{i,1,ord+1,1}],Table[Band[{1,i}]->Diagonal[Hdense,i-1],{i,2,ord+1,1}]}],1]];

OVdense = sym[Table[integrateNew[i,j,bsa,bsda,p0,ord,0],{i,2,ntot-1},{j,2,ntot-1}]];
OVsparse=SparseArray[Flatten[Join[{Table[Band[{i,1}]->Diagonal[OVdense,i-1],{i,1,ord+1,1}],Table[Band[{1,i}]->Diagonal[OVdense,i-1],{i,2,ord+1,1}]}],1]];

shift=5;

{v,v2}=Eigensystem[{Hsparse+shift*OVsparse,OVsparse}];

v-shift//Sort

Table[-1/2/n^2,{n,1,10,1}]//N

sol=v2[[-2]];
norm=sol.(OVsparse.sol)//Sqrt;
Plot[{sol.Table[bsp[i,x,ord,knots],{i,2,ntot-1}]/norm},{x,0,50},PlotRange->All,PlotStyle->{{Red,Thick},{Blue,Dashed}}]

Coulomb potential, l = 3

pot[x_]:=-1/x+3*(3+1)/2/x^2

Plot[pot[x],{x,1*^-8,30}]

μred=1;

ord=6;

{bsa,bsda,p1,p0,gqw,knots,ntot}=setupBsplines[ord,100,1*^-8,100,3,"Quadratic",pot];

Hdense=sym[ParallelTable[1/(2*μred)*integrateNew[i,j,bsa,bsda,p0,ord,1]+integrateNew[i,j,bsa,bsda,p1,ord,0],{i,2,ntot-1},{j,2,ntot-1}]];
Hsparse=SparseArray[Flatten[Join[{Table[Band[{i,1}]->Diagonal[Hdense,i-1],{i,1,ord+1,1}],Table[Band[{1,i}]->Diagonal[Hdense,i-1],{i,2,ord+1,1}]}],1]];

OVdense = sym[Table[integrateNew[i,j,bsa,bsda,p0,ord,0],{i,2,ntot-1},{j,2,ntot-1}]];
OVsparse=SparseArray[Flatten[Join[{Table[Band[{i,1}]->Diagonal[OVdense,i-1],{i,1,ord+1,1}],Table[Band[{1,i}]->Diagonal[OVdense,i-1],{i,2,ord+1,1}]}],1]];

shift=5;

{v,v2}=Eigensystem[{Hsparse+shift*OVsparse,OVsparse}];

v-shift//Sort

Table[-1/2/n^2,{n,4,10,1}]//N

sol=v2[[-3]];
norm=sol.(OVsparse.sol)//Sqrt;
Plot[{sol.Table[bsp[i,x,ord,knots],{i,2,ntot-1}]/norm},{x,0,100},PlotRange->All,PlotStyle->{{Red,Thick},{Blue,Dashed}}]

Resonances: setup

(*we want to solve for the scattering phase shifts, resonance positions, etc of the following potential [this is a common reference one apparently]*)

μred=1;

pot[x_]:=7.5*x^2*Exp[-x]

Plot[pot[x],{x,0,20}]

Resonance solver #1: stabilization.
reference: https://journals.aps.org/pra/pdf/10.1103/PhysRevA.49.2470

(*our goal here is just to repeatedly diagonalize the Hamiltonian, obtaining bound state eigenvalues as a function of the outer box position. 
Those which 'stabilize' as a function of this distance are resonance positions. *)

ord=6;

stabloop=Table[{rmax1,
{bsa,bsda,p1,p0,gqw,knots,ntot}=setupBsplines[ord,50,1*^-8,rmax1,3,"Quadratic",pot];
Hdense=sym[Table[1/(2*μred)*integrateNew[i,j,bsa,bsda,p0,ord,1]+integrateNew[i,j,bsa,bsda,p1,ord,0],{i,2,ntot-1},{j,2,ntot-1}]];
Hsparse=SparseArray[Flatten[Join[{Table[Band[{i,1}]->Diagonal[Hdense,i-1],{i,1,ord+1,1}],Table[Band[{1,i}]->Diagonal[Hdense,i-1],{i,2,ord+1,1}]}],1]];
OVdense = sym[Table[integrateNew[i,j,bsa,bsda,p0,ord,0],{i,2,ntot-1},{j,2,ntot-1}]];
OVsparse=SparseArray[Flatten[Join[{Table[Band[{i,1}]->Diagonal[OVdense,i-1],{i,1,ord+1,1}],Table[Band[{1,i}]->Diagonal[OVdense,i-1],{i,2,ord+1,1}]}],1]];
v=Eigenvalues[{Hsparse,OVsparse}]},{rmax1,25,30,0.1}];//AbsoluteTiming

(*This is a big list of (r0,ej), where r0 is the outer position and ej are the eigenvalues. We save time on mathematica's histogram by looking only at energies <20*)

plotvals=Flatten[Table[Table[{stabloop[[i,1]],stabloop[[i,2,j]]},{i,1,Length[stabloop]}],{j,1,Length[stabloop[[1,2]]]}],1];
pv={};
Do[If[plotvals[[i,2]]<20,pv=Append[pv,plotvals[[i,2]]]],{i,1,Length[plotvals]}]

Histogram[pv,1000]

Show[Histogram[pv,1000],PlotRange->{{3,4},{0,50}}]

(*from https://journals.aps.org/pra/pdf/10.1103/PhysRevA.49.2470 (Eq. 7) - 2/Γ is given by Pi*max(hist)/I, where I is #bins/energy / # r values.*)

Show[Histogram[pv,1000],PlotRange->{{3.38,3.5},{0,25}}]

maxhist=24;
binsperenergy=1/0.02;
rvalues=Length[stabloop];

Gamover2=1/(Pi*maxhist/(binsperenergy/rvalues))

(*this is the usual sort of stabilization plot*)

ListPlot[plotvals,PlotRange->{3,4}]

Resonance solver #2: NDsolve.
reference (for time delay): https://journals.aps.org/pra/pdf/10.1103/PhysRevA.45.1587 and https://iopscience.iop.org/article/10.1088/1361-6455/aceb28/pdf

(*With this method, we solve directly the Schroedinger equation using NDSolve, and extract the asymptotic phase by comparing it to the asymptotic solutions*)

solu[l_,en_,rmin_,rmax_]:=ψ/.NDSolve[{-1/2/μred*D[ψ[r],{r,2}]+(pot[r]+l*(l+1)/2/μred/r^2-en)*ψ[r]==0,ψ[rmin]==0.,ψ'[rmin]==Sqrt[2*en*μred]*rmin},ψ,{r,rmin,rmax}][[1]]

tanδ[l_,en_,rmin_,rmax_]:=Module[{
ψs,ψsp,ψse,ψspe,f,fp,g,gp},
ψs=solu[l,en,rmin,rmax];
ψsp=ψs';
ψse=ψs[rmax];
ψspe=ψsp[rmax];

f=SetPrecision[Sqrt[2*μred*en]*rmax*SphericalBesselJ[l,Sqrt[2*μred*en]*rmax],20];
g=SetPrecision[Sqrt[2*μred*en]*rmax*SphericalBesselY[l,Sqrt[2*μred*en]*rmax],20];
fp=SetPrecision[D[Sqrt[2*μred*en]*r*SphericalBesselJ[l,Sqrt[2*μred*en]*r],r]/.r->rmax,20];
gp=SetPrecision[D[Sqrt[2*μred*en]*r*SphericalBesselY[l,Sqrt[2*μred*en]*r],r]/.r->rmax,20];

SetPrecision[(fp*ψse-f*ψspe)/(gp*ψse-g*ψspe),20]
]

phase=Table[{en,ArcTan[tanδ[0,en,0.0001,20]]},{en,0.00001,5,0.001}];

ListLinePlot[phase,PlotRange->All]

Do[If[phase[[i+1,2]]-phase[[i,2]]>=Pi*0.9,phase[[i+1;;Length[phase],2]]-=Pi],{i,1,Length[phase]-1}]

ListLinePlot[phase,PlotRange->All]

(*the time delay is given by this, twice the energy derivative of δ(E)*)

timedelay=Table[{phase[[i,1]],2*(phase[[i+1,2]]-phase[[i-1,2]])/(phase[[i+1,1]]-phase[[i-1,1]])},{i,2,Length[phase]-1}];

ListPlot[timedelay,PlotRange->{{3,4},{-300,500}}]

(*in turn, this equals Eq. 10 of https://journals.aps.org/pra/pdf/10.1103/PhysRevA.45.1587. From this we see that Γ/2 is given by the max height divided by 2, inverse. *)

2/159.

Resonance solver #3: Siegert states.
reference: https://journals.aps.org/pra/pdf/10.1103/PhysRevA.58.2077

(*Now we are gonna set up the Siegert problem using a outer boundary of rmax1 *)

ord=12;
rmax1=25;
{bsa,bsda,p1,p0,gqw,knots,ntot}=setupBsplines[ord,80,1*^-8,rmax1,4,"Linear",pot];
Hdense=sym[Table[1/(2*μred)*integrateNew[i,j,bsa,bsda,p0,ord,1]+integrateNew[i,j,bsa,bsda,p1,ord,0],{i,2,ntot},{j,2,ntot}]];
Lmat=Table[Table[bsp[i,rmax1,ord,knots]*bsp[j,rmax1,ord,knots],{i,2,ntot}],{j,2,ntot}];
OVdense = sym[Table[integrateNew[i,j,bsa,bsda,p0,ord,0],{i,2,ntot},{j,2,ntot}]];

A=2*Hdense;
B=-Lmat;

lhs=ArrayFlatten[{{0*IdentityMatrix[ntot-1],IdentityMatrix[ntot-1]},{-A,-B}}];
rhs=ArrayFlatten[{{IdentityMatrix[ntot-1],0*IdentityMatrix[ntot-1]},{0*IdentityMatrix[ntot-1],OVdense}}];

kval=Eigenvalues[{lhs,rhs}]/I;

{en,gamby2}={Re[kval^2/2],-Im[kval^2/2]};

so=Table[{en[[i]],gamby2[[i]]},{i,1,Length[en]}]

ListPlot[so]

ListPlot[Table[{Re[kval[[i]]],Im[kval[[i]]]},{i,1,Length[kval]}],PlotRange->{{-10,10},{-1,0}}]

S[k_]:=Exp[-2*I*k*rmax1]*Product[(kval[[i]]+k)/(kval[[i]]-k),{i,1,Length[kval]}]

sinsqδ[k_]:=Sin[Log[S[k]]/2/I]^2

ListLinePlot[Table[{kv,sinsqδ[10^kv]},{kv,-1,1,0.001}],PlotRange->{0,1},PlotStyle->Black]

Show[{,,,}]

δ[k_]:=Log[S[k]]/2/I

phases=Table[{kv^2/2,Re[δ[kv]]},{kv,Sqrt[0.01],Sqrt[2*5],0.001}];

ListLinePlot[phases,PlotStyle->Black]

Do[If[phases[[i+1,2]]-phases[[i,2]]>=Pi*0.9,phases[[i+1;;Length[phases],2]]-=Pi],{i,1,Length[phases]-1}]

Show[{ListLinePlot[phases,PlotStyle->Black],},PlotRange->All]

Resonance solver #4: Complex absorbing potential
reference: https://iopscience.iop.org/article/10.1088/0953-4075/26/23/021/pdf and https://journals.aps.org/pra/pdf/10.1103/PhysRevA.74.034701

First, for an optimized parameter set

rc=8;
η=0.76*^-2;

pot2[x_]:=HeavisideTheta[(x-rc)]*(x-rc)^2

pot2[10.2]

ord=12;
rmax1=25;
{bsa,bsda,p1,p0,gqw,knots,ntot}=setupBsplines[ord,80,1*^-8,rmax1,4,"Linear",pot];
Hdense=sym[Table[1/(2*μred)*integrateNew[i,j,bsa,bsda,p0,ord,1]+integrateNew[i,j,bsa,bsda,p1,ord,0],{i,2,ntot-1},{j,2,ntot-1}]];
OVdense = sym[Table[integrateNew[i,j,bsa,bsda,p0,ord,0],{i,2,ntot-1},{j,2,ntot-1}]];
{bsa,bsda,p1,p0,gqw,knots,ntot}=setupBsplines[ord,80,1*^-8,rmax1,4,"Linear",pot2];
W=-I*η*Table[integrateNew[i,j,bsa,bsda,p1,ord,0],{i,2,ntot-1},{j,2,ntot-1}];

vv=Eigenvalues[{Hdense+W,OVdense}]

ListPlot[Table[{Re[vv[[i]]],Im[vv[[i]]]},{i,1,Length[vv]}],PlotRange->{{0,20},{-3,0}}]

Now for something less optimized

rc=8;
η=0.76*^-1;

pot2[x_]:=HeavisideTheta[(x-rc)]*(x-rc)^2

pot2[10.2]

ord=12;
rmax1=25;
{bsa,bsda,p1,p0,gqw,knots,ntot}=setupBsplines[ord,80,1*^-8,rmax1,4,"Linear",pot];
Hdense=sym[Table[1/(2*μred)*integrateNew[i,j,bsa,bsda,p0,ord,1]+integrateNew[i,j,bsa,bsda,p1,ord,0],{i,2,ntot-1},{j,2,ntot-1}]];
OVdense = sym[Table[integrateNew[i,j,bsa,bsda,p0,ord,0],{i,2,ntot-1},{j,2,ntot-1}]];
{bsa,bsda,p1,p0,gqw,knots,ntot}=setupBsplines[ord,80,1*^-8,rmax1,4,"Linear",pot2];
W=-I*η*Table[integrateNew[i,j,bsa,bsda,p1,ord,0],{i,2,ntot-1},{j,2,ntot-1}];

vv=Eigenvalues[{Hdense+W,OVdense}]

ListPlot[Table[{Re[vv[[i]]],Im[vv[[i]]]},{i,1,Length[vv]}],PlotRange->{{0,20},{-3,0}}]

Now for something far less optimized

rc=8;
η=10.76;

pot2[x_]:=HeavisideTheta[(x-rc)]*(x-rc)^2

pot2[10.2]

ord=12;
rmax1=25;
{bsa,bsda,p1,p0,gqw,knots,ntot}=setupBsplines[ord,80,1*^-8,rmax1,4,"Linear",pot];
Hdense=sym[Table[1/(2*μred)*integrateNew[i,j,bsa,bsda,p0,ord,1]+integrateNew[i,j,bsa,bsda,p1,ord,0],{i,2,ntot-1},{j,2,ntot-1}]];
OVdense = sym[Table[integrateNew[i,j,bsa,bsda,p0,ord,0],{i,2,ntot-1},{j,2,ntot-1}]];
{bsa,bsda,p1,p0,gqw,knots,ntot}=setupBsplines[ord,80,1*^-8,rmax1,4,"Linear",pot2];
W=-I*η*Table[integrateNew[i,j,bsa,bsda,p1,ord,0],{i,2,ntot-1},{j,2,ntot-1}];

vv=Eigenvalues[{Hdense+W,OVdense}]

ListPlot[Table[{Re[vv[[i]]],Im[vv[[i]]]},{i,1,Length[vv]}],PlotRange->{{0,20},{-3,0}}]

Resonance solver #5: R-matrix method. 
reference: https://journals.aps.org/pra/pdf/10.1103/PhysRevA.58.2077

first: do it using the streamlined version via linear solve

(*Now we are gonna set up the R-matrix problem using a outer boundary of rmax1 *)

ord=12;
rmax1=25;
{bsa,bsda,p1,p0,gqw,knots,ntot}=setupBsplines[ord,80,1*^-8,rmax1,4,"Linear",pot];
H=sym[Table[1/(2*μred)*integrateNew[i,j,bsa,bsda,p0,ord,1]+integrateNew[i,j,bsa,bsda,p1,ord,0],{i,2,ntot},{j,2,ntot}]];
OV = sym[Table[integrateNew[i,j,bsa,bsda,p0,ord,0],{i,2,ntot},{j,2,ntot}]];

n=ntot

ΛOO =bsp[ntot,rmax1,ord,knots]*bsp[ntot,rmax1,ord,knots];
ΛCC=Table[bsp[i,rmax1,ord,knots]*bsp[j,rmax1,ord,knots],{i,2,n-1},{j,2,n-1}];

ΓOO[energy_]:= (2*μred)*(energy*OV[[n-1,n-1]]-H[[n-1,n-1]]);
ΓOC[energy_]:=(2*μred)(energy*OV[[n-1,1;;n-2]]-H[[n-1,1;;n-2]]);
ΓCO[energy_]:= (2*μred)*(energy*OV[[1;;n-2,n-1]]-H[[1;;n-2,n-1]]);
ΓCC[energy_]:= (2*μred)*(energy*OV[[1;;n-2,1;;n-2]]-H[[1;;n-2,1;;n-2]]);
tanδ[b_,kv1_,l_] := (b*SphericalBesselJ[l,kv1*rmax1 ]*rmax1+D[SphericalBesselJ[l,kv1*y]*y,y]/.y->rmax1)/(b*SphericalBesselY[l,kv1*rmax1 ]rmax1+D[SphericalBesselY[l,kv1*y]*y,y]/.y->rmax1);

ψrOutside [r_,kv1_,shiftval_] :=(r*Cos[shiftval]*SphericalBesselJ[0,kv1*r]-Sin[shiftval]*r*SphericalBesselY[0,kv1*r])/(rmax1*Cos[shiftval]*SphericalBesselJ[0,kv1*rmax1]-Sin[shiftval]*rmax1*SphericalBesselY[0,kv1*rmax1]);
ψrInside[r_,kv1_,shift_,b1_]:=(Total[Table[b1[[i-1]]*bsp[i,r,ord,knots],{i,2,n-1}]]+bsp[n,r]);
kv[e_]:=Sqrt[2*μred*e]
X[e_]:=LinearSolve[ΓCC[e], -ΓCO[e]];
b[e_]:=(ΓOC[e].X[e]+ΓOO[e])/ΛOO;
shift[e_]:=tanδ[b[e],kv[e],0];

Energy=Table[e,{e,0,5,0.01}];

kvT=Sqrt[Energy*2*μred];

XT=ParallelTable[LinearSolve[ΓCC[Energy[[e]]], -ΓCO[Energy[[e]]]],{e,1,Length[Energy]}];//AbsoluteTiming
bT = ParallelTable[(ΓOC[Energy[[e]]].XT[[e]]+ΓOO[Energy[[e]]])/ΛOO,{e,1,Length[Energy]}];//AbsoluteTiming
dT = ParallelTable[ArcTan[tanδ[bT[[e]],kvT[[e]],0]], {e,1,Length[Energy]}];//AbsoluteTiming
Do[Do[If[dT[[i+1]]-dT[[i]]>0.95Pi,dT[[i+1;;Length[dT]]]-=Pi],{i,1,Length[dT]-1}];,{q,1,1}];
ListPlot[Table[{Energy[[i]],dT[[i]]-Pi},{i,1,Length[Energy]}]]

Plot[b[e],{e,0,5}]

Show[{,}]

second: do it using the streamlined version via diagonalization of Hcc

(*Now we are gonna set up the R-matrix problem using a outer boundary of rmax1 *)

ord=12;
rmax1=25;
{bsa,bsda,p1,p0,gqw,knots,ntot}=setupBsplines[ord,80,1*^-8,rmax1,4,"Linear",pot];
H=sym[Table[1/(2*μred)*integrateNew[i,j,bsa,bsda,p0,ord,1]+integrateNew[i,j,bsa,bsda,p1,ord,0],{i,2,ntot},{j,2,ntot}]];
OV = sym[Table[integrateNew[i,j,bsa,bsda,p0,ord,0],{i,2,ntot},{j,2,ntot}]];

n=ntot;
rmax=rmax1;

ΛOO =bsp[ntot,rmax1,ord,knots]*bsp[ntot,rmax1,ord,knots];
ΛCC=Table[bsp[i,rmax1,ord,knots]*bsp[j,rmax1,ord,knots],{i,2,n-1},{j,2,n-1}];

(*the reason we had to do a linear solve in the previous version is because we use a non-orthogonal B-spline basis. 
If an orthogonal basis is used, we can get all of the needed information from diagonalizing Hcc once. See Eqs. 3.44-3.50 of the orange review.*)

(*we just need to orthogonalize our basis for the closed orbitals then! that is accomplished below, and the we can diagonalize the closed Hamiltonian in this representation...*)

{ov,ovv}=Eigensystem[OV[[1;;-2,1;;-2]]];
ovv=Transpose[ovv];
realovv=IdentityMatrix[ntot-1];
Do[realovv[[i,j]]=ovv[[i,j]],{i,1,ntot-2},{j,1,ntot-2}];
realov=Append[ov,1];

OVdiagrep=DiagonalMatrix[1/Sqrt[realov]].Transpose[realovv].OV.realovv.DiagonalMatrix[1/Sqrt[realov]];

OVdiagrep//Chop//MatrixForm

matdum=Table[Table[If[i==j==ntot,1,0],{i,2,ntot}],{j,2,ntot}];

test2=DiagonalMatrix[1/Sqrt[realov]].Transpose[realovv].matdum.realovv.DiagonalMatrix[1/Sqrt[realov]];

test2//MatrixForm

Hdiagrep=DiagonalMatrix[1/Sqrt[realov]].Transpose[realovv].H.realovv.DiagonalMatrix[1/Sqrt[realov]];

(*diagonalize the closed part*)

{vH,evH}=Eigensystem[Hdiagrep[[1;;n-2,1;;n-2]]];
evH=Transpose[evH];

OnewOC=OVdiagrep[[n-1,1;;n-2]].evH;
HnewOC=Hdiagrep[[n-1,1;;n-2]].evH;

evH//Dimensions

(*Eq. 3.50 of Orange review*)

Ω[ee_]:=2*μred*(ee*OVdiagrep[[n-1,n-1]]-Hdiagrep[[n-1,n-1]])-(2*μred)*Total[Table[(ee*OnewOC[[λ]]-HnewOC[[λ]])*(ee*OnewOC[[λ]]-HnewOC[[λ]])/(ee-vH[[λ]]),{λ,1,Length[vH]}]];

tanδ[b_,kv1_,l_] := (b*SphericalBesselJ[l,kv1*rmax1 ]*rmax1+D[SphericalBesselJ[l,kv1*y]*y,y]/.y->rmax1)/(b*SphericalBesselY[l,kv1*rmax1 ]rmax1+D[SphericalBesselY[l,kv1*y]*y,y]/.y->rmax1);
ψrOutside [r_,kv1_,shiftval_] :=(r*Cos[shiftval]*SphericalBesselJ[0,kv1*r]-Sin[shiftval]*r*SphericalBesselY[0,kv1*r])/(rmax1*Cos[shiftval]*SphericalBesselJ[0,kv1*rmax1]-Sin[shiftval]*rmax1*SphericalBesselY[0,kv1*rmax1]);
ψrInside[r_,kv1_,shift_,b1_]:=(Total[Table[b1[[i-1]]*bsp[i,r,ord,knots],{i,2,n-1}]]+bsp[n,r]);
kv[e_]:=Sqrt[2*μred*e]
b1[e_]:=Ω[e]/ΛOO;

Energy1=Table[e,{e,0,5,0.01}];

kvT1=Sqrt[Energy1*2*μred];

bT1 = Table[b1[Energy1[[e]]],{e,1,Length[Energy1]}];//AbsoluteTiming
dT1 = ParallelTable[ArcTan[tanδ[bT1[[e]],kvT1[[e]],0]], {e,1,Length[Energy1]}];//AbsoluteTiming
Do[Do[If[dT[[i+1]]-dT[[i]]>0.95Pi,dT[[i+1;;Length[dT]]]-=Pi],{i,1,Length[dT]-1}];,{q,1,1}];
ListPlot[Table[{Energy[[i]],dT[[i]]-Pi},{i,1,Length[Energy]}]]

third: a better variation of version 2, which will still be able to take advantage of diagonalizing sparse matrices.

(*Now we are gonna set up the R-matrix problem using a outer boundary of rmax1 *)

ord=12;
rmax1=25;
{bsa,bsda,p1,p0,gqw,knots,ntot}=setupBsplines[ord,80,1*^-8,rmax1,4,"Linear",pot];
H=sym[Table[1/(2*μred)*integrateNew[i,j,bsa,bsda,p0,ord,1]+integrateNew[i,j,bsa,bsda,p1,ord,0],{i,2,ntot},{j,2,ntot}]];
OV = sym[Table[integrateNew[i,j,bsa,bsda,p0,ord,0],{i,2,ntot},{j,2,ntot}]];

n=ntot;
rmax=rmax1;

ΛOO =bsp[ntot,rmax1,ord,knots]*bsp[ntot,rmax1,ord,knots];
ΛCC=Table[bsp[i,rmax1,ord,knots]*bsp[j,rmax1,ord,knots],{i,2,n-1},{j,2,n-1}];

(*the reason we had to do a linear solve in the previous version is because we use a non-orthogonal B-spline basis. 
If an orthogonal basis is used, we can get all of the needed information from diagonalizing Hcc once. See Eqs. 3.44-3.50 of the orange review.*)

(*we just need to orthogonalize our basis for the closed orbitals then! that is accomplished below, and the we can diagonalize the closed Hamiltonian in this representation...*)

{ov,ovv}=Eigensystem[OV[[1;;-2,1;;-2]]];
ovv=Transpose[ovv];
realovv=IdentityMatrix[ntot-1];
Do[realovv[[i,j]]=ovv[[i,j]],{i,1,ntot-2},{j,1,ntot-2}];
realov=Append[ov,1];

OVdiagrep=DiagonalMatrix[1/Sqrt[realov]].Transpose[realovv].OV.realovv.DiagonalMatrix[1/Sqrt[realov]];

Hdiagrep=DiagonalMatrix[1/Sqrt[realov]].Transpose[realovv].H.realovv.DiagonalMatrix[1/Sqrt[realov]];

(*diagonalize the closed part, using the NONORTHOGONAL basis representation here and still benefiting from sparse matrices (I think?)*)

{vH,evH}=Eigensystem[{H[[1;;n-2,1;;n-2]],OV[[1;;n-2,1;;n-2]]}];
evH=Transpose[evH];

(*now here we have to transform both to the representation that makes Hcc diagonal, but also the one tha tmakes OV diagonal...*)

OnewOC=OVdiagrep[[n-1,1;;n-2]].evH.(realovv.DiagonalMatrix[1/Sqrt[realov]])[[1;;n-2,1;;n-2]];
HnewOC=Hdiagrep[[n-1,1;;n-2]].evH.(realovv.DiagonalMatrix[1/Sqrt[realov]])[[1;;n-2,1;;n-2]];

(*Eq. 3.50 of Orange review*)

Ω[ee_]:=2*μred*(ee*OVdiagrep[[n-1,n-1]]-Hdiagrep[[n-1,n-1]])-(2*μred)*Total[Table[(ee*OnewOC[[λ]]-HnewOC[[λ]])*(ee*OnewOC[[λ]]-HnewOC[[λ]])/(ee-vH[[λ]]),{λ,1,Length[vH]}]];

tanδ[b_,kv1_,l_] := (b*SphericalBesselJ[l,kv1*rmax1 ]*rmax1+D[SphericalBesselJ[l,kv1*y]*y,y]/.y->rmax1)/(b*SphericalBesselY[l,kv1*rmax1 ]rmax1+D[SphericalBesselY[l,kv1*y]*y,y]/.y->rmax1);
ψrOutside [r_,kv1_,shiftval_] :=(r*Cos[shiftval]*SphericalBesselJ[0,kv1*r]-Sin[shiftval]*r*SphericalBesselY[0,kv1*r])/(rmax1*Cos[shiftval]*SphericalBesselJ[0,kv1*rmax1]-Sin[shiftval]*rmax1*SphericalBesselY[0,kv1*rmax1]);
ψrInside[r_,kv1_,shift_,b1_]:=(Total[Table[b1[[i-1]]*bsp[i,r,ord,knots],{i,2,n-1}]]+bsp[n,r]);
kv[e_]:=Sqrt[2*μred*e]
b1[e_]:=Ω[e]/ΛOO;

Energy1=Table[e,{e,0,5,0.01}];

kvT1=Sqrt[Energy1*2*μred];

bT1 = Table[b1[Energy1[[e]]],{e,1,Length[Energy1]}];//AbsoluteTiming
dT1 = ParallelTable[ArcTan[tanδ[bT1[[e]],kvT1[[e]],0]], {e,1,Length[Energy1]}];//AbsoluteTiming
Do[Do[If[dT[[i+1]]-dT[[i]]>0.95Pi,dT[[i+1;;Length[dT]]]-=Pi],{i,1,Length[dT]-1}];,{q,1,1}];
ListPlot[Table[{Energy[[i]],dT[[i]]-Pi},{i,1,Length[Energy]}]]