COSX Reformulation for Negative Grid Weights

[davpoolechem]

Description

This PR is a follow-up to #2906, and what can be considered an official solution to the issue discussed in #2890. The current issue is that COSX does not work with certain grid configurations - specifically, it does not work with grids that have negative grid weights, due to the use of an intermediate matrix in COSX that uses the square root of negative grid weights. #2906 "fixed" this issue by having COSX throw an exception when grids with negative weights were encountered. This PR provides a COSX reformulation that allows COSX to work with negative grid weights.

The reformulation does two primary things:

• The $X$ matrix (Eq. 4 in Neese 2009) is redefined as $X_{\kappa g} = \sqrt{|w_{g}|} \kappa(r_{g})$. In words, $X_{\kappa g}$ now uses the square root of the magnitude of the weights, instead of the square root of the raw weights.
• To correct for the above, the computation of the $G$ matrix (Eq. 7 in Neese 2009) is performed as $G_{\nu g} = \sum_{\tau} \text{sign}(w_{g}) A_{\nu \tau} (r_{g}) F_{\tau g}$ . In words, the sign of the corresponding grid weights are included in the formation of $G_{\nu g}$, when $A_{\nu \tau}$ and $F_{\tau g}$ are contracted.

One other thing comes out as a consequence of this:

• For overlap fitting, the numerical overlap matrix $S_{N}$ (Eq. 13 in Izsák 2011) is computed as $S_{N} = X * X_{\text{alt}}^{T}$, where $X_{\kappa g} = \sqrt{|w_{g}|} \kappa(r_{g})$ (the same as the first bullet point in the reformulation above), and $X_{\text{alt},\kappa g} = \text{sign}(w_{g}) \sqrt{|w_{g}|} \kappa(r_{g})$. Essentially, $S_{N}$ now uses two variants of the $X$ matrix in its formulation, one of which folds the grid weight sign into itself.

User API & Changelog headlines

• COSX can now be used with a wider variety of grids, as a bug preventing COSX to be used with specific grid configurations has been fixed.

Dev notes & details

• The COSX implementation has been reformulated to enable calculations with grids containing negative grid weights.

Questions

• Is the current way of handling computation of the numerical overlap matrix acceptable? Currently, two X matrices are used to form the numerical overlap matrix. One goal of folding the grid weight sign into the formation of G was to prevent the use of multiple X matrices. However, I don't think that can be done here without removing the call to linalg::doublet in the compute_numeric_overlap function. One could probably do a slightly-modified, manually-implemented matrix multiply with the grid weight folded in, and construct $S_{N}$ with a single $X$ matrix that way. But that comes at the cost of not utilizing BLAS. What does everyone consider preferable?
• What, exactly, is the best way to test this? I would like to add a test to test_dfjcosk.py for a new grid weight option with negative weights. Unfortunately, it seems that other codes that have COSX just don't have grid options that lead to negative grid weights. For now, I've been comparing against Psi4 COSX calculations using similar grid sizes and ensuring the calculations have similar energies, but I certainly wouldn't consider that robust.

Checklist

• Tests added for any new features
• All or relevant fraction of full tests run

Status

• Ready for review
• Ready for merge

[susilehtola]

I think this should not be a WARNING but INFO (does this exist in Psi4?)

[davpoolechem]

I'm happy to switch to INFO, but Psi4 seems to trend towards using WARNING instead of INFO for printout like this, at least from my observation.

[davpoolechem]

Added! I'll leave the conversation open, though, in case others want to comment on INFO vs. WARNING.

[loriab]

I think this should not be a WARNING but INFO (does this exist in Psi4?)

There's no proper logging from c-side. :-( Only PSIEXCEPTION that stops execution and Printf("WARNING/INFO/ETC that goes to the output file. Admittedly not great. I don't think the situation has been seriously examined for improvement since pybind11 (and having python on the outside of psi4) was young.

in case others want to comment on INFO vs. WARNING.

I agree it's probably INFO at this point (rounding the grid up to down to a pure-positive one would be WARNING), but there's not standardization.

[davpoolechem]

I think this should not be a WARNING but INFO (does this exist in Psi4?)

There's no proper logging from c-side. :-( Only PSIEXCEPTION that stops execution and Printf("WARNING/INFO/ETC that goes to the output file. Admittedly not great. I don't think the situation has been seriously examined for improvement since pybind11 (and having python on the outside of psi4) was young.

in case others want to comment on INFO vs. WARNING.

I agree it's probably INFO at this point (rounding the grid up to down to a pure-positive one would be WARNING), but there's not standardization.

All good to know! Maybe such logging is something that could be more seriously looked at in the future. It wouldn't be a bad functionality to have.

[susilehtola]

Looks pretty good, please fix the issues above

[davpoolechem]

Looks pretty good, please fix the issues above

Thanks for the review! I think everything you suggested should be covered now.

[susilehtola]

Suggested change

	X_block_nosignp[p][k] = point_values->get(p, k) * std::sqrt(std::fabs(w[p]));
	X_block_nosignp[p][k] = point_values->get(p, k) * std::sqrt(std::abs(w[p]));

[susilehtola]

I'm good with the way of computing the overlap matrix. One does want to use BLAS and this just means that you need an extra copy of the basis function value matrix.

[davpoolechem]

I'm good with the way of computing the overlap matrix. One does want to use BLAS and this just means that you need an extra copy of the basis function value matrix.

Yeah, I was thinking it would be preferable to use BLAS, as well.

[davpoolechem]

So with #2955 (SplitJK) also being a thing, I think it would be preferable to get that one in first, and then rebase this one to adjust for the changes.