Atomic Density PT

[kgasperich]

PR Type

• [ ] Breaking change
• [x] Feature
• [ ] Patch

Brief Description

This PR adds a new PT for obtaining atomic densities. related PRs: https://github.com/NWChemEx-Project/NWChemEx/pull/99 https://github.com/NWChemEx-Project/SCF/pull/182

This initial implementation just takes an AtomicBasisSet<double> and returns a std::vector<double>. I had considered returning something with type chemist::PeriodicTable::atom_dm_t so that this PT would always remain consistent with however the things are stored in the PeriodicTable. I can change this if anyone thinks it would be better that way.

When using a module with this PT to retrieve data from a PeriodicTable, a PeriodicTable is also needed, so in the module where I implement that I've added the PeriodicTable as another input.

For the atomic SCF module, we'll probably want the total (possibly non-integer) number of electrons (or equivalently a total charge) as another input. We might also want a PeriodicTable in some cases (to get some reference atomic configuration). I could just add the PeriodicTable to the PT inputs (rather than adding it in both of these module implementations), but it seems more general/flexible to not include it (although I could have it be an optional input that defaults to a default-constructed PeriodicTable).

Not In Scope

PR Checklist

• [ ] Documentation (comments, API, user, design, and developer as applicable)
  • API
  • User
  • Design
  • Developer
• [ ] Testing (unit, integration, performance, as applicable)
  • Unit
  • Integration
  • Performance
• [ ] Code adheres to organization standards
  • C++
  • Python
  • Doxygen Conventions
  • reST Conventions

TODOs

• [ ] Task 1
• [ ] Task 2

[ryanmrichard]

Is this supposed to abstract the whole process of getting the initial AO density, or just the AOs to initial AO guess for your atomic SCF?

If the former, I recommend taking an electronic Hamiltonian and returning a TensorWrapper. Reasons:

• H_e has the number of electrons, AOBasisSet doesn't
• AOBasisSet doesn't have the atomic numbers of the atoms, H_e will (IIRC at present H_e has their charges, which is not necessarily the same thing, but for our initial use cases should work; @jwaldrop107 and I talked yesterday about a slight tweak to how we define point charges/nuclei which should alleviate this problem)
• For @yzhang-23's module, the AOBasisSet you pass in to the overall SCF will in general not be the basis the densities come back in. That places you in an awkward spot of actually building the basis on the fly.
• There's supposed to be a module (I think @jwaldrop107 is writing this too) which goes from Molecule to AOBasisSet. Inside your atomic SCF guess you can use that module to get the basis for the atomic SCF from H_e. You can imagine scenarios where the user may want you to actually do this step in a smaller basis than the one wanted for the entire SCF. The module approach allows the user this control.
• My understanding is the return of this module is going to be used to build the initial Fock operator. If that's the case the TensorWrapper return makes this easier.

[kgasperich]

Is this supposed to abstract the whole process of getting the initial AO density, or just the AOs to initial AO guess for your atomic SCF?

I was thinking that it would just return the AO density for some basis, but I should have thought through the details a bit more. I should have marked this as a draft/work-in-progress. This initial iteration was primarily aimed at just reproducing what @yzhang-23 had done to retrieve the atomic density, but now inside a submodule call. Now that I'm taking another look at it, it clearly doesn't pass enough information through the API to handle the atomic SCF (or cases where we don't have a saved DM for this atom/basis).

For retrieving a stored DM, we need more info about the basis. If we don't have the desired basis, we need to be able to choose one of the bases for which we have a stored DM, and then either (1) project the density back onto the desired basis and return that density or (2) return the DM along with the basis in which it's represented so that we can do the Fock build correctly (based on discussions from the last two Wednesday meetings, I think the plan is to do option (2)). chemist::Density already has the density and a basis, so we could return an instance of that class (I'd just need to wrap the flattened vector as a tensor with the right shape).

For the atomic SCF, we'll need: atomic number; AO basis; number of electrons of each angular momentum (currently, we can get this via PeriodicTable if we have the total charge or total number of electrons). This is pretty much the same input as AtomicDensityLookup, but with an additional double for the charge or number of electrons. This doesn't include any information about pseudopotentials (or any kind of external potential), so if we want to be able to include that then the Hamiltonian would be a better input option.

The AtomicBasisSet has an atomic number, but I haven't looked at the constructor to see if this is necessarily the same as the charge of the ion in that position (in something like a BSSE counterpoise calculation, there will be centers corresponding to particular atoms, but without an actual nucleus at that center, so I don't know whether the atomic number from the basis set will be the number of the corresponding atom or zero (or in some cases I've seen negative numbers used to represent dummy centers, e.g. -8 would be a ghost atom with oxygen's basis functions)).

I can make another pass through this and try to make a more general API using H_e and returning a Density.

[ryanmrichard]

For retrieving a stored DM, we need more info about the basis. If we don't have the desired basis, we need to be able to choose one of the bases for which we have a stored DM, and then either (1) project the density back onto the desired basis and return that density or (2) return the DM along with the basis in which it's represented so that we can do the Fock build correctly (based on discussions from the last two Wednesday meetings, I think the plan is to do option (2)). chemist::Density already has the density and a basis, so we could return an instance of that class (I'd just need to wrap the flattened vector as a tensor with the right shape).

You're right. We need the basis set to build the density in. So I'd recommend taking H_e and that basis.

The AtomicBasisSet has an atomic number, but I haven't looked at the constructor to see if this is necessarily the same as the charge of the ion in that position (in something like a BSSE counterpoise calculation, there will be centers corresponding to particular atoms, but without an actual nucleus at that center, so I don't know whether the atomic number from the basis set will be the number of the corresponding atom or zero (or in some cases I've seen negative numbers used to represent dummy centers, e.g. -8 would be a ghost atom with oxygen's basis functions)).

IIRC we added this as a work around for something; I'd be inclined to remove it since it becomes ill-defined for things like mid-bond functions, which we should try to support. Atomic numbers should be coming from the operators involving nuclei.