# NIntegrate fails to converge to desired accuracy

I am calculating an integral with Mathematica. This is the code I am using:

d = 1.401; RA = {0, 0, 0}; RB = {d, 0, 0};
a1 = 0.416613; a2 = 0.664404447997;

g1[x_, y_, z_] := (2 a1/Pi)^(3/4) Exp[-a1 ({x, y, z} - RA).({x, y, z} - RA)];
g2[x_, y_, z_] := (2 a1/Pi)^(3/4) Exp[-a1 ({x, y, z} - RB).({x, y, z} - RB)];

chi1[x_, y_, z_] := (g1[x, y, z] + g2[x, y, z])/(Sqrt[2] Sqrt[1 + a2])
chi2[x_, y_, z_] := (g1[x, y, z] - g2[x, y, z])/(Sqrt[2] Sqrt[1 - a2])


This is the integral I am trying to calculate

NIntegrate[(
chi1[x1, y1, z1])^2*(chi1[x2, y2, z2])^2*1/Sqrt[(x1 - x2)^2 + (y1 - y2)^2 + (z1 -z2)^2],
{x1, -Infinity, Infinity}, {y1, -Infinity, Infinity}, {z1, -Infinity, Infinity},
{x2, -Infinity, Infinity}, {y2, -Infinity, Infinity}, {z2, -Infinity, Infinity}]


Unforunately, the integral does not converge. I know that the result should be 0.677391. Can someone help me out? I tried using MonteCarlo methods, but in principle I would like to obtain the exact solution. Maybe I should add an Exclusion region?

• This integral can be calculated analytically. Are you particularly interested in the Monte Carlo method or are you looking for the exact result? The exact result is (Sqrt[a1]*(d+2d*E^(-a1*d^2)+(Sqrt[π](8E^(-a1*d^2/2)*Erf[Sqrt[a1]*d/2]+Erf[Sqrt[a1]*d]))/(2*Sqrt[a1])))/((1+a2)^2*d*Sqrt[π]) Jun 5 at 19:52
• Might be overkill but how's this: NIntegrate[(chi1[x1, y1, z1])^2*(chi1[x2, y2, z2])^2*1/ Sqrt[(x1 - x2)^2 + (y1 - y2)^2 + (z1 - z2)^2], {x1, -Infinity, 0, d, 2 d, Infinity}, {y1, -Infinity, -1, 0, 1, Infinity}, {z1, -Infinity, -1, 0, 1, Infinity}, {x2, -Infinity, 0, d, 2 d, Infinity}, {y2, -Infinity, -1, 0, 1, nfinity}, {z2, -Infinity, -1, 0, 1, Infinity}, Method -> {"MonteCarlo"}, MaxPoints -> 1000000] Jun 5 at 19:52
• Substitute $x_2=x_1+r \sin\theta\cos\phi$, $y_2=y_1+r \sin\theta\sin\phi$, $z_2=z_1+r \cos\theta$ and simplify the integrand, then integrate one coordinate after another. Notice that the denominator only depends on $r$ and so all other integrations are fairly simple. Don't forget the Jacobian of $r^2\sin\theta$ that comes from the transformation to spherical coordinates! Jun 5 at 20:00
• @Roman If it is not too much to ask, do you think you could share the code? By the way, your manual on using Mathematica for quantum mechanics is great! I used it in the past Jun 5 at 20:05
• @Riccardo Please have a look at this answer. If you can help me figure out my last roadblock for using pygments/minted, then I'd recommend that. Jun 6 at 8:58

My usual answer for the numerical calculation of high-dimensional integrals is: try to do as many dimensions as possible analytically, and then use numerical integration for the remaining dimensions. This is true even if some of the analytical integrations result in special function (hypergeometric, Bessel, etc.) because the algorithms for evaluating these special functions are much more efficient than numerical integration.

In the present case, a judicious coordinate transformation allows us to do all of the integrations analytically.

RA = {0, 0, 0};
RB = {d, 0, 0};
g1[x_, y_, _] = (2 a1/Pi)^(3/4) Exp[-a1 ({x, y, z} - RA) . ({x, y, z} - RA)];
g2[x_, y_, z_] = (2 a1/Pi)^(3/4) Exp[-a1 ({x, y, z} - RB) . ({x, y, z} - RB)];
chi1[x_, y_, z_] = (g1[x, y, z] + g2[x, y, z])/(Sqrt[2] Sqrt[1 + a2]);
chi2[x_, y_, z_] = (g1[x, y, z] - g2[x, y, z])/(Sqrt[2] Sqrt[1 - a2]);


The integrand:

J = (chi1[x1, y1, z1])^2*(chi1[x2, y2, z2])^2*1/Sqrt[(x1 - x2)^2 + (y1 - y2)^2 + (z1 - z2)^2] // FullSimplify
(*    (2 a1^3 E^(-2 a1 (2 d^2 + x1^2 + x2^2 + y1^2 + y2^2 + z1^2 +
z2^2)) (E^(a1 d^2) + E^(2 a1 d x1))^2 (E^(a1 d^2) +
E^(2 a1 d x2))^2)/((1 + a2)^2 π^3 Sqrt[(x1 - x2)^2 +
(y1 - y2)^2 + (z1 - z2)^2])                                     *)


Substitute relative coordinates: $$x_2=x_1+dx$$, $$y_2=y_1+dy$$, $$z_2=z_1+dz$$, and notice that the denominator is now independent of $$x_1,y_1,z_1$$, which makes these integrals simple (analytic):

Assuming[a1 > 0, Integrate[J /. {x2 -> x1 + dx, y2 -> y1 + dy, z2 -> z1 + dz},
{x1, -∞, ∞}, {y1, -∞, ∞}, {z1, -∞, ∞}]]

(*    (a1^(3/2) E^(-a1 ((d + dx)^2 + dy^2 + dz^2)) (1 + 4 E^(2 a1 d dx) +
E^(4 a1 d dx) + 2 E^(a1 d (d + 2 dx)) + 4 E^(1/4 a1 d (d + 4 dx)) +
4 E^(1/4 a1 d (d + 12 dx))))/(4 (1 + a2)^2 Sqrt[dx^2 + dy^2 + dz^2] π^(3/2))    *)


Go to spherical coordinates: express the vector $$(dx,dy,dz)$$ in spherical coordinates $$(r,\theta,\phi)$$,

Assuming[r > 0, % /. {dx -> r Sin[θ] Cos[φ],
dy -> r Sin[θ] Sin[φ],
dz -> r Cos[θ]} // FullSimplify]
(*    (a1^(3/2) E^(-a1 (d^2 + r^2 + 2 d r Cos[φ] Sin[θ])) (1 +
4 E^(2 a1 d r Cos[φ] Sin[θ]) + E^(
4 a1 d r Cos[φ] Sin[θ]) +
2 E^(a1 d (d + 2 r Cos[φ] Sin[θ])) +
4 E^(1/4 a1 d (d + 4 r Cos[φ] Sin[θ])) +
4 E^(1/4 a1 d (d + 12 r Cos[φ] Sin[θ]))))/(4 (1 +
a2)^2 π^(3/2) r)                                      *)


spherical integration (don't forget the Jacobian $$r^2\sin\theta$$):

Assuming[a1 > 0 && d > 0, Integrate[% r^2 Sin[θ], {r, 0, ∞}, {θ, 0, π}, {φ, 0, 2 π}]]
(*    (Sqrt[a1] (d + 2 d E^(-a1 d^2) + (Sqrt[π] (8 E^(-((a1 d^2)/2)) Erf[(Sqrt[a1] d)/2] +
Erf[Sqrt[a1] d]))/(2 Sqrt[a1])))/((1 + a2)^2 d Sqrt[π])    *)


Check numerical value:

% /. {d -> 1.401, a1 -> 0.416613, a2 -> 0.664404447997}
(*    0.677391    *)

• The code shared by Roman produces the correct solution for this integral Jun 11 at 18:41
• @Riccardo, if you feel that Roman's answer has been satisfactory, do click on the check mark at the left of his answer. Jun 11 at 18:55