Calculate the ESP Grid based on Density Grid

Question

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 = 
  Compile[{{x, _Real, 1}, {y, _Real, 1}}, 
     (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:

• Initial density file (H$_2$): http://www.file-upload.net/download-11494733/density.cub.html
• Resulting ESP file: http://www.file-upload.net/download-11494735/totesp.cub.html

---

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]}] &)]

Answer 1

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]}, 
   Total[charges/distances]];
cub = Import["~/Downloads/furandensity_coarse.cube", "Table"];
{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]];
values = Flatten[cub[[nAtoms + 7 ;;]]];
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];
f[xyz_] := 
  vnuc[molecule[[All, 2]], molecule[[All, 3 ;; 5]], xyz] - 
   vnuc[values, xyzVals, xyz] dx dy dz;
pot = ParallelMap[f, xyzVals];

AbsoluteTime[] - start

24.107253

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} = 
  ReadCube["~/Downloads/furandensity_coarse.cube"];
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, 
     xgrid2}];
nuclearPotential = 
  Sum[ElementData[atomlist[[n]], "AtomicNumber"] Map[
     1/(Norm[# - atomcoords[[n]]] + 10^-7) &, rgrid, {3}], {n, 
    Length@atomlist}];
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;

4.788925

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),

CubePlot@ReadCube["~/Downloads/furanpotential_coarse.cube"]

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}]

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

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

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.