integration of a list of gaussian type orbitals (GTO)

I am currently trying to integrate various components of a current density (j={jx,jy,jz}) over space. The problem is that the current j is extracted from a quantum chemical computation and is a sum of many simple functions (gaussian functions). Let us focus on jX :

jX = mofunc1 gradmo2[] - mofunc2 gradmo1[];


Where mofunc1 (mofunc2) is a weighted sum of 239 atomic orbitals (that are a sum of up to ten gaussians coef Exp[-expo (x^2+y^2+z^2)] with various coefficients and centered at different points in space). And gradmo1[] is the x derivative along x of mofunc1.

Now, the integration : I tried to compile the function and integrate it in the following way :

jCX = Compile[{{x1, _Real}, {y1, _Real}, {z1, _Real}},
Evaluate@jX /. {x -> x1, y -> y1, z -> z1}, CompilationTarget -> "C"];

jXarg[x1_?NumericQ, y1_?NumericQ, z1_?NumericQ] := jCX[x1, y1, z1]

NIntegrate[jXarg[x1, y1, z1] , {x1, -bd, +bd}, {y1, -bd, +bd}, {z1, -bd, +bd},
Method -> "QuasiMonteCarlo"]


I set the boundary bd (=5) to a value where I am quite sure that jx=0 with a good accuracy. I am stuck at this point as the code will never end or not converge if I set a number for MaxPoints. Other integration methods were not succesful too so far. I also tried unsuccesfully to Expand and/or Simplify before giving the function to NIntegrate but it also take forever.

Some informations about evaluation timings for j on a specific value :

jX /. {x -> 1., y -> 1., z -> 1.} // Timing
jCX[1., 1., 1.] // Timing
jXarg[1., 1., 1.] // Timing

{0.014889, 0.000161008}
{0.02063, 0.000161008}
{0.014801, 0.000161008}


So it looks it make things only marginally better to go through those compilation. Or maybe I am not doing it right... Please let me know if there is some optimization that I could do to make that integral converge.

Thank you Michael E2 for your reply. In fact, I changed my strategy and decided to evaluate all the integrals symbolically. For that, I modified my parser of the molden file that I get out the quantum chem program in such a way to have all the molecular orbitals (MO) in the following structure : $\{\{\text{coef},k,l,m,{x_p,y_p,z_p},\alpha\},...\}$ Where each element represent a gaussian orbital contribution to the MO of the form : $$\text{coef} \ (x-x_p)^k (y-y_p)^l (z-z_p)^m e^{- \alpha ((x-x_p)^2+(y-y_p)^2+(z-z_p)^2)}$$ Then, I made a routing to do the integral symbolically by running over the list of gaussian orbtials.

  integrate[listmod_] := Block[{list = listmod, k, l, m, i, alpha},
int = 0;
For[i = 1, i <= Length[list], i++,
k = list[[i, 2]];
l = list[[i, 3]];
m = list[[i, 4]];
alpha = list[[i, 6]];
val = list[[i, 1]] 1/
8 (1 + (-1)^k) (1 + (-1)^l) (1 + (-1)^m) alpha^(
1/2 (-3 - k - l - m))
Gamma[(1 + k)/2] Gamma[(1 + l)/2] Gamma[(1 + m)/2];
int += val;
];
Return[int];
]


I also made routines that makes a gradient of a list of gaussian orbitals and store it as a similar list and the same goes for a product (so I can calculate the current density operator). It is probably not well optimized but it does the job. The codes are a bit lengthy so I do not know if it is the place to share more of it here...

• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – Michael E2 Oct 13 '15 at 2:05
• Probably the complete code will be needed, or a minimal example with the same problem. One thing I can spot is that in Compile, you want Evaluate[jX /. {x -> x1, y -> y1, z -> z1}]. – Michael E2 Oct 13 '15 at 2:07