Do you have a cluster?

that you can run the jobs on independently?

you are probably right ... :)

yes, I can work on a cluster

Because I find it interesting that you can specify the grid points manually for cubegen...

gaussian.com/g_tech/g_ur/u_cubegen.htm

no, that is my own idea that I want to work on

so you could split up the 10-million grid points into, say, 100 chunks

and in the end have 100 cube files

that you combine in MMA

that would be ...

I've never used that functionality

but I would test it on a smaller molecule, like furan,

and then once you've gotten the syntax right, scale it up to the big molecule

I would not start with the big one, yes.

:-D

Ja, sorry to state the obvious

but in this way you could parallelize cubegen over multiple jobs on the cluster

indeed, that is probably the best solution so far

I found a way to speed it up in MMA, using the fact that you shouldn't have to recalculate the distances each time, there is a finite number of distances (it's a grid....)

so you calculate the distances once, store it as an array

I thought of this, too, but could not come up with a solution.

I also found some errors with my EucDist function. I replaced it with EuclideanDistance

then for each grid point to get the potential you multiply the density array by a portion of the distance array (really 1 divided by it), and use Total to take the sum

Shall I post my current version?

here in chat

I have something working, but still multiplying these huge arrays and taking the Total is more than .01 seconds, so slower than cubegen

let me see how fast I can get my version going, but in the end you are probably better off using cubegen on chunks.

Yeah post here so I can verify that my function gives the same as yours and see if it does end up being faster

vnuc[charges_, molecule_, xyz_] :=
  With[{distances =
     Map[EuclideanDistance[xyz, #] + 10^-7 &, molecule]},
   Total[charges/distances]];
ESP3[cub_] := Block[
  (* Definition der Variablen *)
  {nAtoms, xStart, yStart, zStart,
   xStep, dx, yStep, dy, zStep, dz, molecule, values, data, f, xVals,
   yVals, zVals, xyzVals},
  (* Einlesen der Parameter des Cube-
  Files *)
  {nAtoms, xStart, yStart, zStart} = cub[[3]];
  {xStep, dx} = cub[[4, 1 ;; 2]];
  {yStep, dy} = {cub[[5, 1]], cub[[5, 3]]};

sry ... text is german :)

:-)

living here 2 years and still I don't speak the language

but in Hamburg, almost everyone speaks english anyway

^^

also - do you have a better solution than my 1xN-vector?

Okay, will try it and post an answer after I eat, but it probably won't be faster than cubegen

that was the easiest version that I came up with, instead of using the "real" grid

the version above takes 0.0039918 seconds per point

Why do you add 10^-7 to all of the distances? the distances from the cube file are in atomic units..

Velec will else have to divide through zero

and MMA is not happy with that

right now, I'm at ~0.001 s/point ... wondering what the distances-matrix might bring ... and whatever you also might have been able to improve

at least for my silly test cube

I don't get nearly as fast as you are

and I do ParallelMap[f, xyzVals[[;; 10]]]; // AbsoluteTiming

just on the first 10 data points,

it takes ~20 seconds

@JasonB maybe you are ... I was running my function on a smaller cube

how many cores do you have for parallisation?

@JasonB you can again change vnuc ... I think the EuclideanDistance is a little bit too slow

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

I'm curious about the formula for the grid-based version

Did you write it or find it somewhere?

it seems that that trick with the small distance, can't be right

why that?

if you make that addition smaller, then the potential at any grid point will essentially be due only to the denisty at that point

if you made it 10^-15, then the contribution to the potential at that point would be huge

but the 10^-7 was arbitrary

arbitrary

you could have just as well chosen 10^-15

yes

I thought that 10⁻7 was small enough

small ... little ... ?

right, but the smaller you make it, the closer you are to the divergence you were calculating......

Do you have a source for that formula, or a description anywhere of how to calculate this? I'm having trouble finding a formula online to look at

I wrote it myself

but I had someone confirming it ...

when I first saw it I was assuming that you just skipped over the atom in question, that is, that the summand j was over all grid points except i

@JasonB are you comparing your results with the esp file?

@JasonB I think, when you write ParallelMap[f, xyzVals[[;; 10]]] the results will not be the same as the first ten values from the given esp file, because you were neglecting much of the grid during the calculation as you only enter the first ten values of the grid.

if you want, I can give you a smaller cubefile, then you can use the whole grid for the calculation

It is same though, because what it is doing is determining the potential at those first ten grid points due to all the other grid points. Notice that in vnuc you still give it all the other grid points in the first two arguments. Essentially f[xyz] will always just give the potential at a single point.

I was thinking about it, and I really feel like you need to skip that particular point when calculating the potential, as that's the only way to avoid the infinities. The potential is defined as the force that a positive test charge would feel at that point, so if the positive test charge is at that point then there shouldn't be elctron density at that same point.

but that's just a handwaving argument, so it doesn't seem right

but what I can say is that if you make that arbitrary small value small enough, say 10^-15, then the potential at any point would depend **only ** on the density at that point

The CCL may be a decent place to ask this question

but those are nuclei, they don't spread out to take up the entire space....

touché

You could calculate them separately, the "self potential" and the potential due to other grid points, they would be additive

Currently I try to think about a clever way to use the -5 option for cubegen

@JasonB GaussView does not seem to able to work with the meshes that cubgen produces with "-5" :/