Potential functions

Quantum Mechanics
                     K. Konishi, G. Paffuti

Write a Mathematica procedure to compute the effective potentials appearing in the Hartree and Hartree-Fock equations.

ClearAll["Global`*"];
Off[General::spell, General::spell1];

The problem

In solving Hartree-Fock  equations effective potentials as

(Y^k)[a,b,r] = (∫_0^r)((((s/r)))^k)P[a,s]P[b,s]ds + (∫_r^∞)((((r/s)))^(k+1))P[a,s]P[b,s]ds

must be computed very often. P[a,r] are the reduced radial wave functions for orbital a,b etc.

The computation of the integrals is simplified by a method due to Hartree. It is convenient to introduce an auxiliary function (Z^k)and transform the original problem in a couple of differential equations:

(d/dr)(Z^k)[a,b,r]= P[a,r]P[b,r]- (k/r)(Z^k)[a,b,r];
(d/dr)(Y^k)[a,b,r]=(1/r)( (k+1)(Y^k)[a,b,r] - (2k+1)(Z^k)[a,b,r] );
(Z^k)[a,b,0]=0;   (Y^k)[a,b,r] -> (Z^k)[a,b,r]    per    r -> ∞

or:

(d/dr)(r^k) (Z^k)[a,b,r]=(r^k) P[a,r]P[b,r];
(d/dr)(((Y^k)[a,b,r])/(r^(k+1)))= - ((2k+1)/(r^(k+2)))(Z^k)[a,b,r];

We leave the (easy) proof to the reader.  It is necessary to take the derivative of (), substitute in the result the definition of Y and use the solution of eq.() for Z.

To have an accurate solution it is useful to exploit the behavior of the wave functions for r→0.

It can be proved (see below) that the functions (Z^k) as r→0 have the expansion ( Z in the rhs is the atomic number)

(Z^k)[a,b,r] ≃ (1/((Subscript[L, a])+(Subscript[L, b])+k+3))r P[a,r]P[b,r](1 + ((Z r)/((Subscript[L, a])+(Subscript[L, b])+k+4))((1/((Subscript[L, a])+1))+(1/((Subscript[L, b])+1))))

(Subscript[L, a]) and (Subscript[L, b]) are the orbital angular momenta of the orbitals.

Proof

For small r reduced wave functions behave as

f[a_,r_,L_]:=c1[a] (r^(L+1))+c2[a](r^(L+2))+c3[a](r^(L+3))

Acting with the kinetic and the coulombic part of the Schrödinger Hamiltonian one has

ris1=Simplify[-(1/2)((∂_r)((∂_r) f[a,r,L])) - (Z/r)f[a,r,L] + ((L(L+1))/(2 (r^2)))f[a,r,L]]

The rhs of the Schrödinger behaves as E f[r] ∼ r^(L+1), then the term proportional to r^L must vanish in the lhs.

ris2=Normal[Series[-ris1/r^L,{r,0,0}]];
sost1=First[Solve[ris2== 0,c2[a]]]

The coefficient c2 is fixed by the Schrödinger equation and the expansion takes the form

g[a_,r_,s_]:=c1[a] r^(s+1)-((Z c1[a])/(1+s))r^(s+2)

Let us compute for small values of r integrals of g[a] g[b] r^k

ris3 = Integrate[g[a,x,sa]g[b,x,sb]*(x^k),{x,0,r},
Assumptions-> {sa>0,sb>0, k>0} ];

At next to leading order

ris4=Normal[Series[ris3,{r,0,4}]]

In the result we can express our coefficients in terms of the behavior at the origin of the wave functions

sost2=First[Solve[{g[a,r,sa] == Pa,g[b,r,sb] == Pb},{c1[a],c1[b]}]];
sost3=Simplify[Normal[Series[sost2,{r,0,1}]]]

The result is finally :

risA=Simplify[((ris4/r^k)/.sost3)];
risA = Normal[Series[risA,{r,0,2}]]

This result is the same as the one written in ()

risIneq4 = (1/(sa+sb+k+3))r Pa Pb(1 + ((Z r)/(sa+sb+k+4))((1/(sa+1))+(1/(sb+1))));

Simplify[risIneq4-risA== 0]

Numerical Integration

The couple of equations () can be directly integrated. To perform numerical integrations we will use (see also notebook [NB7] )

(∫_(Subscript[x, j-1])^(Subscript[x, j+1]))f[x]dx = 2h(f[j] +(1/6)(δ^2)f[j]-(1/180)(δ^4)f[j])+O((h^7))=
                                (h/90)(-f[j+2] + 34 f[j+1] + 114 f[j] + 34 f[j-] - f[j-2])+O((h^7));

In the first interval we adopt the usual Simpson rule

(∫_(Subscript[x, 0])^(Subscript[x, 2]))f[x]dx = (h/3)(f[0]+4 f[1] + f[2])

Explicitly the form of the discretized form is

s[j_]:= (1/90)(-f[j+2]+34f[j+1]+114 f[j] + 34 f[j-1]-f[j-2]);
Expand[Simplify[(1/3)(f[0]+4f[1]+f[2])+Sum[s[k],{k,3,30,2}]]]

ClearAll["Global`*"];

Recurrence relations

Consider now the differential equations ().

First equation

In the first equation let us call (Subscript[z, j]) the value of (r^k) (Z^k) on grid points, and let f be the function in the rhs. The integration formula () allows to write a recursion formula for even and odd sites, separately:

(Subscript[z, j]) - (Subscript[z, j-2]) =  (h/90)(-(Subscript[f, j+1]) + 34 (Subscript[f, j]) +114 (Subscript[f, j-1]) + 34 (Subscript[f, j-2]) - (Subscript[f, j-3])).

For the first and last intervals the simplest Simpson formula is used. In the following we assume an even number of internal grid points.

The asymptotic condition can be refined for k=0 by imposing the asymptotic value  Z=1 for a=b and 0 otherwise (this is the normalization integral of radial functions).

Second equation

In the same way the equation for Y can be reduced to (here f is the rhs of the equation for Y):

(Subscript[y, j-2]) = (Subscript[y, j]) +  (h/90)(-(Subscript[f, j+1]) + 34 (Subscript[f, j]) +114 (Subscript[f, j-1]) + 34 (Subscript[f, j-2]) - (Subscript[f, j-3])).

Exponential grid

Wave functions in a multi-electrons atom can have quite different scales, and in any case the characteristic lenghts vary with Z, the atomic number. It is convenient in the numerical solution of the problem to have a grid which an reasonably works for every Z and for the different electrons in the same atom.

We will use an exponential grid defined by the change of variables

r = x/Z;
x = (Subscript[x, 0])*(Exp[ξ]-1);   we will put   ρ≡ (Subscript[x, 0])Exp[ξ]= x+(Subscript[x, 0]).

In these variables the scalar product of the reduced radial wave functions becomes

P[a,x] = (Sqrt[ρ])ψ[a,ξ];   (Subscript[δ, ab])=∫P[a,x]P[b,x]dx =∫(ρ^2)ψ[a,x]ψ[b,x]dξ

while equations () take the form

(d/(dξ))(x^k) (Z^k)[a,b,x]=(ρ^2)(r^k) ψ[a,x]ψ[b,x];
(d/(dξ))(((Y^k)[a,b,x])/(x^(k+1)))= - ((2k+1)/(r^(k+2)))ρ (Z^k)[a,b,x];

We see that these equations are identical in form to equations (), only the rhs changes, then the recurrence relations () and () are valid.

Let us consider now the behavior of (Z^k)as r→0.  The change of variables  r = x/Z produce an overall factor Z which is reabsorbed by the change in the normalizationof the wave function P = (Sqrt[Z])P[x]. The term Zr in bracket becomes x, then:

(Z^k)[a,b,x] ∼ c1 x Pa Pb(1 + c2 x)

Passing to ξ variable and using the notation fH = (ρ^2) ψ[a,x]ψ[b,x] for the term independent of k we have

(Z^k) ≃ c1 x ρ (Subscript[ψ, a]) (Subscript[ψ, b])(1 + c2 x) = c1 x fH/ρ (1 + c2 x) = c1 x fH/(x+x0) (1 + c2 x)

Integration of the potentials

ClearAll["Global`*"];

Numerical functions for Y and Z

The procedure described in the first section is translated in Mathematica language in two functions funZ and funY. The caling sequence for these functions is

funZ[k,La,Lb,fa,fb,x,x0,h];         funY[k,La,Lb,fa,fb,fZ,x,x0,h];

where

(* Input functions are rescaled rad. wave funct. ψ = P/(Sqrt[ρ]) *)
funZ=
Compile[{{valk,_Integer},{La,_Integer},{Lb,_Integer},
{fa,_Real,1},{fb,_Real,1},{x,_Real,1},{x0,_Real},{h,_Real}},
Module[{nL,fH,fun,cost1,cost2,fZ,j,k,kind},
(nL = Length[fa];
If[valk== 0,fH = (x+x0)^2*fa*fb,fH = x^valk*(x+x0)^2*fa*fb];
If[La == Lb && fa[[3]] == fb[[3]] && fa[[5]] == fb[[5]] && valk==0,
asint =1.;,asint = 0.];
fZ=Table[0.,{j,1,nL}];
cost1= N[(1./(La+Lb+valk+3.))];
cost2=N[(1./(La+Lb+valk+4.))((1/(La+1.))+(1/(Lb+1.)))];
(* qui va r* Sqrt[rho] ψ_a Sqrt[rho] ψ_b = r fun /(r+r0) *)
fZ[[1]] = cost1*x[[1]]*fH[[1]]/(x[[1]]+x0)*(1+cost2*x[[1]]);
fZ[[2]] = cost1*x[[2]]*fH[[2]]/(x[[2]]+x0)*(1+cost2*x[[2]]);
fZ[[3]] = fZ[[1]]+ h/3.*(fH[[1]]+4.*fH[[2]]+fH[[3]]);
Table[fZ[[j]] = fZ[[j-2]]+ h/90.*(-fH[[j+1]] + 34.fH[[j]] +114.*fH[[j-1]] +
34.*fH[[j-2]] - fH[[j-3]]);,{j,5,nL-1,2}];
Table[fZ[[j]] = fZ[[j-2]]+ h/90.*(-fH[[j+1]] + 34.fH[[j]] +114.*fH[[j-1]] +
34.*fH[[j-2]] - fH[[j-3]]);,{j,4,nL-2,2}];
fZ[[nL]] = fZ[[nL-2]]+ h/3.*(fH[[nL-2]]+4.*fH[[nL-1]]+fH[[nL]]);
If[asint > 0.1,
fZ[[Range[1,nL-1,2]]] = fZ[[Range[1,nL-1,2]]] + asint - fZ[[nL-1]];
fZ[[Range[2,nL,2]]] = fZ[[Range[2,nL,2]]] + asint - fZ[[nL]]; ,];
If[valk== 0,fZ = fZ ;,fZ = fZ/x^valk];
fZ)]];

funY=
Compile[{{valk,_Integer},{La,_Integer},{Lb,_Integer},
{fa,_Real,1},{fb,_Real,1},{fZ,_Real,1},{x,_Real,1},{x0,_Real},{h,_Real}},
Module[{nL,fH,j,corr,k,fY},
nL = Length[fa];
fH = (2.*valk+1.)*(x+x0)*fZ/x^(valk+2);
fY=Table[0.,{j,1,nL}];
fY[[nL]] = fZ[[nL]]/x[[nL]]^(valk+1);
fY[[nL-1]]=  fZ[[nL-1]]/x[[nL-1]]^(valk+1);
Table[fY[[j-2]] = fY[[j]]+ 
h/90.*(-fH[[j+1]] + 34.fH[[j]] +114.*fH[[j-1]] +
34.*fH[[j-2]] - fH[[j-3]]);,{j,nL-1,5,-2}];
fY[[1]] = fY[[3]] + h/3.*(fH[[1]]+4.*fH[[2]]+fH[[3]]);
fY[[nL-2]] = fY[[nL]] +h/3.*(fH[[nL-2]]+4.*fH[[nL-1]]+fH[[nL]]);
Table[fY[[j-2]] = fY[[j]]+
 h/90.*(-fH[[j+1]] + 34.fH[[j]] +114.*fH[[j-1]] +
34.*fH[[j-2]] - fH[[j-3]]);,{j,nL-2,4,-2}];
fY = fY*x^(valk+1);
fY]];

An utility function which can be useful is a function which performs the integral of a function defined on the exponential grid, i.e. in the ξ variable. As input we have to provide the grid values and the parameter h.

integ = Compile[{{f,_Real,1},{h,_Real}},
Module[{nL,f0,val},
nL = Length[f];
f0=0.;val=0.;
val = h*(119./90.*f[[1]] + 32./45.f[[2]] + 113./90.f[[3]] +
113./90. f[[nL-1]] + 17./45.f[[nL]]);
val = val + h/90.*68.*Total[f[[Range[4,nL-2,2]]]] +
 h/90.*112.*Total[f[[Range[5,nL-3,2]]] ];
val]];

The grid

Here is an axample of a function which produces an exponential grid. We use the same parameters that will be used in the solution of Hartree-Fock equations. The parameter Z appearing in the function will be the atomic number in the Hartree-Fock equation, while numrad will play the role of radial quantum number.

faiGriglia[h_,x0_]:=
Module[{nP,Z,numrad,x,rho},
Z = 1; numrad = Range[4];
xmax = 20.*Max[numrad]^2 * Z;
nP = Ceiling[ Log[xmax/x0]/h];
nP = nP + Mod[nP,2];
x = x0*(Exp[h*Range[nP]] - 1.);
rho = x+x0; 
{x,rho}];

Numerical checks

In this section we present some example of calculations for funZ and funY. We will use Hydrogen functions and we will compare the results with analytical values.

Hydrogen functions

These are radial reduced wave functions for Hydrogen

χ[r_, n_, L_] := r* (2/(n^2))(Sqrt[(((n-L-1)!)/((n+L)!))])*
(((((2 r)/n)))^L) * (E^(-(r/n)))* LaguerreL[n-L-1, 2 L + 1,((2 r)/n)] ;

The analitic expression for functions Y and Z is

y[k_,n1_,L1_,n2_,L2_,r_]:=
Integrate[((((s/r)))^k)χ[s, n1, L1] χ[s, n2, L2] ,{s,0,r},GenerateConditions->False] +
Integrate[((((r/s)))^(k+1))χ[s, n1, L1] χ[s, n2, L2] ,{s,r,Infinity},GenerateConditions->False]

z[k_,n1_,L1_,n2_,L2_,r_]:=
Integrate[((((s/r)))^k)χ[s, n1, L1] χ[s, n2, L2] ,{s,0,r},GenerateConditions->False]

Tests

Test for Zk functions

This define the grid

h = 1./10; x0 = 0.1;
{x,rho} = faiGriglia[h,x0];

Here we write the reduced wave function, divided by (Sqrt[ρ]). The function chosen is the 1s radial function.

psi =Chop[ N[χ[x, 1, 0]]/Sqrt[rho]];

The analytical value computed in grid points

testZ1s = z[0,1,0,1,0,zz]/.zz-> x;

The numerical value computed with funZ

calcZ1s=funZ[0,0,0,psi,psi,x,x0,h];

Max[Abs[calcZ1s - testZ1s]]

Test for Yk functions

The same kind of test can be done for funY. Using the previous computed value for Z0

testY1s = y[0,1,0,1,0,zz]/.zz-> x;
calcY1s=funY[0,0,0,psi,psi,calcZ1s,x,x0,h];
Max[Abs[calcY1s - testY1s]]

Test for several functions with different k values

Here we organize a more systematic test

testZ = Range[100];
testY = Range[100];
numeri = Range[100];
kind=0;
Do[
kind=kind+1;
numeri[[kind]] = {k,na,la,nb,lb};
testZ[[kind]]=z[k,na,la,nb,lb,zz]/.zz-> x;
testY[[kind]]=y[k,na,la,nb,lb,zz]/.zz-> x;,
{na,1,2},{nb,1,na},{la,0,na-1},{lb,0,nb-1},{k,Abs[lb-la],la+lb,2}];
testZ=Take[testZ,kind];
testY=Take[testY,kind];
numeri=Take[numeri,kind];

calcZ = Range[100];
calcY = Range[100];
kind=0;
Do[
kind=kind+1;
Qa =Chop[ N[χ[x, na, la]]/Sqrt[rho]];
Qb =Chop[ N[χ[x, nb, lb]]/Sqrt[rho]];
calcZ[[kind]]=funZ[k,la,lb,Qa,Qb,x,x0,h];
calcY[[kind]]=funY[k,la,lb,Qa,Qb,calcZ[[kind]],x,x0,h];,
{na,1,2},{nb,1,na},{la,0,na-1},{lb,0,nb-1},{k,Abs[lb-la],la+lb,2}];
calcZ=Take[calcZ,kind];
calcY=Take[calcY,kind];

diffZ=Table[Max[Abs[calcZ[[kind]]-testZ[[kind]]]],{kind,1,Length[testZ]}];
diffY=Table[Max[Abs[calcY[[kind]]-testY[[kind]]]],{kind,1,Length[testY]}];

TableForm[Table[Join[numeri[[i]],{diffZ[[i]],diffY[[i]]}],{i,1,Length[numeri]}],
TableHeadings->{{},{k,"na","la","nb","lb","diffZ","diffY"}},TableSpacing->{1,4}]

Test on a finer grid

We change the values of the parameter h

h = 1./20; x0 = 0.1;
{x,rho} = faiGriglia[h,x0];

and repeat the test

testZ = Range[100];
testY = Range[100];
numeri = Range[100];
kind=0;
Do[
kind=kind+1;
numeri[[kind]] = {k,na,la,nb,lb};
testZ[[kind]]=z[k,na,la,nb,lb,zz]/.zz-> x;
testY[[kind]]=y[k,na,la,nb,lb,zz]/.zz-> x;,
{na,1,2},{nb,1,na},{la,0,na-1},{lb,0,nb-1},{k,Abs[lb-la],la+lb,2}];
testZ=Take[testZ,kind];
testY=Take[testY,kind];
numeri=Take[numeri,kind];

calcZ = Range[100];
calcY = Range[100];
kind=0;
Do[
kind=kind+1;
Qa =Chop[ N[χ[x, na, la]]/Sqrt[rho]];
Qb =Chop[ N[χ[x, nb, lb]]/Sqrt[rho]];
calcZ[[kind]]=funZ[k,la,lb,Qa,Qb,x,x0,h];
calcY[[kind]]=funY[k,la,lb,Qa,Qb,calcZ[[kind]],x,x0,h];,
{na,1,2},{nb,1,na},{la,0,na-1},{lb,0,nb-1},{k,Abs[lb-la],la+lb,2}];
calcZ=Take[calcZ,kind];
calcY=Take[calcY,kind];

diffZ=Table[Max[Abs[calcZ[[kind]]-testZ[[kind]]]],{kind,1,Length[testZ]}];
diffY=Table[Max[Abs[calcY[[kind]]-testY[[kind]]]],{kind,1,Length[testY]}];

TableForm[Table[Join[numeri[[i]],{diffZ[[i]],diffY[[i]]}],{i,1,Length[numeri]}],
TableHeadings->{{},{k,"na","la","nb","lb","diffZ","diffY"}},TableSpacing->{1,4}]

Test with h=1/100

With small h the agreement is excellent

h = 1./100; x0 = 0.05;
{x,rho} = faiGriglia[h,x0];

testZ = Range[100];
testY = Range[100];
numeri = Range[100];
kind=0;
Do[
kind=kind+1;
numeri[[kind]] = {k,na,la,nb,lb};
testZ[[kind]]=z[k,na,la,nb,lb,zz]/.zz-> x;
testY[[kind]]=y[k,na,la,nb,lb,zz]/.zz-> x;,
{na,1,2},{nb,1,na},{la,0,na-1},{lb,0,nb-1},{k,Abs[lb-la],la+lb,2}];
testZ=Take[testZ,kind];
testY=Take[testY,kind];
numeri=Take[numeri,kind];

calcZ = Range[100];
calcY = Range[100];
kind=0;
Do[
kind=kind+1;
Qa =Chop[ N[χ[x, na, la]]/Sqrt[rho]];
Qb =Chop[ N[χ[x, nb, lb]]/Sqrt[rho]];
calcZ[[kind]]=funZ[k,la,lb,Qa,Qb,x,x0,h];
calcY[[kind]]=funY[k,la,lb,Qa,Qb,calcZ[[kind]],x,x0,h];,
{na,1,2},{nb,1,na},{la,0,na-1},{lb,0,nb-1},{k,Abs[lb-la],la+lb,2}];
calcZ=Take[calcZ,kind];
calcY=Take[calcY,kind];

diffZ=Table[Max[Abs[calcZ[[kind]]-testZ[[kind]]]],{kind,1,Length[testZ]}];
diffY=Table[Max[Abs[calcY[[kind]]-testY[[kind]]]],{kind,1,Length[testY]}];

TableForm[Table[Join[numeri[[i]],{diffZ[[i]],diffY[[i]]}],{i,1,Length[numeri]}],
TableHeadings->{{},{k,"na","la","nb","lb","diffZ","diffY"}},TableSpacing->{1,4}]