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]}] &)]
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$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 evaluatevnuc[values, xyzVals, {-6., -6., -6.}] dx dy dz
then I get{15.0976, 15.0976, 0.199034}
in return $\endgroup$