Hylleraas Three-Electron Integral Recursion Relations

This Mathematica notebook implements the recursion relations for the Hylleraas three-electron integrals commonly used in atomic physics.  The relevant article is by Pachucki, Puchalski and Remiddi, and the DOI is 10.1103/PhysRevA.77.032511.  Their paper can be found at http://link.aps.org/doi/10.1103/PhysRevA.77.032511.  They have also graciously posted their Fortran code at http://www.fuw.edu.pl/~krp.

This notebook was written by Denton Woods (denton.woods@unt.edu) and is available at http://www.dentonwoods.com.  Please email me if you have any questions about using this. The last update to this notebook was 6/7/2014 with the “Alternate Method of Accessing” section.

This notebook is licensed under a Creative Commons Attribution - ShareAlike 3.0 License.  For information on this, please see http://creativecommons.org/licenses/by-sa/3.0/.  If this license does not work for you, contact me, and we should be able to work something out.

Initialization

verbose=false;
(* The recursion relations are analytical, but this controls the digits of the final output value. *)
precision=30;

Γ Definition

Assuming for the moment that we won't have any need for other cases (not using the recursion relations for Γ yet)

I have to do this in a sort of roundabout way where I define a temporary function, since Mathematica gives an error about attempting to take a derivative of something that is not a valid variable (i.e. trying to take the derivative w.r.t. a constant).

Γ[n1_,n2_,n3_,α1_,α2_,α3_]:=Module[{Α1,Α2,Α3,γ1,γ2,γ3},
If[n1≥0&&n2≥0&&n3≥0,
γ1[Α1_,Α2_,Α3_]=Γ1[n1,n2,n3,Α1,Α2,Α3];
ret=γ1[α1,α2,α3];
Return[ret];
];

If[n1==-1&&n2≥0&&n3≥0,
γ2[Α1_,Α2_,Α3_]=Γ2[n1,n2,n3,Α1,Α2,Α3];
ret=γ2[α1,α2,α3];
Return[ret];
];
If[n1≥0&&n2==-1&&n3≥0,
γ2[Α2_,Α1_,Α3_]=Γ2[n2,n1,n3,Α2,Α1,Α3];
ret=γ2[α2,α1,α3];
Return[ret];
];
If[n1≥0&&n2≥0&&n3==-1,
γ2[Α3_,Α2_,Α1_]=Γ2[n3,n2,n1,Α3,Α2,Α1];
ret=γ2[α3,α2,α1];
Return[ret];
];

If[n1==-1&&n2==-1&&n3≥0,
γ3[Α1_,Α2_,Α3_]=Γ3[n1,n2,n3,Α1,Α2,Α3];
ret=γ3[α1,α2,α3];
Return[ret];
];
If[n1≥0&&n2==-1&&n3==-1,
γ3[Α1_,Α2_,Α3_]=Γ3[n2,n3,n1,Α3,Α2,Α1];
ret=γ3[α1,α2,α3];
Return[ret];
];
If[n1==-1&&n2≥0&&n3==-1,
γ3[Α1_,Α2_,Α3_]=Γ3[n1,n3,n2,Α1,Α3,Α2];
ret=γ3[α1,α2,α3];
Return[ret];
];

If[verbose==true,Print["Input values out of range for Γ function: ", n1, " ", n2," ",n3]];
Return[1];
]

(* Used when l,m,n ≥ 0 *)
Γ1[n1_,n2_,n3_,α1_,α2_,α3_]:=(((-1))^(n1+n2+n3))(D[D[D[(1/((α1+α2)(α2+α3)(α3+α1))),{α1,n1}],{α2,n2}],{α3,n3}])

(* Used when we have -1,m,n *)
Γ2[n1_,n2_,n3_,α1_,α2_,α3_]:=Module[{r},
If[n2==0&&n3==0&&α2==α3,Return[(1/((α1+α2)α2))]];
r=(((-1))^(n2+n3))(D[D[((Log[α1+α2]-Log[α1+α3])/((α2-α3)(α2+α3))),{α2,n2}],{α3,n3}]);
Return[r];
]

(* Used when we have -1,-1,n3 *)
Γ3[n1_,n2_,n3_,α1_,α2_,α3_]:=(1/(4α3))(((-1))^n3)D[((Log[((α1+α3)/(α2+α3))])^2)+2(PolyLog[2,((α1+α2)/(α1+α3))]+PolyLog[2,((α1+α2)/(α2+α3))])+((π^2)/3),{α3,n3}]

Testing Γ definition for some cases

Γ[0,0,0,1,1/2,3/2]
Γ[-1,0,0,1,1/2,3/2]
Γ[-2,0,0,1,1/2,3/2]
Γ[0,-1,0,1,1/2,3/2]
Γ[0,0,-1,1,1/2,3/2]

Redefining Γ1 in terms of recursion relations (Korobov paper)

Γ000[α_,β_,γ_]:=(2/((α+β)(β+γ)(γ+α)))
Γ1rec[l_,m_,n_,α_,β_,γ_]:=Module[{},
(* This shouldn't happen, because we have restricted the calls to Γ1 to be non-negative lmn. *)
If[l<0||m<0||n<0,Return[0];];
r=(1/(α+β))(l*Γ1rec[l-1,m,n,α,β,γ]+m Γ1rec[l,m-1,n,α,β,γ]+B[l,m,n,α,β,γ]);
Return[r];
];

B[l_,m_,n_,α_,β_,γ_]:=Module[{},
r=(1/(α+γ))(l*B[l-1,m,n,α,β,γ]+n*B[l,m,n-1,α,β,γ]+A[l-1,m,n,α,β,γ]);
Return[r];
];

A[l_,m_,n_,α_,β_,γ_]:=δ[l]((2(m+n)!)/(((β+γ))^(m+n+1)))

Defining the Starred Functions in Equation 8

f1star[n2_,n3_,n4_,n5_,n6_,w1_,w2_,w3_]:=Module[{},
Return[Γ[n5+n6-1,n4,n3+n2-1,w2+w3,w1,0]];
]
f2star[n1_,n3_,n4_,n5_,n6_,w1_,w2_,w3_]:=Module[{},
Return[Γ[n4+n6-1,n5,n1+n3-1,w1+w3,w2,0]];
]
f3star[n1_,n2_,n4_,n5_,n6_,w1_,w2_,w3_]:=Module[{},
Return[Γ[n4+n5-1,n6,n1+n2-1,w1+w2,w3,0]];
]

The next three are special cases, because the n4, n5 and n6 must be 0 (omitting the starred term)

f4star[n1_,n2_,n3_,n5_,n6_,w1_,w2_,w3_]:=Module[{},
If[n5≠0||n6≠0,Print["Error in input values for f4*"];Return[0];];
Return[Γ[n3-1,n2-1,n1,w2,w3,0]];
]
f5star[n1_,n2_,n3_,n4_,n6_,w1_,w2_,w3_]:=Module[{},
If[n4≠0||n6≠0,Print["Error in input values for f5*"];Return[0];];
Return[Γ[n1-1,n3-1,n2,w3,w1,0]];
]
f6star[n1_,n2_,n3_,n4_,n5_,w1_,w2_,w3_]:=Module[{},
If[n4≠0||n5≠0,Print["Error in input values for f6*"];Return[0];];
Return[Γ[n2-1,n1-1,n3,w1,w2,0]];
]

Defining the Initial Functions in Equation 7

The number next to the f shows which of the parameters is equal to 1 (f0 has all parameters equal to 0).

f0[w1_,w2_,w3_]:=-(1/(2 w1 w2 w3))(Log[(w3/(w1+w2))]Log[1+(w3/(w1+w2))]+PolyLog[2,-(w3/(w1+w2))]+PolyLog[2,1-(w3/(w1+w2))]+Log[(w2/(w3+w1))]Log[1+(w2/(w1+w3))]+PolyLog[2,-(w2/(w1+w3))]+PolyLog[2,1-(w2/(w1+w3))]+Log[(w1/(w2+w3))]Log[1+(w1/(w2+w3))]+PolyLog[2,-(w1/(w2+w3))]+PolyLog[2,1-(w1/(w2+w3))])

f1[w1_,w2_,w3_]:=-(1/((w2^2)(w3^2)))Log[((w1(w1+w2+w3))/((w1+w2)(w1+w3)))]
f2[w1_,w2_,w3_]:=-(1/((w1^2)(w3^2)))Log[((w2(w1+w2+w3))/((w2+w3)(w2+w1)))]
f3[w1_,w2_,w3_]:=-(1/((w1^2)(w2^2)))Log[((w3(w1+w2+w3))/((w3+w1)(w3+w2)))]

f12[w1_,w2_,w3_]:=(1/(w1 w2(w1+w2)(w3^2)))
f13[w1_,w2_,w3_]:=(1/(w1 w3(w1+w3)(w2^2)))
f23[w1_,w2_,w3_]:=(1/(w2 w3(w2+w3)(w1^2)))

f123[w1_,w2_,w3_]:=(1/((w1^2)(w2^2)(w3^2)))

First Set of Recursion Relations

δ[x_]:=KroneckerDelta[x,0]

fn1plus2[n11_,n2_,n3_,n4_,n5_,n6_,w1_,w2_,w3_]:=Module[{n1},
n1=n11-2;
r=((1+n1)/2)((1/(w3^2))(((n2(n2-1))/(n1+1))f[n1+2,n2-2,n3,0,0,0,w1,w2,w3]+(n3+2n2+n1+2)f[n1,n2,n3,0,0,0,w1,w2,w3]+(1/(n1+1))f6star[n1+2,n2,n3,0,0,w1,w2,w3]+((n2(n2-1))/(n3+1))f[n1,n2-2,n3+2,0,0,0,w1,w2,w3]+(1/(n3+1))f4star[n1,n2,n3+2,0,0,w1,w2,w3]-((δ[n2])/(n1+1))f2star[n1+2,n3,0,0,0,w1,w2,w3]-((δ[n2])/(n3+1))f2star[n1,n3+2,0,0,0,w1,w2,w3])
+(1/(w2^2))(((n3(n3-1))/(n1+1))f[n1+2,n2,n3-2,0,0,0,w1,w2,w3]+(2n3+n2+n1+2)f[n1,n2,n3,0,0,0,w1,w2,w3]+(1/(n1+1))f5star[n1+2,n2,n3,0,0,w1,w2,w3]+((n3(n3-1))/(n2+1))f[n1,n2+2,n3-2,0,0,0,w1,w2,w3]+(1/(n2+1))f4star[n1,n2+2,n3,0,0,w1,w2,w3]-((δ[n3])/(n1+1))f3star[n1+2,n2,0,0,0,w1,w2,w3]-((δ[n3])/(n2+1))f3star[n1,n2+2,0,0,0,w1,w2,w3])
-((w1^2)/((w3^2)(w2^2)))(((n1(n1-1))/(n2+1))f[n1-2,n2+2,n3,0,0,0,w1,w2,w3]+(n3+n2+2n1+2)f[n1,n2,n3,0,0,0,w1,w2,w3]+(1/(n2+1))f6star[n1,n2+2,n3,0,0,w1,w2,w3]+((n1(n1-1))/(n3+1))f[n1-2,n2,n3+2,0,0,0,w1,w2,w3]+(1/(n3+1))f5star[n1,n2,n3+2,0,0,w1,w2,w3]-((δ[n1])/(n2+1))f1star[n2+2,n3,0,0,0,w1,w2,w3]-((δ[n1])/(n3+1))f1star[n2,n3+2,0,0,0,w1,w2,w3]));
Return[r];
]

fn2plus2[n1_,n21_,n3_,n4_,n5_,n6_,w1_,w2_,w3_]:=Module[{n2},
n2=n21-2;
r=((1+n2)/2)((1/(w1^2))(((n3(n3-1))/(n2+1))f[n1,n2+2,n3-2,0,0,0,w1,w2,w3]+(n1+2n3+n2+2)f[n1,n2,n3,0,0,0,w1,w2,w3]+(1/(n2+1))f4star[n1,n2+2,n3,0,0,w1,w2,w3]+((n3(n3-1))/(n1+1))f[n1+2,n2,n3-2,0,0,0,w1,w2,w3]+(1/(n1+1))f5star[n1+2,n2,n3,0,0,w1,w2,w3]-((δ[n3])/(n2+1))f3star[n1,n2+2,0,0,0,w1,w2,w3]-((δ[n3])/(n1+1))f3star[n1+2,n2,0,0,0,w1,w2,w3])
+(1/(w3^2))(((n1(n1-1))/(n2+1))f[n1-2,n2+2,n3,0,0,0,w1,w2,w3]+(2n1+n3+n2+2)f[n1,n2,n3,0,0,0,w1,w2,w3]+(1/(n2+1))f6star[n1,n2+2,n3,0,0,w1,w2,w3]+((n1(n1-1))/(n3+1))f[n1-2,n2,n3+2,0,0,0,w1,w2,w3]+(1/(n3+1))f5star[n1,n2,n3+2,0,0,w1,w2,w3]-((δ[n1])/(n2+1))f1star[n2+2,n3,0,0,0,w1,w2,w3]-((δ[n1])/(n3+1))f1star[n2,n3+2,0,0,0,w1,w2,w3])
-((w2^2)/((w1^2)(w3^2)))(((n2(n2-1))/(n3+1))f[n1,n2-2,n3+2,0,0,0,w1,w2,w3]+(n1+n3+2n2+2)f[n1,n2,n3,0,0,0,w1,w2,w3]+(1/(n3+1))f4star[n1,n2,n3+2,0,0,w1,w2,w3]+((n2(n2-1))/(n1+1))f[n1+2,n2-2,n3,0,0,0,w1,w2,w3]+(1/(n1+1))f6star[n1+2,n2,n3,0,0,w1,w2,w3]-((δ[n2])/(n3+1))f2star[n1,n3+2,0,0,0,w1,w2,w3]-((δ[n2])/(n1+1))f2star[n1+2,n3,0,0,0,w1,w2,w3]));
Return[r];
]

fn3plus2[n1_,n2_,n31_,n4_,n5_,n6_,w1_,w2_,w3_]:=Module[{n3},
n3=n31-2;
r=((1+n3)/2)((1/(w1^2))(((n2(n2-1))/(n3+1))f[n1,n2-2,n3+2,0,0,0,w1,w2,w3]+(n1+2n2+n3+2)f[n1,n2,n3,0,0,0,w1,w2,w3]+(1/(n3+1))f4star[n1,n2,n3+2,0,0,w1,w2,w3]+((n2(n2-1))/(n1+1))f[n1+2,n2-2,n3,0,0,0,w1,w2,w3]+(1/(n1+1))f6star[n1+2,n2,n3,0,0,w1,w2,w3]-((δ[n2])/(n3+1))f2star[n1,n3+2,0,0,0,w1,w2,w3]-((δ[n2])/(n1+1))f2star[n1+2,n3,0,0,0,w1,w2,w3])
+(1/(w2^2))(((n1(n1-1))/(n3+1))f[n1-2,n2,n3+2,0,0,0,w1,w2,w3]+(2n1+n2+n3+2)f[n1,n2,n3,0,0,0,w1,w2,w3]+(1/(n3+1))f5star[n1,n2,n3+2,0,0,w1,w2,w3]+((n1(n1-1))/(n2+1))f[n1-2,n2+2,n3,0,0,0,w1,w2,w3]+(1/(n2+1))f6star[n1,n2+2,n3,0,0,w1,w2,w3]-((δ[n1])/(n3+1))f1star[n2,n3+2,0,0,0,w1,w2,w3]-((δ[n1])/(n2+1))f1star[n2+2,n3,0,0,0,w1,w2,w3])
-((w3^2)/((w1^2)(w2^2)))(((n3(n3-1))/(n2+1))f[n1,n2+2,n3-2,0,0,0,w1,w2,w3]+(n1+n2+2n3+2)f[n1,n2,n3,0,0,0,w1,w2,w3]+(1/(n2+1))f4star[n1,n2+2,n3,0,0,w1,w2,w3]+((n3(n3-1))/(n1+1))f[n1+2,n2,n3-2,0,0,0,w1,w2,w3]+(1/(n1+1))f5star[n1+2,n2,n3,0,0,w1,w2,w3]-((δ[n3])/(n2+1))f3star[n1,n2+2,0,0,0,w1,w2,w3]-((δ[n3])/(n1+1))f3star[n1+2,n2,0,0,0,w1,w2,w3]));
Return[r];
]

Second Set of Recursion Relations

From the fn6plus1, fn4plus1 is given by 1↔3 and 4↔6.  I am unsure about whether we have to switch the inputs to the f functions for w1, w2 and w3.  They are not switched right now.

fn4plus1[n1_,n2_,n3_,n41_,n5_,n6_,w1_,w2_,w3_]:=Module[{n4},
n4=n41-1;
r=(1/(w3 w2 w1))((n3-1)n3 n6 f[n1,n2,n3-2,n4,n5+1,n6-1,w1,w2,w3]+(n2-1)n2 n5 f[n1,n2-2,n3,n4,n5-1,n6+1,w1,w2,w3]-(n1-1)n1 n6 f[n1-2,n2,n3,n4,n5+1,n6-1,w1,w2,w3]-(n1-1)n1 n5 f[n1-2,n2,n3,n4,n5-1,n6+1,w1,w2,w3]+n6 n5 (n3+n2-n1+n4+1)f[n1,n2,n3,n4,n5-1,n6-1,w1,w2,w3]-(n3-1)n3 f[n1,n2,n3-2,n4,n5+1,n6,w1,w2,w3]w3+(n1-1)n1 f[n1-2,n2,n3,n4,n5+1,n6,w1,w2,w3]w3-n5(n3+n2-n1+n4+1)f[n1,n2,n3,n4,n5-1,n6,w1,w2,w3]w3-(n2-1)n2 f[n1,n2-2,n3,n4,n5,n6+1,w1,w2,w3]w2+(n1-1)n1 f[n1-2,n2,n3,n4,n5,n6+1,w1,w2,w3]w2-n6(n3+n2-n1+n4+1)f[n1,n2,n3,n4,n5,n6-1,w1,w2,w3]w2-n6 n5 f[n1,n2,n3,n4+1,n5-1,n6-1,w1,w2,w3]w1+(n3+n2-n1+n4+1)f[n1,n2,n3,n4,n5,n6,w1,w2,w3]w3 w2+n6 f[n1,n2,n3,n4+1,n5,n6-1,w1,w2,w3]w2 w1+n5 f[n1,n2,n3,n4+1,n5-1,n6,w1,w2,w3]w1 w3-δ[n3]n6 f3star[n1,n2,n4,n5+1,n6-1,w1,w2,w3]+δ[n3] f3star[n1,n2,n4,n5+1,n6,w1,w2,w3]w3-δ[n2] n5 f2star[n1,n3,n4,n5-1,n6+1,w1,w2,w3]+δ[n2]f2star[n1,n3,n4,n5,n6+1,w1,w2,w3]w2+δ[n1]n6 f1star[n2,n3,n4,n5+1,n6-1,w1,w2,w3]+δ[n1]n5 f1star[n2,n3,n4,n5-1,n6+1,w1,w2,w3]-δ[n1]f1star[n2,n3,n4,n5+1,n6,w1,w2,w3]w3-δ[n1]f1star[n2,n3,n4,n5,n6+1,w1,w2,w3]w2);
Return[r];
]

From the fn6plus1, fn5plus1 is given by 2↔3 and 5↔6.  I am unsure about whether we have to switch the inputs to the f functions for w1, w2 and w3.  They are not switched right now.

fn5plus1[n1_,n2_,n3_,n4_,n51_,n6_,w1_,w2_,w3_]:=Module[{n5},
n5=n51-1;
r=(1/(w1 w3 w2))((n1-1)n1 n4 f[n1-2,n2,n3,n4-1,n5,n6+1,w1,w2,w3]+(n3-1)n3 n6 f[n1,n2,n3-2,n4+1,n5,n6-1,w1,w2,w3]-(n2-1)n2 n4 f[n1,n2-2,n3,n4-1,n5,n6+1,w1,w2,w3]-(n2-1)n2 n6 f[n1,n2-2,n3,n4+1,n5,n6-1,w1,w2,w3]+n4 n6 (n1+n3-n2+n5+1)f[n1,n2,n3,n4-1,n5,n6-1,w1,w2,w3]-(n1-1)n1 f[n1-2,n2,n3,n4,n5,n6+1,w1,w2,w3]w1+(n2-1)n2 f[n1,n2-2,n3,n4,n5,n6+1,w1,w2,w3]w1-n6(n1+n3-n2+n5+1)f[n1,n2,n3,n4,n5,n6-1,w1,w2,w3]w1-(n3-1)n3 f[n1,n2,n3-2,n4+1,n5,n6,w1,w2,w3]w3+(n2-1)n2 f[n1,n2-2,n3,n4+1,n5,n6,w1,w2,w3]w3-n4(n1+n3-n2+n5+1)f[n1,n2,n3,n4-1,n5,n6,w1,w2,w3]w3-n4 n6 f[n1,n2,n3,n4-1,n5+1,n6-1,w1,w2,w3]w2+(n1+n3-n2+n5+1)f[n1,n2,n3,n4,n5,n6,w1,w2,w3]w1 w3+n4 f[n1,n2,n3,n4-1,n5+1,n6,w1,w2,w3]w3 w2+n6 f[n1,n2,n3,n4,n5+1,n6-1,w1,w2,w3]w2 w1-δ[n1]n4 f1star[n2,n3,n4-1,n5,n6+1,w1,w2,w3]+δ[n1] f1star[n2,n3,n4,n5,n6+1,w1,w2,w3]w1-δ[n3] n6 f3star[n1,n2,n4+1,n5,n6-1,w1,w2,w3]+δ[n3]f3star[n1,n2,n4+1,n5,n6,w1,w2,w3]w3+δ[n2]n4 f2star[n1,n3,n4-1,n5,n6+1,w1,w2,w3]+δ[n2]n6 f2star[n1,n3,n4+1,n5,n6-1,w1,w2,w3]-δ[n2]f2star[n1,n3,n4,n5,n6+1,w1,w2,w3]w1-δ[n2]f2star[n1,n3,n4+1,n5,n6,w1,w2,w3]w3);
Return[r];
]

This is given directly in the paper as equation 5.

fn6plus1[n1_,n2_,n3_,n4_,n5_,n61_,w1_,w2_,w3_]:=Module[{n6},
n6=n61-1;
r=(1/(w1 w2 w3))((n1-1)n1 n4 f[n1-2,n2,n3,n4-1,n5+1,n6,w1,w2,w3]+(n2-1)n2 n5 f[n1,n2-2,n3,n4+1,n5-1,n6,w1,w2,w3]-(n3-1)n3 n4 f[n1,n2,n3-2,n4-1,n5+1,n6,w1,w2,w3]-(n3-1)n3 n5 f[n1,n2,n3-2,n4+1,n5-1,n6,w1,w2,w3]+n4 n5 (n1+n2-n3+n6+1)f[n1,n2,n3,n4-1,n5-1,n6,w1,w2,w3]-(n1-1)n1 f[n1-2,n2,n3,n4,n5+1,n6,w1,w2,w3]w1+(n3-1)n3 f[n1,n2,n3-2,n4,n5+1,n6,w1,w2,w3]w1-n5(n1+n2-n3+n6+1)f[n1,n2,n3,n4,n5-1,n6,w1,w2,w3]w1-(n2-1)n2 f[n1,n2-2,n3,n4+1,n5,n6,w1,w2,w3]w2+(n3-1)n3 f[n1,n2,n3-2,n4+1,n5,n6,w1,w2,w3]w2-n4(n1+n2-n3+n6+1)f[n1,n2,n3,n4-1,n5,n6,w1,w2,w3]w2-n4 n5 f[n1,n2,n3,n4-1,n5-1,n6+1,w1,w2,w3]w3+(n1+n2-n3+n6+1)f[n1,n2,n3,n4,n5,n6,w1,w2,w3]w1 w2+n4 f[n1,n2,n3,n4-1,n5,n6+1,w1,w2,w3]w2 w3+n5 f[n1,n2,n3,n4,n5-1,n6+1,w1,w2,w3]w3 w1-δ[n1]n4 f1star[n2,n3,n4-1,n5+1,n6,w1,w2,w3]+δ[n1] f1star[n2,n3,n4,n5+1,n6,w1,w2,w3]w1-δ[n2] n5 f2star[n1,n3,n4+1,n5-1,n6,w1,w2,w3]+δ[n2]f2star[n1,n3,n4+1,n5,n6,w1,w2,w3]w2+δ[n3]n4 f3star[n1,n2,n4-1,n5+1,n6,w1,w2,w3]+δ[n3]n5 f3star[n1,n2,n4+1,n5-1,n6,w1,w2,w3]-δ[n3]f3star[n1,n2,n4,n5+1,n6,w1,w2,w3]w1-δ[n3]f3star[n1,n2,n4+1,n5,n6,w1,w2,w3]w2);
Return[r];
]

Calculate the Integral

Probably don't need to even check whether n4 == Max[...], etc.  Just check whether they are greater than 0 (or 1 for n1, n2 and n3).

f[n1_,n2_,n3_,n4_,n5_,n6_,w1_,w2_,w3_]:=Module[{},
If[n1<0||n2<0||n3<0||n4<0||n5<0||n6<0,
If[verbose==true,Print["Input to f out of range: ",n1," ",n2," ",n3," ",n4," ",n5," ",n6," "]];
(* Doesn't really matter what the output here is, because a coefficient will always be 0 in this case *)
Return[1]];
(* At the end of the recursion - initial terms *)
If[n1==0&&n2==0&&n3==0&&n4==0&&n5==0&&n6==0,Return[f0[w1,w2,w3]]];
If[n1==1&&n2==0&&n3==0&&n4==0&&n5==0&&n6==0,Return[f1[w1,w2,w3]]];
If[n1==0&&n2==1&&n3==0&&n4==0&&n5==0&&n6==0,Return[f2[w1,w2,w3]]];
If[n1==0&&n2==0&&n3==1&&n4==0&&n5==0&&n6==0,Return[f3[w1,w2,w3]]];
If[n1==1&&n2==1&&n3==0&&n4==0&&n5==0&&n6==0,Return[f12[w1,w2,w3]]];
If[n1==1&&n2==0&&n3==1&&n4==0&&n5==0&&n6==0,Return[f13[w1,w2,w3]]];
If[n1==0&&n2==1&&n3==1&&n4==0&&n5==0&&n6==0,Return[f23[w1,w2,w3]]];
If[n1==1&&n2==1&&n3==1&&n4==0&&n5==0&&n6==0,Return[f123[w1,w2,w3]]];
(* Reduce n4, n5 and n6 down to 0 first *)
If[Max[n4,n5,n6]>0&&n4==Max[n4,n5,n6],Return[fn4plus1[n1,n2,n3,n4,n5,n6,w1,w2,w3]]];
If[Max[n4,n5,n6]>0&&n5==Max[n4,n5,n6],Return[fn5plus1[n1,n2,n3,n4,n5,n6,w1,w2,w3]]];
If[Max[n4,n5,n6]>0&&n6==Max[n4,n5,n6],Return[fn6plus1[n1,n2,n3,n4,n5,n6,w1,w2,w3]]];
(* Then reduce n1, n2 and n3 down to the initial terms *)
If[Max[n1,n2,n3]>1&&n1==Max[n1,n2,n3],Return[fn1plus2[n1,n2,n3,n4,n5,n6,w1,w2,w3]]];
If[Max[n1,n2,n3]>1&&n2==Max[n1,n2,n3],Return[fn2plus2[n1,n2,n3,n4,n5,n6,w1,w2,w3]]];
If[Max[n1,n2,n3]>1&&n3==Max[n1,n2,n3],Return[fn3plus2[n1,n2,n3,n4,n5,n6,w1,w2,w3]]];
(* Shouldn't get here - could print a warning message. *)
Return[0];
]

Tests

Testing first set with n = 0

f[0,0,0,0,0,0,1,1,1]
N[%,precision]

Testing first set with n = 1

f[1,0,0,0,0,0,1,1,1]
N[%,precision]

Testing second set with n = 0

f[0,0,0,0,0,0,1,1,1]
N[%,precision]

Testing second set with n = 1

f[0,0,0,1,0,0,1,1,1]
N[%,precision]

Testing first set with n = 2

f[2,0,0,0,0,0,1,1,1]
N[%,precision]

Testing n4 = 2

f[0,0,0,2,0,0,1,1,1]
N[%,precision]

Testing n3 = 2

f[0,0,2,0,0,0,(5/10),(6/10),1]
N[%,precision]

Recreating Table I in their paper

For[n=0,n≤9,n++,Print[n,": ",N[f[n,0,0,0,0,0,1,1,1],precision ]]]

For[n=0,n≤9,n++,Print[n,": ",N[f[0,0,0,n,0,0,1,1,1],precision ]]]

Alternate Method of Accessing

Instead of calling the f function with the typical list of parameters, we can pass the integrand to HylleraasIntegral. This gives a much nicer and easier way to view exactly what we are doing here. The HylleraasIntegral and HylleraasIntegralSub functions aren’t pretty, and there is likely to be a better way to do this. The integrand must be passed in the expected form, or else there will be no evaluation of this integral.

HylleraasIntegral[Exp[α_ r1+ β_ r2+ γ_ r3]c_]:=Module[{d=c},
If[α≥0||β≥0||γ≥0,Print["-α, -β and -γ must be < 0"];Return[0];];
(* The pattern matching does not match against 1 ((r^0)) for missing polynomials, so we fake it by choosing -100 as the power. The powers are nowhere near this singular, so HylleraasIntegralSub can recognize this and change them to 0. *)
If[MatchQ[c,g_ r1^j1_.]==False,d=d/(r1^100)];
If[MatchQ[c,g_ r2^j2_.]==False,d=d/(r2^100)];
If[MatchQ[c,g_ r3^j3_.]==False,d=d/(r3^100)];
If[MatchQ[c,g_ r12^j12_.]==False,d=d/(r12^100)];
If[MatchQ[c,g_ r23^j23_.]==False,d=d/(r23^100)];
If[MatchQ[c,g_ r31^j31_.]==False,d=d/(r31^100)];
HylleraasIntegralSub[Exp[α r1+ β r2+ γ r3]d]
];

HylleraasIntegralSub[Exp[α_ r1+ β_ r2+ γ_ r3] (r1^(j1_.))(r2^(j2_.))(r3^(j3_.))(r12^(j12_.))(r23^(j23_.))(r31^(j31_.)) ]:=Module[{j1c=j1,j2c=j2,j3c=j3,j12c=j12,j23c=j23,j31c=j31},
(* Any (Subscript[r, i]) or (Subscript[r, ij]) with power of -100 is actually 0, as mentioned in the HylleraasIntegral comment. *)
If[j1c==-100,j1c=0];
If[j2c==-100,j2c=0];
If[j3c==-100,j3c=0];
If[j12c==-100,j12c=0];
If[j23c==-100,j23c=0];
If[j31c==-100,j31c=0];
(*Print[-α," ",-β," ",-γ,": ", j1c," ",j2c," ",j3c," ",j12c," ",j23c," ",j31c];*)
N[f[j23c+1,j31c+1,j12c+1,j1c+1,j2c+1,j3c+1,-α,-β,-γ],precision]
]

Recreating Table I

For[n=0,n≤9,n++,Print[n,": ",HylleraasIntegral[(E^(-r1-r2-r3))(r23^(n-1))(r31^(-1))(r12^(-1))(r1^(-1))(r2^(-1))(r3^(-1))]]]

For[n=0,n≤9,n++,Print[n,": ",HylleraasIntegral[(E^(-r1-r2-r3))(r23^(-1))(r31^(-1))(r12^(-1))(r1^(n-1))(r2^(-1))(r3^(-1))]]]