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[[1]] - mofunc2 gradmo1[[1]];
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[[1]]
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...
Compile
, you wantEvaluate[jX /. {x -> x1, y -> y1, z -> z1}]
. $\endgroup$