Porcelain

[FarnazH]

• A function which takes atomic numbers, atomic coordinates, and some description of grid and returns an instance of MolecularGrid class. (Note: here is where tv-*.txt points stored in needed; check Atomic Grid Class #11)

[PaulWAyers]

Also:

• Something to expand a function into its spherical atomic composition. (Needed for the next issues.)
• A function that converts a function on a grid to its potential (Poisson solve) or converts a potential on a grid to a function (differentiation). (For arbitrary functions.)
• Second-order differential equations solve. (Already existing, needs to port over.)
• Multipole moments to arbitrary order for each atom, and also for the molecule with an arbitrary center.

These things may use splines (C++) or else could be done with scipy splines (slower probably, less flexible for sure).

A lot of this should be array operations in Python. Probalby no C required anywhere though compatibility with Toon's (in progress) cubic splines should be maintained.

[PaulWAyers]

For the ODE solver, scipy.integrate.odeint

You need to convert higher-order ODE into a first-order one to use this but that can always be done.

See, for example,
https://math.stackexchange.com/questions/70993/how-to-express-a-2nd-order-ode-as-1st-order-odes

[tovrstra]

Regarding order of implementation: the poisson solver also needs the ODE solver, so it needs to be implemented first. I would suggest to split the ODE, Poisson solver, multipole decomposition and evaluation into separate issues or pull requests. Doing it all together is far too much for one PR.

[FarnazH]

Regarding order of implementation: the poisson solver also needs the ODE solver, so it needs to be implemented first. I would suggest to split the ODE, Poisson solver, multipole decomposition and evaluation into separate issues or pull requests. Doing it all together is far too much for one PR.

@tczorro looking at this closed issue again, I was wondering whether multipole decomposition and evaluation is implemented?

[PaulWAyers]

• I don't see evaluating the Laplacian here. That's important as the "inverse" of the Poisson solver.
• Spherical average. Must be around an atomic center.
• Multipole moments (for an arbitrary center)

1. The Laplacian evaluation is relatively straightforward and may be implemented, but I couldn't find it.
2. Spherical Average notes
3. Moments notes

[PaulWAyers]

Spherical average notes updated, hat tip especially @Ali-Tehrani and @FarnazH for clarifications.
Sphericalaverage.md
Sphericalaverage.pdf

[PaulWAyers]

To generate spherical moments, we can use the code in denspart (which generates spherical harmonics).
https://github.com/theochem/denspart/blob/082cf5df6e0438cef06589cb9d371ca34dbee621/denspart/properties.py#L129

An alternative is to generate the real spherical harmonics directly using scipy. We can do this easily, and it makes the code a lot shorter (and more robust?) but we don't benefit from the recursion when evaluating the spherical harmonics. This is probably not a huge deal performance-wise, as evaluating the spherical harmonics is much faster than any function we are likely to be evaluating the moments of.

[PaulWAyers]

Note on multipole moments are attached here.
MomentsGrid.md
MomentsGrid.pdf

[PaulWAyers]

I'm confused, though. Looking through things it seems we had already translated the old moments code into Python for grid.
https://github.com/theochem/grid/blob/713527cd6e2bfdb1d467ae8f58e962b90b49d160/grid/moments.py

@tczorro am I missing something? Can you explain why this file is no longer part of the repo?

[tczorro]

@PaulWAyers I think most those changes are done by Toon. So I am very aware of those.

[FarnazH]

@tczorro and @PaulWAyers it seems like grid/grid/moment.py was removed when the grid directory ported from HORTON was removed. I might be wrong, but I remember that @tczorro started the new grid in grid/src/grid (instead of modifying the one from HORTON), and at some point, we removed the files extracted from HORTON because all features were implemented. In any event, considering that we had/have a version of moments, I hope reimplementing it would be fast.

[PaulWAyers]

Talking with @tovrstra @evohringer and @FarnazH yesterday, it seems best to add moments as an integration option, in the base class. The only caveat is that we would need to be very careful with periodic boundary conditions. That's something that can be, I think, handled in some elegant way (that I haven't figured out or even really thought about yet). (It's a standard problem scientifically, but code-wise it is trickier.)

Also, we should perhaps take spherical harmonics (by duplication) from the new denspart. The reason is that scipy spherical harmonics seem to not use a recursion. That is not a huge deal for the moments (where you do not compute to very high order) but when we are doing interpolation/differential-eq-solving, we compute spherical harmonics up to very high order (degree of the angular grid // 2). Note that the normalization used there is appropriate for moments (good for that) but not great for other bits of the code.
https://github.com/theochem/denspart/blob/082cf5df6e0438cef06589cb9d371ca34dbee621/denspart/properties.py#L129

Another alternative is to use shtools, which may be faster (Fortran code with Python wrapper). They will compute all the spherical harmonics up to an order, but don't seem to be vectorized (you can only compute one point at a time). They support real or complex spherical harmonics with a wide variety of normalizations.
https://shtools.github.io/SHTOOLS/pyspharm.html

A final alternative would to be to code the recursion and use some sort of dynamic programming, but I think the memory obligations from the caching would be crazy.

[PaulWAyers]

Reminder to self @PaulWAyers to figure out the algorithm in denspart for solid harmonics so we can document it.

[PaulWAyers]

For moments, we can use as tests a function (e.g. a Gaussian) that is centered at the origin and also one centered off the origin. The change in solid-harmonic coefficients when you shift the origin of a moment expansion is analytically computable.

• There is code in Julia that shifts the origin of a solid harmonic expansion

There are also some analytic formulas available.

• see page 66
• Ruedenberg found a classical result, which is simplified in Eq. (8) of this reference
• There is also a result in this reference
• A tabulation of the coefficients for translation/rotation is available in the supplementary material of this reference

[PaulWAyers]

Talking with @tovrstra @evohringer and @FarnazH yesterday, it seems best to add moments as an integration option, in the base class. The only caveat is that we would need to be very careful with periodic boundary conditions. That's something that can be, I think, handled in some elegant way (that I haven't figured out or even really thought about yet). (It's a standard problem scientifically, but code-wise it is trickier.)

Also, we should perhaps take spherical harmonics (by duplication) from the new denspart. The reason is that scipy spherical harmonics seem to not use a recursion. That is not a huge deal for the moments (where you do not compute to very high order) but when we are doing interpolation/differential-eq-solving, we compute spherical harmonics up to very high order (degree of the angular grid // 2). Note that the normalization used there is appropriate for moments (good for that) but not great for other bits of the code. https://github.com/theochem/denspart/blob/082cf5df6e0438cef06589cb9d371ca34dbee621/denspart/properties.py#L129

Another alternative is to use shtools, which may be faster (Fortran code with Python wrapper). They will compute all the spherical harmonics up to an order, but don't seem to be vectorized (you can only compute one point at a time). They support real or complex spherical harmonics with a wide variety of normalizations. https://shtools.github.io/SHTOOLS/pyspharm.html

A final alternative would to be to code the recursion and use some sort of dynamic programming, but I think the memory obligations from the caching would be crazy.

We implemented a recursion for spherical harmonics from the literature.
Journal of Geodesy 76.5 (2002): 279-299.