I want to build a fast way to convert a density grid file into a grid file of the respective electrostatic potential (ESP; also known as MEP, the molecular electrostatic potential).

This calculation is based on the formula: $$ V_\mathrm{tot}(\mathbf{r}) = V_\mathrm{nuc}(\mathbf{r}) + V_{ele}(\mathbf{r}) = \sum_A \frac{Z_A}{|\mathbf{r}-\mathbf{R}_A|} - \int \frac{\rho(\mathbf{r}^\prime)}{|\mathbf{r}-\mathbf{r}^\prime|} \mathrm d \mathbf{r}^\prime $$

with $Z_A$ as atom A's charge, $\mathbf{R}_A$ as it's point in space and $\rho(\mathbf{r})$ as the density at any point $\mathbf{r}$.

As I work with a grid and not with a function for $\rho$, the equation turns into this: $$ V_\mathrm{tot}(\mathbf{r}_i) = \sum_A \frac{Z_A}{|\mathbf{r}_i-\mathbf{R}_A|} - \Delta x \Delta y \Delta z \sum_j \frac{\rho(\mathbf{r}_j)}{|\mathbf{r}_i-\mathbf{r}_j|} $$

where the integral was converted into a sum that is multiplied with the distance between two points in x-, y, or z-direction. (How cube files are formatted can be found here (by Paul Bourke) and here (from Gaussian).)

The best version that I could produce so far looks like this:

EucDist = 
     (Sqrt[Total[(x - #)^2]] + 10.^-7) & /@ y, 
   CompilationTarget -> "C", RuntimeAttributes -> {Listable}];
vnuc[charges_, molecule_, xyz_] := 
  With[{distances = EucDist[xyz, molecule]}, Total[charges/distances]];
ESP[cub_] := Block[
  (* Definition of the variables *)
  {nAtoms, xStart, yStart, zStart, xStep, dx, yStep, dy, zStep, dz, 
   molecule, values, data, f, xVals, yVals, zVals, xyzVals},
  (* Input of the fube-file's parameters *)
  {nAtoms, xStart, yStart, zStart} = cub[[3]];
  {xStep, dx} = cub[[4, 1 ;; 2]];
  {yStep, dy} = {cub[[5, 1]], cub[[5, 3]]};
  {zStep, dz} = {cub[[6, 1]], cub[[6, 4]]};
  molecule = cub[[7 ;; nAtoms + 6]];
  (* Definition of the grid points *)
  xyzVals = 
   Flatten[Table[{xStart + (i - 1) dx, yStart + (j - 1) dy, 
      zStart + (k - 1) dz}, {i, 1, xStep}, {j, 1, yStep}, {k, 1, 
      zStep}], 2];
  (* Input of the grid and conversion into a list *)
  values = Flatten[cub[[nAtoms + 7 ;;]]];
  (* Definition of the V_tot function as: f(xyz)=Subscript[V, nuc](charges, 
  atom positions, grid point)-Subscript[V, elec](density values, 
  complete grid, grid point) *)
  f[xyz_] := 
   vnuc[molecule[[All, 2]], molecule[[All, 3 ;; 5]], xyz] - 
    vnuc[values, xyzVals, xyz] dx dy dz;
  (* mapping f on the grid points *)
  ParallelMap[f, xyzVals][[1]]

But at the current state it takes way too long to be evaluated. As I'd like to use it on much bigger grids, this function needs to be more efficient.

Now my question is: How can I come up with a function that is a lot faster than my current version?

Test files:

P.s. Up to now I did not care about exporting the output as the calculation takes too much time to be useful. But nonetheless I already have a function to turn Mathematica-style numbers into the fortran-style numbers.

g[n_] := EngineeringForm[N[n], NumberFormat -> (Row[{#1, "E", If[#3 == "", "0", #3]}] &)]
  • 1
    $\begingroup$ Re: CompiledFunction::cfta: error: The variable specification should be {{x, _Real, 1}, {y, _Real, 1}} as {-6.,-6.,-6.703122} is a list of depth 1 and the third argument is the depth, not the length of the argument. $\endgroup$ – shrx Apr 18 '16 at 10:05
  • $\begingroup$ @shrx thank you! Everyone else - don't hesitate to tell me if something is unclear. $\endgroup$ – pH13 - Yet another Philipp Apr 19 '16 at 7:55
  • 1
    $\begingroup$ @pH13-YetanotherPhilipp, why does your vnuc function return a list of 3 numbers and not just one number? I would expect the MEP at a particular point to be a single number, but if I evaluate vnuc[values, xyzVals, {-6., -6., -6.}] dx dy dz then I get {15.0976, 15.0976, 0.199034} in return $\endgroup$ – Jason B. Apr 25 '16 at 10:30

Your cube file had a very large grid ( 117*117*130 = 1779570), and 2 million points is just far too many for testing a function.

So I created cube files for the electron density and electrostatic potential for the molecule furan, using a much sparser grid (around 8000 grid points instead). Here they are:

Density cube file

Potential cube file

Now that we have them downloaded, let's try your code and see how long this takes,

start = AbsoluteTime[];
EucDist = 
  Compile[{{x, _Real, 1}, {y, _Real, 1}}, Sqrt[Total[(x - y)^2]]];
vnuc[charges_, molecule_, xyz_] := 
  With[{distances = Map[EucDist[xyz, #] + 10^-7 &, molecule]}, 
cub = Import["~/Downloads/furandensity_coarse.cube", "Table"];
xyzVals = 
     zStep}], 2];
f[xyz_] := 
   vnuc[values, xyzVals, xyz] dx dy dz;
pot = ParallelMap[f, xyzVals];

AbsoluteTime[] - start


So 24 seconds from start to finish. Now let's try my method, which uses the functions included in this file, and are detailed here and here.

start = AbsoluteTime[];
<< "ReadWriteCube.m"
{cubeDat, xgrid, ygrid, zgrid, xyzText, headerTxt} = 
angstromToBohr = 1.889725989;
atomlist = First /@ ImportString[xyzText, "Table"][[3 ;;]];
atomcoords = Rest /@ ImportString[xyzText, "Table"][[3 ;;]];
atomcoords = atomcoords angstromToBohr;
densitydata = Transpose[cubeDat, {3, 2, 1}];
rgrid = Table[{x, y, z}, {z, zgrid}, {y, ygrid}, {x, xgrid}];
{dx, dy, dz} = First@*Differences /@ {xgrid, ygrid, zgrid};
{xl, yl, zl} = Length /@ {xgrid, ygrid, zgrid};
{xgrid2, ygrid2, zgrid2} = 
  Table[x, {x, -(#1 - 1) #2, (#1 - 1) #2, #2}] & @@@ {{xl, dx}, {yl, 
     dy}, {zl, dz}};
rgrid2c = 
  1/Table[Norm[{x, y, z}] + 10^-7, {z, zgrid2}, {y, ygrid2}, {x, 
nuclearPotential = 
  Sum[ElementData[atomlist[[n]], "AtomicNumber"] Map[
     1/(Norm[# - atomcoords[[n]]] + 10^-7) &, rgrid, {3}], {n, 
epfunc[i_, j_, k_] := 
  dx dy dz Flatten[
     rgrid2c[[zl - i + 1 ;; 2 zl - i, yl - j + 1 ;; 2 yl - j, 
       xl - k + 1 ;; 2 xl - k]]].Flatten[densitydata];
electricPotential = 
   Array[epfunc, Dimensions@densitydata]; // AbsoluteTiming
pot2 = nuclearPotential - electricPotential;
AbsoluteTime[] - start;


There is a similar speedup if parallel versions of the Array, Map, etc. are used. The speedup is even greater if the grid is larger.

The reason for the speedup is that I only calculate the distances once. The various r values are on a grid, so their distances follow a pattern. You can verify for yourself that the potentials calculated in this manner are identical. The potentials are in a different format, since yours is a 1-dimensional list and mine is a 3 dimensional grid,

Flatten[Transpose[pot2, {3, 2, 1}]] - pot // MinMax
(* {-3.65544*10^-8, 2.23517*10^-8} *)

So they are the same with a small numerical difference probably related to the compiled function.

So this speeds you up somewhat, but if you want to do this on a grid with 10 million points, well, it doesn't speed it up enough.

But this is all superfluous anyway, since the resulting potential is not in any way similar to the potential given by the electronic structure code. Here is the potential generated by Gaussian, plotted with isosurfaces at $\pm 0.02$ (and here's the real reason to answer this question, to show off the CubePlot function),


enter image description here

See the negative region around the oxygen and above the center of the ring? Let's compare this with the potential derived directly from the density (transposing the potential to put it back in the form from the Gaussian input):

CubePlot[{Transpose[pot2,{3,2,1}], xgrid, ygrid, zgrid, xyzText}]

Mathematica graphics

Even using a larger value for the isosurface value doesn't help.

CubePlot[{pot2, xgrid, ygrid, zgrid, xyzText}, Contours -> {-.5, .5}]

Mathematica graphics

So in conclusion, you can speed up your function by pre-computing the finite number of inter-grid distances, but the basic formula does not appear to be sound.

  • $\begingroup$ Is .cube import (with the result as maybe an InterpolatingFunction) coming to Mathematica as a built-in one of these days? Of course I don't think formats can be pushed as paclets which makes it somewhat less useful for me, but there are theoretically people for whom it would be useful. $\endgroup$ – b3m2a1 Aug 6 '17 at 6:37

Your Answer

