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[π]) Commented Jun 5, 2022 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] Commented Jun 5, 2022 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! Commented Jun 5, 2022 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 Commented Jun 5, 2022 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. Commented Jun 6, 2022 at 8:58

1 Answer

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.

Start with analytic definitions (leaving all numerical constants undefined for now):

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 Commented Jun 11, 2022 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. Commented Jun 11, 2022 at 18:55