2
$\begingroup$

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...

$\endgroup$
1
  • 3
    $\begingroup$ 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}]. $\endgroup$
    – Michael E2
    Commented Oct 13, 2015 at 2:07

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.