###################################################################### # # # MAPLE V programs for analytic evaluation of three-electron # # correlated atomic integral using Slater wave functions # # # # Frank E. Harris, Department of Chemistry, University of Utah # # Salt Lake City, UT 84112 (e-mail harris@dirac.chem.utah.edu) # # # # Reference: paper submitted to Phys. Rev. A. June 1996 by # # Frank E. Harris # # # # These programs last modified November 2, 1996 # # # ###################################################################### # # # Programs included: # # # # AA(a1,a2,a3,a12,a23,a31): Analytically computes integral # # over 3-D variables r1,r2,r3 of (r1*r2*r3*r12*r23*r31)^-1 # # * exp(-a1*r1-a2*r2-a3*r3-a12*r12-a23*r23-a31*r31) # # # # RM(a1,a2,a3): Computes integral AA for a12=a23=a31=0 # # by formula given by Remiddi, Phys. Rev. A 44, 5492 (1991) # # --Eq. (11) of referenced paper. # # # # BH(b1,b2,b12): Analytically computes integral over 3-D # # variables r1,r2 of (r1*r2)^-2 * r12^-1 * exp(-b1*r1-b2*r2 # # -b12*r12); special case of AA in limit of infinite a3 # # See Fromm and Hill, Phys. Rev. A 36, 1013 (1987), Eq. # # (2-47) and (A-2) through (A-6). # # # # BFH(b1,b2,b12): Fromm and Hill's formula (Appendix A) for BH. # # # # DL(x): Computes the real part of the dilogarithm of argument # # x, using transformation if needed to make |x| <= 0.5. # # See Fromm and Hill. # # # # CL(x): Computes Clausen's function for argument x--Eqs. (60) # # and (61) of referenced paper. # # # # CHIX(x): Computes Chi(x)/x, with leading (constant) term # # omitted, for real x--Chi is Euler's chi function-- # # Eq. (72) of referenced paper. # # # # CHIY(y): Computes Chi(iy)/(iy), with leading (constant) # # term omitted, for real y--also Eq. (72). # # # # V(x): Computes V(x), Harris' redefinition of function v(x), # # for real x, using Eqs. (41) and (3) of referenced paper # # # # VN(D,x): Computes V(x), for real x, where D is the product # # of factors comprising (1-|x|)/(1+|x|), or a reduced # # subset thereof, using Eq. (43) of referenced paper # # # # VS(x): Computes v(x)/x, omitting leading (constant) term, # # for real x, by Maclaurin series expansion--Eq. (63) of # # referenced paper. # # # # VSI(y): Computes v(iy)/iy, omitting leading (constant) term, # # for real y, by Maclaurin series expansion--also Eq. (63). # # # # VCL(x): Computes v(x) for real x by series analogous to that # # for Clausen function--Eq. (62) of referenced paper. # # # ###################################################################### Digits := 14; # Change if desire different precision AA := proc(a1,a2,a3,a12,a23,a31) local mu,SS,s,g,G,p,d,D,q,T,X,Y1,Y2,DD, j,i,i1,i2,i3,a; X := 0.1; # "small factors" have absolute value < X Y1 := 0.0125; # Positive sigma^2 < Y1 for "small sigma" Y2 := -0.002; # Negative sigma^2 > Y2 for "small sigma" a := array(0..3,0..3); # a[j,i] is aji mu := array(0..3, 0..3); # mu[j,i] is mu sub i super j g := array(0..3,0..3); # g[j,i] is gamma sub i super j G := array(1..3); # G[j] is Gamma sub j p := array(0..3,0..3); # Factors for making products d := array(0..3,0..3); # d[j,i] is product for g[j,i] D := array(1..3); # D[j] is product for G[j] q := array(1..3); # q[j], Eq. 22 # Place arguments in array a a[0,1]:=a1;a[0,2]:=a2;a[0,3]:=a3;a[1,2]:=a12;a[2,3]:=a23;a[1,3]:=a31; for j from 1 to 3 do for i from 0 to j-1 do a[j,i]:=a[i,j] od od; # Build arrays mu and gamma for j from 0 to 3 do i1:=1; i2:=2; i3:=3; if j=1 then i1:=0 elif j=2 then i2:=0 elif j=3 then i3:=0 fi; mu[j,j] := 2*a[i1,i2]*a[i1,i3]*a[i2,i3]; mu[j,i1] := a[i2,i3]*(a[i1,i2]^2+a[i1,i3]^2-a[j,i1]^2); mu[j,i2] := a[i1,i3]*(a[i2,i1]^2+a[i2,i3]^2-a[j,i2]^2); mu[j,i3] := a[i1,i2]*(a[i3,i1]^2+a[i3,i2]^2-a[j,i3]^2); g[j,j] := sum('mu[j,i]','i'=0..3); for i from 0 to 3 do if i <> j then g[j,i] := g[j,j]-2*(mu[j,j]+mu[j,i]) fi od od; # Build q[j] and G[j] for j = 1..3 for j from 1 to 3 do i1:=1; i2:=2; if j=1 then i1:=3 elif j=2 then i2:=3 fi; q[j] := a[0,i1] + a[0,i2] + a[j,i1] + a[j,i2]; G[j] := (a[0,j]^2 + (a[0,i1]+a[0,i2])*(a[j,i1]+a[j,i2])) * (a[i1,i2]^2 + (a[0,i1]+a[j,i1])*(a[0,i2]+a[j,i2])) / q[j] - q[j] * (a[0,i1]*a[j,i2] + a[0,i2]*a[j,i1]) od; # Make SS=sigma^2 and s=|sigma| SS := a1^2*a23^2*(+a1^2-a2^2-a3^2-a12^2+a23^2-a31^2) + a2^2*a31^2*(-a1^2+a2^2-a3^2-a12^2-a23^2+a31^2) + a3^2*a12^2*(-a1^2-a2^2+a3^2+a12^2-a23^2-a31^2) + a1^2*a2^2*a3^2+a1^2*a12^2*a31^2+a2^2*a23^2*a12^2+a3^2*a31^2*a23^2: s := sqrt(abs(SS)); # print(`sigma-squared, sigma:`,SS,s); # restore if wanted if (SS > Y1) then # This is branch for real sigma not too close to zero # Computed by singularity-cancelling method unless alternate # code (now commented out) at end of this branch is used for j from 0 to 3 do i1:=1; i2:=2; i3:=3; if j=1 then i1:=0 elif j=2 then i2:=0 elif j=3 then i3:=0 fi; p[j,j] := a[j,i1] + a[j,i2] + a[j,i3]; p[j,i1] := -a[j,i1] + a[j,i2] + a[j,i3]; p[j,i2] := a[j,i1] - a[j,i2] + a[j,i3]; p[j,i3] := a[j,i1] + a[j,i2] - a[j,i3]; od; T := 0; # Accumulation cell for small-factor add-on terms # Sets "small factors" less than X to unity and computes add-on for T for j from 0 to 3 do i1:=1; i2:=2; i3:=3; if j=1 then i1:=0 elif j=2 then i2:=0 elif j=3 then i3:=0 fi; if abs(p[j,i1]) < X then T := T + AD(p,q,s,G,g,j,i1,i2,i3); p[j,i1]:=1 fi; if abs(p[j,i2]) < X then T := T + AD(p,q,s,G,g,j,i2,i1,i3); p[j,i2]:=1 fi; if abs(p[j,i3]) < X then T := T + AD(p,q,s,G,g,j,i3,i1,i2); p[j,i3]:=1 fi; od; # Computes products of factors for constructing (1-|x|)/(1+|x|): for j from 0 to 3 do i1:=1; i2:=2; i3:=3; if j=1 then i1:=0 elif j=2 then i2:=0 elif j=3 then i3:=0 fi; d[j,j] := p[i1,i1]*p[i2,i2]*p[i3,i3]*p[i1,j]*p[i2,j]*p[i3,j]/(abs(g[j,j])+s)^2; d[j,i1] := p[i1,i1]*p[i1,j]*p[i2,i1]*p[i3,i1]*p[i2,i3]*p[i3,i2]/(abs(g[j,i1])+s)^2; d[j,i2] := p[i2,i2]*p[i2,j]*p[i1,i2]*p[i3,i2]*p[i1,i3]*p[i3,i1]/(abs(g[j,i2])+s)^2; d[j,i3] := p[i3,i3]*p[i3,j]*p[i1,i3]*p[i2,i3]*p[i1,i2]*p[i2,i1]/(abs(g[j,i3])+s)^2; od; DD:=product('p[j,j]','j'=0..3); for j from 1 to 3 do i1:=1; i2:=2; if j=1 then i1:=3 elif j=2 then i2:=3 fi; D[j] := DD*p[0,j]*p[j,0]*p[i1,i2]*p[i2,i1]/(q[j]*(abs(G[j])+s))^2 od; # If want to use less well conditioned method, comment out the next # two lines and remove comment flags from the two following lines evalf(16*Pi^3*(Pi^2/2 - 2*sum('VN(D[j],G[j]/s)','j'=1..3) + sum('sum('VN(d[j,i],g[j,i]/s)','i'=0..3)','j'=0..3) - T)/s) # evalf(16*Pi^3*(Pi^2/2 - 2*sum('V(G[j]/s)','j'=1..3) # + sum('sum('V(g[j,i]/s)','i'=0..3)','j'=0..3))/s) elif (SS > 0) then # This branch is for sigma real but small DD := -8*a12*a23*a31/(g[0,0]*g[0,1]*g[0,2]*g[0,3]) - 8*a2*a3*a23/(g[1,0]*g[1,1]*g[1,2]*g[1,3]) + 1/(G[1]*g[3,3]*g[2,3]) + 1/(G[1]*g[2,2]*g[3,2]) + 1/(G[2]*g[0,0]*g[2,0]) + 1/(G[2]*g[1,1]*g[3,1]) + 1/(G[3]*g[1,1]*g[2,1]) + 1/(G[3]*g[0,0]*g[3,0]); evalf(16*Pi^3*(-2*sum('FD(G[j],s)','j'=1..3) + 2*SS*DD*(1-log(2*s)) + sum('sum('FD(g[j,i],s)','i'=0..3)','j'=0..3))) elif (SS = 0) then # This branch is for sigma a precise zero evalf(16*Pi^3*(-4*sum('log(abs(G[j]))/G[j]','j'=1..3) + 2*sum('sum('log(abs(g[j,i]))/g[j,i]','i'=0..3)','j'=0..3))) elif (SS > Y2) then # This branch is for small imaginary sigma DD := -8*a12*a23*a31/(g[0,0]*g[0,1]*g[0,2]*g[0,3]) - 8*a2*a3*a23/(g[1,0]*g[1,1]*g[1,2]*g[1,3]) + 1/(G[1]*g[3,3]*g[2,3]) + 1/(G[1]*g[2,2]*g[3,2]) + 1/(G[2]*g[0,0]*g[2,0]) + 1/(G[2]*g[1,1]*g[3,1]) + 1/(G[3]*g[1,1]*g[2,1]) + 1/(G[3]*g[0,0]*g[3,0]); evalf(16*Pi^3*(-2*sum('FDI(G[j],s)','j'=1..3) + 2*SS*DD*(1-log(2*s)) +sum('sum('FDI(g[j,i],s)','i'=0..3)','j'=0..3))) else # This branch is for imaginary sigma not close to zero -16*evalf(Pi^3*(-2*sum('C(G[k]/s)','k'=1..3) + sum('sum('C(g[j,i]/s)','i'=0..3)','j'=0..3))/s) fi end; ########### Procedures C, AD, FD, FDI internal to AA: ################ C := proc(x); CL(evalf(Pi + 2*arctan(x))) end; AD := proc(p,q,s,G,g,j,i1,i2,i3) local J; if abs(p[j,i1]) > 10^(-16) then J := min(i1+j,6-i1-j); sign(G[J])*log(abs(p[j,i1]))*log((q[J]*(abs(G[J])+s)/(q[J]-p[j,i1]))^2 *(abs(g[i2,i1])+s)*(abs(g[i3,i1])+s)/ ((abs(g[i1,j])+s)*(abs(g[i1,i1])+s)*(abs(g[i2,i3])+s)*(abs(g[i3,i2])+s))) else 0 fi end; FD := proc(g,s) local n,y; y := s/g; (1/g)*(2*log(abs(g)) - 2*log(abs(y))*sum(y^(2*n)/(2*n+1),n=1..10) + VS(y) + 2*CHIX(y)) end; FDI := proc(g,s) local n,y; y := s/g; (1/g)*(2*log(abs(g)) - 2*log(abs(y))*sum((-y^2)^n/(2*n+1),n=1..10) + VSI(y) + 2*CHIY(y)) end; ################# Two-electron integral procedures BH, BFH ########## BH := proc(b1,b2,b12) local K12,K1,K2; K12 := evalf((-4*Pi^2/b12)*log((b1-b12)/(b1+b12))*log((b2-b12)/(b2+b12))); K1 := evalf((4*Pi^2/b12)*(DL(-(b2-b12)/(b1+b12))-DL(-(b2+b12)/(b1-b12)))); K2 := evalf((4*Pi^2/b12)*(DL(-(b1-b12)/(b2+b12))-DL(-(b1+b12)/(b2-b12)) +(1/2)*log((b2+b12)/(b1+b12))^2-(1/2)*log((b2-b12)/(b1-b12))^2)); K12+K1+K2; end; BFH := proc(b1,b2,b12) local B12,T; B12 := abs(b12); T := sign(b12)*evalf(+ log((b1+b12)*(b2+b12))^2 - log((b1+b12)*(b2-b12))^2 - log((b1-b12)*(b2+b12))^2 + log((b1-b12)*(b2-b12))^2); evalf((4*Pi^2/B12)*(Pi^2/2 + 2*V(b2/B12) + 2*V(b1/B12) - 2*V((b12^2+b1*b2)/(B12*(b1+b2))) - (1/2)*T)); end; ################# Other procedures listed in preamble ############### ####### Remiddi's formula ####### RM := proc(a1,a2,a3) local w1,w2,w3; w1:=(a2+a3)/a1; w2:=(a1+a3)/a2; w3:=(a1+a2)/a3; evalf((32*Pi^3/(a1*a2*a3))*(log(w1)*log(1+1/w1)-DL(-1/w1)-DL(1-1/w1) +log(w2)*log(1+1/w2)-DL(-1/w2)-DL(1-1/w2) +log(w3)*log(1+1/w3)-DL(-1/w3)-DL(1-1/w3))) end; ####### Dilogarithm ###### DL := proc(x) local w; if (abs(x) <= 1/2) then RETURN(evalf(DLseries(x))) fi; if (x < -1) then w:=1/(1-x); RETURN(evalf(DLseries(w)+log(w)*log(w*x^2)/2-Pi^2/6)) fi; if (x < -1/2) then w:=x/(x-1); RETURN(evalf(-DLseries(w)-log(-w/x)^2/2)) fi; if (evalf(x)=1.) then RETURN(evalf(Pi^2/6)) fi; if (x < 1) then w:=1-x; RETURN(evalf(-DLseries(w) -log(x)*log(w)+Pi^2/6)) fi; if (x <= 2) then w:=1-1/x; RETURN(evalf(DLseries(w)-log(x)*log(w)-log(x)^2/2+Pi^2/6)) fi; evalf(-DLseries(1/x)-log(x)^2/2+Pi^2/3) end; DLseries := proc(x) local n; # series expansion of dilogarithm if (abs(x) > 1/2) then RETURN('procname(args)') fi; sum(x^n/n^2,n=1..30) end; ####### Clausen function ####### CL := proc(y) local n,z; if y <= 2.0944 then y*(1-log(abs(y))+sum((-1)^(n-1)*BB[n]*y^(2*n)/(2*n*(2*n+1)!),n=1..15)) elif y > 4.18879 then -CL(evalf(2*Pi)-y) else z := y-evalf(Pi); z*(-evalf(log(2))+sum((-1)^(n-1)*BB[n]*(2^(2*n)-1)*z^(2*n)/ (2*n*(2*n+1)!),n=1..15)) fi end; ####### Chi function ######## CHIX := proc(y) local n; # Chi, real argument sum(y^(2*n)/(2*n+1)^2,n=1..10) end; CHIY := proc(y) local n; # Chi, imaginary argument sum(y^(2*n)*(-1)^n/(2*n+1)^2,n=1..10) end; ####### V function codes ######## V := proc(X) local x,y,z; x := evalf(X); if not type(x,numeric) then RETURN('procname(args)') fi; if (x = 0) then 0 elif (x < 0) then -V(-x) elif (x > 1) then V(1/x) + DL(1-x)+DL(-x)+evalf(log(x)*log(1+x)+Pi^2/12) elif (evalf(x)=1.) then evalf(-log(2)^2/4-Pi^2/12) else y := (1-x)/2; z := (1+x)/2; (DL(y)-DL(z))/2 - evalf(log(y)^2-log(z)^2)/4 fi end; VN := proc(d,g); sign(g)*evalf(-log(abs(d))^2/4-Pi^2/12+DL((1-abs(g))/2)+log((1+abs(g))/2)^2/2) end; VS := proc(y) local n; sum(CV[n]*y^(2*n-2),n=2..10) end; VSI := proc(y) local n; sum(CV[n]*y^(2*n-2)*(-1)^(n-1),n=2..10) end; VCL := proc(y) local w; w := 2*arctanh(y); -w*(evalf(log(2)) + sum(BB[n]*(2^(2*n)-1)*w^(2*n)/(2*n*(2*n+1)!),n=1..15)) end; ############# Bernoulli numbers and coefficients for VS, VSI ########### BB := array(1..15); # BB[n] is Bernoulli number B(2n) BB[1] := 1/6; BB[2] := -1/30; BB[3] := 1/42; BB[4] := -1/30; BB[5] := 5/66; BB[6] := -691/2730; BB[7] := 7/6; BB[8] := -3617/510; BB[9] := 43867/798; BB[10] := -174611/330; BB[11] := 854513/138; BB[12] := -236364091/2730; BB[13] := 8553103/6; BB[14] := -23749461029/870; BB[15] := 8615841276005/14322; CV:=array(1..10); # Coefficients for Maclaurin series for v(z) CV[1]:=-2*log(2); # CV[n] is coefficient of z^(2n-1) CV[2]:=-(2/3)*log(2) - 1/3; CV[3]:=-(2/5)*log(2) - 3/10; CV[4]:=-(2/7)*log(2) - 11/42; CV[5]:=-(2/9)*log(2) - 25/108; CV[6]:=-(2/11)*log(2) - 137/660; CV[7]:=-(2/13)*log(2) - 49/260; CV[8]:=-(2/15)*log(2) - 121/700; CV[9]:=-(2/17)*log(2) - 761/4760; CV[10]:=-(2/19)*log(2) - 7129/47880; ##################That's all, folks #######################