4
$\begingroup$

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?

$\endgroup$
11
  • 1
    $\begingroup$ 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[π]) $\endgroup$
    – Roman
    Jun 5 at 19:52
  • $\begingroup$ 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] $\endgroup$
    – Michael E2
    Jun 5 at 19:52
  • 1
    $\begingroup$ 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! $\endgroup$
    – Roman
    Jun 5 at 20:00
  • 2
    $\begingroup$ @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 $\endgroup$
    – Riccardo
    Jun 5 at 20:05
  • 2
    $\begingroup$ @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. $\endgroup$
    – Roman
    Jun 6 at 8:58

1 Answer 1

11
$\begingroup$

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    *)
$\endgroup$
2
  • $\begingroup$ The code shared by Roman produces the correct solution for this integral $\endgroup$
    – Riccardo
    Jun 11 at 18:41
  • 1
    $\begingroup$ @Riccardo, if you feel that Roman's answer has been satisfactory, do click on the check mark at the left of his answer. $\endgroup$ Jun 11 at 18:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.