Density-based Integral Screening and Incremental Fock Build

[andyj10224]

Description

Implements Density matrix-based integral screening algorithms for Direct SCF, as well as Incremental Fock Build, the process of building a Fock matrix using the difference in the density matrix between SCF iterations.

Reference Paper:
https://onlinelibrary.wiley.com/doi/abs/10.1002/jcc.540100111

Updated Version of PR #2062

Todos

• Provide Massive Speedups to DirectJK Builds for HF and DFT
• Accuracy of Calculations remain uncompromised

Questions

• Is my current method of handling the density matrix for UHF and ROHF ideal?

Checklist

• Passing cc-pVDZ and cc-pVTZ tests for HF and DFT (e_convergence 1.0e-6)
• Passing aug-cc-pVDZ, aug-cc-pVTZ tests for HF and DFT

Status

• Ready for review
• Ready for merge

[JonathonMisiewicz]

Marked as a 1.5 target.

Please fix any merge conflicts. I also recommend that "adding new database files" be moved to a separate PR.

[andyj10224]

Some graphs from running B3LYP on the S22 dimer set (Not S22x7, as it would take too long)

Calculation Wall Times (running the whole database, 8 cores):

Reference: 7186.40 seconds
Density Screening: 4126.07 seconds
Increfock: 3635.78 seconds

[andyj10224]

Marked as a 1.5 target.

Please fix any merge conflicts. I also recommend that "adding new database files" be moved to a separate PR.

Will do, thank you!

[hokru]

I am a bit confused what the graphs show, could you explain them a bit more?

[JonathonMisiewicz]

Also, how have you verified correctness? While I'm not familiar with these methods, "can induce an energy difference of over 2 hartrees even at conservative cutoffs" does not seem safe.

[andyj10224]

I am a bit confused what the graphs show, could you explain them a bit more?

Yes, those graphs represent the distribution of differences expected from using density screening (or density screening + incremental Fock build) compared to not using density screening at all, on the total B3LYP energy of dimers, as well as the interaction energies in the S22 set (not S22x7). Results suggest that there is no harm in using density screening.

[hokru]

What is the number in the top left corner (1e4, 1e-7) and what is plotted on the x-axis?

[andyj10224]

Also, how have you verified correctness? While I'm not familiar with these methods, "can induce an energy difference of over 2 hartrees even at conservative cutoffs" does not seem safe.

I have extensively verified the correctness of this approach, across many different molecules, basis sets, etc. I have never observed an error of "2 Hartrees". The worst I have ever observed is on the order of 1.0e-5 Hartrees, better than density fitting. Here is a link to my "benzene case study".

https://docs.google.com/spreadsheets/d/1HDMZ5PV6GhnK4i68Y1t86VhauTwqMP_X2kII21GCsiY/edit?usp=sharing

I empirically discovered that the best screening threshold given an energy convergence of 1.0e-n is 1.0e-(n+3) for non-augmented basis sets and 1.0e-(n+5) for augmented basis sets.

For example, for an e_convergence of 1.0e-8 in cc-pVTZ, the ideal threshold is 1.0e-11, while it is 1.0e-13 in aug-cc-pVTZ.

[andyj10224]

What is the number in the top left corner (1e4, 1e-7) and what is plotted on the x-axis?

The numbers on the top left corner reference the order of magnitude (e.g. 1e-7 hartrees, or 1.0e-4 kcal/mol).

The values on the x-axis are violin plot represent distribution frequency (It's a violin plot).

[JonathonMisiewicz]

Including the order of magnitude in a corner of the plot is highly non-standard and leads to the confusion I just went through. Choose units that are of the proper order of magnitude, and include any needed decimals on the axis directly. According to your axes, you have 2 hartree error, which is utterly unacceptable. 0.2 microhartree error is much more reasonable.

What you're describing confirms that the error is small. That is different from confirming correctness. Are you able to compare the numbers to some other implementation and show agreement? Is there some rare property that the exact scheme has, which you can numerically reproduce?

[andyj10224]

Including the order of magnitude in a corner of the plot is highly non-standard and leads to the confusion I just went through. Choose units that are of the proper order of magnitude, and include any needed decimals on the axis directly. According to your axes, you have 2 hartree error, which is utterly unacceptable. 0.2 microhartree error is much more reasonable.

What you're describing confirms that the error is small. That is different from confirming correctness. Are you able to compare the numbers to some other implementation and show agreement? Is there some rare property that the exact scheme has, which you can numerically reproduce?

Ah, the energy convergence criteria itself is 1.0e-6, and the errors are well below that number, so I strongly believe that we are safe.

It may not be a good idea to compare to another implementation since every implementation has different tricks thrown in. Rather, if the energy difference is less than the e_convergence criteria, we are definitely safe. Though we could try to reproduce something like MBIS charges to check if we are safe.

[hokru]

Now I also get it :)

Incremental Fock builds can show numerical error accumulation. Adding a reset frequency is sensible, e.g. after 20 incremental builds do a full Fock build from scratch. The actual number should be a (expert) user option.

The incremental Fock build can be compared against a standard build. They have to agree.
Comparing density screening to other programs can probably only be done for HF calculations and even then I would expect that one can only see 'similar errors' and no exact number agreement.

[andyj10224]

Now I also get it :)

Incremental Fock builds can show numerical error accumulation. Adding a reset frequency is sensible, e.g. after 20 incremental builds do a full Fock build from scratch. The actual number should be a (expert) user option.

The incremental Fock build can be compared against a standard build. They have to agree.
Comparing density screening to other programs can probably only be done for HF calculations and even then I would expect that one can only see 'similar errors' and no exact number agreement.

This is already done, as I turn off incremental Fock build after the density drops below 1.0e-5, to reduce the differences, though I could make this a user option. :)

[JonathonMisiewicz]

Wait, what? Then how is Incremental Fock build set to energy tolerance on the order of 1e-9 giving you errors on the order of 1e-7?

[andyj10224]

Wait, what? Then how is Incremental Fock build set to energy tolerance on the order of 1e-9 giving you errors on the order of 1e-7?

Oh no, I mean the SCREENING threshold is set to 1.0e-9, the energy tolerance is 1.0e-6.

[JonathonMisiewicz]

Let me make sure I have this right: you're plotting energy errors on the order of 1e-7, but you only converged the energy to 1e-6? If so, then your plots are pure noise. If not, then what am I getting wrong?

[andyj10224]

Let me make sure I have this right: you're plotting energy errors on the order of 1e-7, but you only converged the energy to 1e-6? If so, then your plots are pure noise. If not, then what am I getting wrong?

Yes sir that is right :) The point of the plots is to show that the errors are irrelavent

[JonathonMisiewicz]

You should not have made those graphs. All you know for sure is that the error is less than 1 micro hartree. You cannot tell the difference between 1.0 e-7 and 1.5 e-7, but people read graphs as if you can tell the difference between your data points.

This is also why you can't just present graphs. You need to be absolutely sure that your labels are clear, and that somebody who isn't you will be able to figure out what the graph means. If you need to add a sentence to explain them, do so.

I request additional benchmarks so we can get more precise estimates on how much error these techniques introduce. Can you increase energy convergence to 1e-10?

Because this PR isn't coming in until 1.5 anyways, I'm going to turn my attention to other things for a while.

[andyj10224]

You should not have made those graphs. All you know for sure is that the error is less than 1 micro hartree. You cannot tell the difference between 1.0 e-7 and 1.5 e-7, but people read graphs as if you can tell the difference between your data points.

This is also why you can't just present graphs. You need to be absolutely sure that your labels are clear, and that somebody who isn't you will be able to figure out what the graph means. If you need to add a sentence to explain them, do so.

I request additional benchmarks so we can get more precise estimates on how much error these techniques introduce. Can you increase energy convergence to 1e-10?

Because this PR isn't coming in until 1.5 anyways, I'm going to turn my attention to other things for a while.

Ah, I am learning a lot. Good teachers like you will prepare me well for grad school :)

[JonathonMisiewicz]

I feel like I need a copy of the original Ahlrichs paper before I can review this PR - I don't understand the "big picture."

While I'm waiting for the library to get that to me, here are some initial comments.

[jturney]

Super nice! Thanks for the feature!

[andysim]

This is really outstanding work - kudos for getting it all figured out, and implemented in some of the most difficult parts of the code in terms of control flow. My concern before was that the list of max densities is an (N^2) quantity and could get large for larger systems, while storing it only for those shell pairs that survive the Schwarz/CSAM screening would not. This is not really a valid concern because we have overlap matrices and other O(N^2) quantities in core anyway, so the implementation is great as-is. I added a suggestion or two to add the CSAM reference because I can never remember which paper it comes from. LGTM!

[JonathonMisiewicz]

Per Lori request, we have a hold on merging this in until David Poole can look this over.

[davpoolechem]

I just finished looking at the code. I left a couple of comments to be reviewed before I approve these changes.

[andyj10224]

I addressed all of David's comments. Ready to merge.

[JonathonMisiewicz]

I'd rather have @davpoolechem officially sign off on this first. It'll be an hour before the test suite check finishes, anyway.

[davpoolechem]

Andy and I talked a bit about the comments that I made, and he has sufficiently addressed all of them. I give my seal of approval to this code!