[REQUEST] Calculation of density matrices, and spin density matrices

[thomaspigeon]

Is your feature request related to a problem? Please describe.

[thomaspigeon]

I use gbasis to evaluate the spin densities of a molecule. The basis informations are first read from a wavefunction file with the help of iodata and the the basis list of tuples is created from the gbasis wrapper with iodata. The thing is that to evaluate the spin densities values on a grid, I need the spin density matrices (called one one_dm in the documentation on how to use the function evaluate_density) so I have to recompute the spin density matrices from the alpha/beta coeffs and occupation numbers. I currently do it in a short function. I am wondering whether it would be useful to add this in gbasis, maybe as an option in the iodata wrapper(or maybe it would be more logical to add it in iodata package). What do you think about it ? I could try to make a pull request if you think this is relevant.

[FarnazH]

Thanks for your comment and suggestions @thomaspigeon. We have thought about this issue and tried multiple things. The one_dm used to be computed by iodata in the past, but it was removed because iodata is intended to be only a parser module for loading/dumping various file formats (i.e., the one_dm might be loaded in one of the formats, but we decided not to include any computations in iodata even for sth as simple as one_dm)... It also doesn't really belong to gbasis, because you need molecular orbitals information to compute density matrices which is not the scope of gabsis... So, we believe the one_dm computation should be handled in the glue/porcelain code that puts iodata and gbasis together. This porcelain code would be part of horton (which is not available yet) and chemtools (which is available, but currently the public chemtools still uses horton2 and it is being changed to use iodata, gbasis, grid & denspart modules). In the meanwhile, I think your function for computing one_dm is the best solution. Hope this helps.

[PaulWAyers]

We already implicitly have the density-matrix and density as "zeroth derivatives" in gbasis. So I'm not sure we shouldn't support it.

It's a high-level design decision, really related to how we are going to piece things together. If HORTON is a library of modules with the porcelain that allows them to talk to each other, it is plausible to put this in gbasis. It would be nice, probably for numerical efficiency as well as for convenience, if gbasis computed the density and its derivatives, and the kinetic energy density and other density-matrix derivatives, and the Fock matrix and such. That's sort of the way it is in the current spec.

If not, then there needs to be a rather substantial porcelain module that would take evaluations of basis functions and their derivatives from gbasis and puts it together with input from other sources to construct the DMs. I guess it boils down to "how finely sliced" the modules are.

I have some preference for having gbasis handling the bulk of the mathematical operations related to manipulating Gaussian basis sets so that chemtools and meanfield and fanpy can focus on other pieces. But this does risk making gbasis large, which means that gbasis needs to be "internally modular" to a great extent too.

In any event, I'm not sure it is clear-cut. But in the new more modular HORTON, there are things that are repeated between modules because of the desire to keep modularity (e.g., overlap integrals in iodata, some Gaussian basis set things in grid, etc.) So the DM can be several places if it is convenient.

[PaulWAyers]

Getting 1DM and 2DM, including spin-resolution, should be a feature of HORTON. Those are then passed to gbasis to compute properties.

[PaulWAyers]

@sanchezw17 needs this for his work, so will tackle implementing it (at the Slater determinant/uncorrelated level for now).