BeginPackage["hylleraas`"]

hmatrix::usage =
   "hmatrix[a,b,c,i1,j1,k1,i2,j2,k2,ndim] returnes three matrix elements\n
 {<i1,j1,k1|T|i2,j2,k2>, <i1,j1,k1|V0|i2,j2,k2>, <i1,j1,k1|V1|i2,j2,k2>, <i1,j1,k1|S|i2,j2,k2>}\n
 of operators T - kinetic energy, V0 = 1/r_1 + 1/r_2, V1 = 1/r_12, S = I\n
 on Hylleraas basis set\n
 |i1,j1,k1> =\n
    [r1^i1 r2^j1 E^(-a r1) E^(-b r2) + r2^i1 r1^j1 E^(-a r2) E^(-b r1)] r12^k1 E^(-c r12)\n
 in the form of expression through Hylleraas integrals.\n
Note that the only dependence on a, b, c is in <i1,j1,k1|T|i2,j2,k2>.\n\n
ndim is dimensionality of the coordinate space\n
  (usually, ndim=3 for the ground or ndim=5 for 2p2 3P-state).\n
Variables a,b,c may be numbers or undefined symbols.\n
Variables i1,j1,k1,i2,j2,k2 should be undefined symbols.\n
ndim should be an odd positive integer larger or equal than 3.\n\n
To perform Hylleraas matrix calculations,\n
it is recommended firstly to calculate the elements in symbolic form by calling hmatrix.\n
Then, Hylleraas integrals are calculated numerically\n
   for particular a, b, c, nm\n
   and finally they are substituted into already known expressions for matrix elements"

hintegrals::usage =
   "hintegrals[a,b,c,ndim,nm] computes numerically the set of Hylleraas integrals\n
 entering expressions for Hylleraas matrix elements\n
 {<i1,j1,k1|T|i2,j2,k2>, <i1,j1,k1|V0|i2,j2,k2>, <i1,j1,k1|V1|i2,j2,k2>, <i1,j1,k1|S|i2,j2,k2>}\n
 of operators T - kinetic energy, V0 = 1/r_1 + 1/r_2, V1 = 1/r_12, S = I\n
 on Hylleraas basis set\n
 |i1,j1,k1> =\n
    [r1^i1 r2^j1 E^(-a r1) E^(-b r2) + r2^i1 r1^j1 E^(-a r2) E^(-b r1)] r12^k1 E^(-c r12)\n
 with i1+j1+k1 <= nm and i2+j2+k2 <= nm.\n\n
ndim is dimensionality of the coordinate space\n
  (usually, ndim=3 for the ground or ndim=5 for 2p2 3P-state).\n
Variables a,b,c may be numbers or undefined symbols.\n
ndim should be an odd positive integer larger or equal than 3.\n
nm should be a non-negative integer.\n\n
Returns nothing."

hmapping::usage =
   "hmapping[nm,choice,ind_,mmax_] maps a number m of the member of the Hylleraas basis set\n
 into the triple index {i,j,k}=ind[m] of the Hylleraas function\n
|i1,j1,k1> =\n
    [r1^i1 r2^j1 E^(-a r1) E^(-b r2) + r2^i1 r1^j1 E^(-a r2) E^(-b r1)] r12^k1 E^(-c r12).\n
If choice!=1, then the standard restriction on indexes is supposed, i+j+k <= nm,\n
 if choice=1, then more strong restriction is supposed, i+j^2+k^2 <= nm.\n
Calculates the size of the basis set mm=mmax[nm] as a function of nm (up to the given nm)\n
 and the mapping function ind[m] (m=1,...,mm).\n
Returnes the basis set size mm=mmax[nm].\n\n
nm should be a non-negative integer.\n
choice should be an integer (1 or something different).\n
ind and mmax should be symbols (used as names of functions)."

Begin["`private`"]

hmatrix[a0_,b0_,c0_,i1_,j1_,k1_,i2_,j2_,k2_,ndim_Integer]:=
(* Body of the package starts *)(
psi1=(r1^i1 r2^j1 r3^k1 E^(-a r1) E^(-b r2)+
        r2^i1 r1^j1 r3^k1 E^(-a r2) E^(-b r1)) E^(-c r3);
psi2=(r1^i2 r2^j2 r3^k2 E^(-a r1) E^(-b r2)+
        r2^i2 r1^j2 r3^k2 E^(-a r2) E^(-b r1)) E^(-c r3);
(* T = <T> *)
Tpsi2=-1/
        2((ndim-1)(1/r1 D[psi2,r1]+1/r2 D[psi2,r2]+2/r3 D[psi2,r3])+
          D[psi2,{r1,2}]+
          D[psi2,{r2,2}]+2 D[psi2,{r3,2}]+(r1^2+r3^2-r2^2)/(r1 r3)*
            D[psi2,r1,r3]+(r2^2+r3^2-r1^2)/(r2 r3) D[psi2,r2,r3]);
ce=Expand[
    psi1 Tpsi2 (r1 r2 r3)^3 (
          r1^2 r2^2+r2^2 r3^2+r3^2 r1^2-r1^4/2-r2^4/2-r3^4/2)^((ndim-3)/2)];
integ=ce/.{
      E^(-2 a r1-2 b r2-2 c r3) r1^ii_ r2^jj_ r3^kk_->int1[ii-3,jj-3,kk-3],
      E^(-2 a r2-2 b r1-2 c r3) r1^ii_ r2^jj_ r3^kk_->int2[ii-3,jj-3,kk-3],
      E^(-a r1-b r1-a r2-b r2-2 c r3) r1^ii_ r2^jj_ r3^kk_->
        ints[ii-3,jj-3,kk-3]};
T=integ/.{a->a0,b->b0,c->c0};
(* V0 = <1/r1 + 1/r2> *)
ce=Expand[
    psi1 (1/r1+1/r2) psi2 (r1 r2 r3)^3 (
          r1^2 r2^2+r2^2 r3^2+r3^2 r1^2-r1^4/2-r2^4/2-r3^4/2)^((ndim-3)/2)];
integ=ce/.{
      E^(-2 a r1-2 b r2-2 c r3) r1^ii_ r2^jj_ r3^kk_->int1[ii-3,jj-3,kk-3],
      E^(-2 a r2-2 b r1-2 c r3) r1^ii_ r2^jj_ r3^kk_->int2[ii-3,jj-3,kk-3],
      E^(-a r1-b r1-a r2-b r2-2 c r3) r1^ii_ r2^jj_ r3^kk_->
        ints[ii-3,jj-3,kk-3]};
V0=integ/.{a->a0,b->b0,c->c0};
(* V1 = <1/r3> *)
ce=Expand[
    psi1 (1/r3) *
      psi2(r1 r2 r3)^3 (
          r1^2 r2^2+r2^2 r3^2+r3^2 r1^2-r1^4/2-r2^4/2-r3^4/2)^((ndim-3)/2)];
integ=ce/.{
      E^(-2 a r1-2 b r2-2 c r3) r1^ii_ r2^jj_ r3^kk_->int1[ii-3,jj-3,kk-3],
      E^(-2 a r2-2 b r1-2 c r3) r1^ii_ r2^jj_ r3^kk_->int2[ii-3,jj-3,kk-3],
      E^(-a r1-b r1-a r2-b r2-2 c r3) r1^ii_ r2^jj_ r3^kk_->
        ints[ii-3,jj-3,kk-3]};
V1=integ/.{a->a0,b->b0,c->c0};
(* S = <1> *)
ce=Expand[
    psi1   psi2 (r1 r2 r3)^3 (
          r1^2 r2^2+r2^2 r3^2+r3^2 r1^2-r1^4/2-r2^4/2-r3^4/2)^((ndim-3)/2)];
integ=ce/.{
      E^(-2 a r1-2 b r2-2 c r3) r1^ii_ r2^jj_ r3^kk_->int1[ii-3,jj-3,kk-3],
      E^(-2 a r2-2 b r1-2 c r3) r1^ii_ r2^jj_ r3^kk_->int2[ii-3,jj-3,kk-3],
      E^(-a r1-b r1-a r2-b r2-2 c r3) r1^ii_ r2^jj_ r3^kk_->
        ints[ii-3,jj-3,kk-3]};
S=integ/.{a->a0,b->b0,c->c0};
{T,V0,V1,S})

hintegrals[a0_,b0_,c0_,ndim_Integer,nm_Integer]:=
(* Body of the package starts *)(
nmax=2 nm;
nmadd=2 (ndim-3);
nmax3=nmax+nmadd+3;
(* Generating integrals  r1^i r2^j r12^k r1^-1 r2^-1 r12^-1 
    E^(-a r1 -b r2 -c r12)
		from the integral r1^-1 r2^-1 r12^-1 E^(-a r1 -b r2 -c r12) *)
int00=2/((a00+b00) (a00+c00) (b00+c00));
                (* <exp(-a r1 -b r2 -c r12) exp(-a r1 -b r2 -c r12) ...> *)
int=int00/.{a00->2 a0+a t,b00->2 b0+b t,c00->2 c0+ t};
intser=Series[int,{t,0,nmax3}];
intc=CoefficientList[intser,t];
Do[intcc=CoefficientList[intc[[n+1]],{a,b}];
	Do[namax=n-nc;
		Do[nb=namax-na;
			intabc=(-1)^(na+nb+nc) na! nb! nc! intcc[[na+1,nb+1]];
			int1[na-1,nb-1,nc-1]=int2[nb-1,na-1,nc-1]=intabc//Together
		,{na,0,namax}],{nc,0,n}],
	{n,0,nmax3}];
                (* <exp(-a r1 -b r2 -c r12) exp(-b r1 -a r2 -c r12) ...> *)
int=int00/.{a00->a0+b0+a t,b00->a0+b0+b t,c00->2 c0+ t};
intser=Series[int,{t,0,nmax3}];
intc=CoefficientList[intser,t];
Do[intcc=CoefficientList[intc[[n+1]],{a,b}];
	Do[namax=n-nc;namax1=Floor[namax/2];
		Do[nb=namax-na;
			intabc=(-1)^(na+nb+nc) na! nb! nc! intcc[[na+1,nb+1]];
			ints[nb-1,na-1,nc-1]=ints[na-1,nb-1,nc-1]=intabc//Together
		,{na,0,namax1}],{nc,0,n}],
	{n,0,nmax3}];
(* Body of the package ends *)
)

hmapping[nm_Integer,choice_Integer,ind_,mmax_]:=
(* Body of the package starts *)(
mm=0;
If[choice==1,
Do[km=Floor[Sqrt[n]]; (*  i + j^2 + k^2 = n  *)
Do[n1=n-k^2;
jm=Floor[Sqrt[n1]];
Do[i=n1-j^2;
mm=mm+1;
ind[mm]={i,j,k},
{j,0,jm}],{k,0,km}];mmax[n]=mm,{n,0,nm}],
Do[km=n; (*  i + j + k = n  *)
Do[n1=n-k;
jm=n1;
Do[i=n1-j;
mm=mm+1;
ind[mm]={i,j,k},
{j,0,jm}],{k,0,km}];mmax[n]=mm,{n,0,nm}]
];
(* Body of the package ends *)
mm)

End[]

EndPackage[]