phi[a_,r_]:=(2a/π)^(3/4)Exp[-a r^2]
P1s[Z_,r_]:=Z^(3/2)/Sqrt[π]Exp[-Z r]
P2s[Z_,r_]:=Z^(3/2)/(4Sqrt[2π])(2-Z r)Exp[-Z r/2]
sij[a1_,a2_]:=(2Sqrt[a1 a2]/(a1+a2))^(3/2)
hij[a1_,a2_]:=(3a1 a2/(a1+a2)-Z 2/Sqrt[π]Sqrt[a1+a2])sij[a1,a2]
gijkl[a1_,a2_,a3_,a4_]:=16/Sqrt[π]Sqrt[a1 a2 a3 a4]^(3/2)/((a1+a2)(a3+a4)Sqrt[a1+a2+a3+a4])