atomic quadrature grid blocking [part 1] #2336

hokru · 2021-10-26T12:11:11Z

Description

This PR introduces an atomic blocking scheme for quadrature grid points. All grid points in a BlockOPoints object belong to a singular parent atom.
Part 1 contains just the basic feature to get it our for people depending on it. Optimization for speed is yet to come.

New options:

DFT_BLOCK_SCHEME = ATOMIC (unique blocks of grid points for each atom)
DFT_REMOVE_DISTANT_POINTS (new flag for existing functionality)

primary C++ feature:

grid->atomic_blocks()[N_ATOM][N_BLOCKS] additionally to grid->blocks()[N_BLOCKS] # provides all grid points for an atom
block->parent_atom() # atom the current block belongs to.

fixes

collocation size estimate in naive gridblocker was wrong.
removes unused index vector

note

Automated formatting of the files with clang-format made unrelated changes.

Questions

The ugly code here is because of an issue with the BlockOPoints object. I'd like advice how to handle this better.

Checklist

Tests added for any new features

Status

Ready for review
Ready for merge

Usage

Normal loop structure

for (size_t Q = 0; Q < grid_->blocks().size(); Q++) {
  size_t parent_atom_ = grid_->blocks()[Q]->parent_atom();
  .
  .
  .
}

Looping over atoms and their blocks.

# openmp note: for best performance the atom and block loop could possibly
# be collapsed into a singular loop using `collapse(2)`
for (size_t i = 0; i < grids_.size(); i++) { # here size = number of atoms
    for (size_t Q = 0; Q < grid_->atomic_blocks()[A].size(); Q++) {

        // access a block
        std::shared_ptr<BlockOPoints> block = grid_->atomic_blocks()[A][Q];
  .
  .
  .
}

JonathonMisiewicz

I don't know the DFT infrastructure at all, so can you give me a big picture overview of what this PR is doing? I know that DFT relies on numerical integration, and I know that a block contains numerical integration points, but I don't know what a block is used for.

I have some structural concerns with this PR, but I don't want to suggest something without knowing more about what this PR is trying to accomplish.

JonathonMisiewicz · 2021-10-26T13:34:52Z

psi4/src/psi4/libfock/cubature.h

@@ -399,6 +417,10 @@ class BlockOPoints {
    size_t index() const { return index_; }
    /// Print a trace of this BlockOPoints
    void print(std::string out_fname = "outfile", int print = 2);
+    /// set parent atom for the current block


Suggested change

/// set parent atom for the current block

/// Set parent atom for the current block

JonathonMisiewicz · 2021-10-26T13:37:13Z

psi4/src/psi4/libfock/cubature.h

@@ -399,6 +417,10 @@ class BlockOPoints {
    size_t index() const { return index_; }
    /// Print a trace of this BlockOPoints
    void print(std::string out_fname = "outfile", int print = 2);
+    /// set parent atom for the current block
+    void set_parent_atom(size_t atom) {parent_atom_ = atom;};
+    /// parent atom if the current block


Suggested change

/// parent atom if the current block

/// Parent atom if the current block

JonathonMisiewicz · 2021-10-26T13:44:30Z

psi4/src/psi4/libfock/cubature.h

@@ -377,6 +392,9 @@ class BlockOPoints {
    /// Bounding radius of the BlockOPoints
    double R_;

+    /// parent atom of this BlockOPoints
+    size_t parent_atom_;


Why are you adding this function, and the associated getters/setters? The only place this is used is for printing in your test.

Are all BlockOPoints atom-centered?

psi4/src/psi4/libfock/cubature.h

psi4/src/psi4/libfock/cubature.cc

JonathonMisiewicz · 2021-10-26T14:11:24Z

psi4/src/psi4/libfock/cubature.cc

            checkBrian();
        }
    }
 #endif

    npoints_ = 0;
-    for (int i = 0; i < grid.size(); ++i) {
-        npoints_ += grid[i].size();
+    for (int i = 0; i < atomic_grids_.size(); ++i) {


Replace this with a range-based for loop.

JonathonMisiewicz · 2021-10-26T14:12:12Z

psi4/src/psi4/libfock/cubature.cc

+            y_[grid_vector_index] = atomic_grids_[i][j].y;
+            z_[grid_vector_index] = atomic_grids_[i][j].z;
+            w_[grid_vector_index] = atomic_grids_[i][j].w;
+            index_[i] = grid_vector_index;  // TODO: It's not used and wrong->remove.


If index_ isn't used, get rid of it now and use range-based for loops.

~~I stand corrected. octree uses it for debug printing. I still think it is wrong, though, or I just don't understand it.~~

JonathonMisiewicz · 2021-10-26T14:14:36Z

psi4/src/psi4/libfock/cubature.cc

                    mp.w *= nuc.computeNuclearWeight(mp, A, stratmannCutoff);  // This ain't gonna fly. Must abate this
                                                                               // mickey mouse a most rikky tikki tavi.


While I'm glad that Rob Parrish had fun, I have no clue what he's talking about. Remove this comment.

Suggested change

mp.w *= nuc.computeNuclearWeight(mp, A, stratmannCutoff); // This ain't gonna fly. Must abate this

// mickey mouse a most rikky tikki tavi.

mp.w *= nuc.computeNuclearWeight(mp, A, stratmannCutoff);

JonathonMisiewicz · 2021-10-26T14:14:47Z

psi4/src/psi4/libfock/cubature.cc

@@ -3836,7 +3841,9 @@ void MolecularGrid::buildGridFromOptions(MolecularGridOptions const &opt) {
                mp.w *= nuc.computeNuclearWeight(
                    mp, A,
                    stratmannCutoff);  // This ain't gonna fly. Must abate this mickey mouse a most rikky tikki tavi.


Suggested change

stratmannCutoff); // This ain't gonna fly. Must abate this mickey mouse a most rikky tikki tavi.

stratmannCutoff);

JonathonMisiewicz · 2021-10-26T14:16:38Z

psi4/src/psi4/libfock/cubature.cc

    postProcess(extents, max_points, min_points, max_radius);
 }

+// REFACTOR NOTE: PS grids are not used. Incomplete use of MolecularGridOptions.


I do not understand the second sentence.

hokru · 2021-10-26T20:39:13Z

I don't know the DFT infrastructure at all, so can you give me a big picture overview of what this PR is doing? I know that DFT relies on numerical integration, and I know that a block contains numerical integration points, but I don't know what a block is used for.

Instead of making a loop over all grid points that calculates basis function values and XC contributions one by one , the points are grouped together into blocks for efficiency reasons.
Currently there is no connection between a single grid point and the atom it originally belongs to. One block of points can even contain grid points from multiple atoms. This connection is lost immediately after the grid is constructed. Adding it back is a bit awkward and e.g. remove_distant_points(extents_->maxR()); needs to be rewritten with the new data structure.

This feature is needed for algorithms or methods that look at contributions from atoms when looping over the grid points. Right now, those new things are not needed, and thus not used anywhere, but they will be used for COSX and ddCOSMO.

So you might want to access for a given atom all blocks of grid points, where a block has only grid points from that atom. This is what atomic_blocks provides. However that leads to a double-loop (see examples) and possibly worse parallelization. For a simple loop over all block, like what is used now, you want to know the atom to which the current block belongs (-> parent_atom) to collect values by atom. Two options the developers can chose from.

This is part 1 to get it the basics out faster and allow COSX and ddCOSMO development to continue.

cubature.cc still needs a big cleanup, but that has to wait.

susilehtola · 2021-10-26T22:55:18Z

This feature is needed for algorithms or methods that look at contributions from atoms when looping over the grid points. Right now, those new things are not needed, and thus not used anywhere, but they will be used for COSX and ddCOSMO.

It's also needed for generalized Pipek-Mezey, which is a better localization method than what is currently implemented in Psi4.

JonathonMisiewicz · 2021-10-27T10:01:17Z

The other point I missed the first time was when you said grid->atomic_blocks()[N_ATOM][N_BLOCKS] was the primary feature. One word, but an important one.

For future, it's clearer to put the unambiguous statement of the primary feature before technical notes (like on clang-format and options) when writing the PR description. I was expecting the notes to be in decreasing order of importance.

psi4/src/psi4/libfock/cubature.cc

JonathonMisiewicz · 2021-10-26T20:58:51Z

psi4/src/psi4/libfock/cubature.cc

+    opt.bench = full_int_options["BENCH"];
+
+    // Until the subroutine is fixed we need to disable the screening for atomic blocking
+    // Otherwise there is little reason to distable it.


Suggested change

// Otherwise there is little reason to distable it.

The comment above is clear enough that we can re-enable it when some bug is fixed.

psi4/src/psi4/libfock/cubature.cc

tests/pytests/test_dft_blocking.py

JonathonMisiewicz · 2021-10-27T10:15:04Z

psi4/src/psi4/libfock/cubature.cc

+                                         "PRINT",
+                                         "DEBUG",
+                                         "BENCH",
+                                         "DFT_REMOVE_DISTANT_POINTS"};


debug, bench, and dft_remove_distant_points should be bool rather than int. The first two were existing keywords, so I can agree to not convert those until a future PR. But dft_remove_distant_points should be made a bool now.

I'll make DFT_REMOVE_DISTANT_POINTS a global options for now. Eventually all grid options should become build-specific. For this one would need to add double and bool std::maps to the grid constructor. A more universal option would be better. I added a note in the code:

JonathonMisiewicz · 2021-10-27T10:54:41Z

psi4/src/psi4/libfock/cubature.cc

+            // below segfaults, cannot use grids_ directly.
+            // auto bop =std::make_shared<BlockOPoints>(block_index, n, &grids_[i][grid_index].x,
+            // &grids_[i][grid_index].y, &grids_[i][grid_index].z, &grids_[i][grid_index].w, extents_);
+            if (bop->local_nbf()) {


What would it even mean for this to be zero?

I took this from the octree code. I think there is a chance that no significant basis functions remain after screening within the BlockOPoints.

psi4/src/psi4/libfock/cubature.cc

JonathonMisiewicz · 2021-10-27T11:00:55Z

psi4/src/psi4/libfock/gridblocker.h

@@ -55,6 +56,8 @@ class GridBlocker {
    double const* z_ref_;
    double const* w_ref_;
    int const* index_ref_;
+    std::vector<std::vector<MassPoint>> grids_;


Members that are only needed by one subclass should only be defined for that subclass.

I seem to need it in the base class for this blocker->set_atomic_grids(atomic_grids_); here because blocker is not actually the subclass object. (as far as I understand the compiler)

JonathonMisiewicz

This looks much better. I have a few comments before I can approve this. 👍 for getting rid of index_ and index_ref_.

psi4/src/psi4/libfock/cubature.cc

psi4/src/psi4/libfock/cubature.h

Stale.

…r distant point removal

JonathonMisiewicz

As somebody new to the DFT code, LGTM. Thanks for the PR! It looks like this will enable some useful things.

jturney

This looks good!

maxscheurer

LGTM. 👍 for the tests!

JonathonMisiewicz reviewed Oct 26, 2021

View reviewed changes

JonathonMisiewicz previously requested changes Oct 27, 2021

View reviewed changes

hokru force-pushed the atomic_grid1 branch 4 times, most recently from f24757a to bfea8b4 Compare November 15, 2021 13:37

susilehtola mentioned this pull request Dec 8, 2021

Move IBO localizer into libmints #2379

Open

hokru force-pushed the atomic_grid1 branch from a16af3d to 36e0a9e Compare December 9, 2021 12:23

initial implementation of atom-wise blocking of quadrature points

d610ea2

hokru force-pushed the atomic_grid1 branch from 36e0a9e to d610ea2 Compare December 9, 2021 12:51

JonathonMisiewicz reviewed Dec 9, 2021

View reviewed changes

review requests. And shrinking atomic_blocks to correct capacity afte…

ea55bc4

…r distant point removal

JonathonMisiewicz approved these changes Dec 10, 2021

View reviewed changes

jturney approved these changes Dec 10, 2021

View reviewed changes

maxscheurer approved these changes Feb 4, 2022

View reviewed changes

maxscheurer merged commit 6141a69 into psi4:master Feb 4, 2022

JonathonMisiewicz added this to the Psi4 1.6 milestone Feb 11, 2022

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

atomic quadrature grid blocking [part 1] #2336

atomic quadrature grid blocking [part 1] #2336

hokru commented Oct 26, 2021 •

edited

Loading

JonathonMisiewicz left a comment

JonathonMisiewicz Oct 26, 2021

JonathonMisiewicz Oct 26, 2021

JonathonMisiewicz Oct 26, 2021

JonathonMisiewicz Oct 26, 2021

JonathonMisiewicz Oct 26, 2021

hokru Oct 28, 2021 •

edited

Loading

JonathonMisiewicz Oct 26, 2021

JonathonMisiewicz Oct 26, 2021

JonathonMisiewicz Oct 26, 2021

hokru commented Oct 26, 2021

susilehtola commented Oct 26, 2021

JonathonMisiewicz commented Oct 27, 2021

JonathonMisiewicz Oct 26, 2021

JonathonMisiewicz Oct 27, 2021

hokru Oct 28, 2021

JonathonMisiewicz Oct 27, 2021

hokru Oct 28, 2021

JonathonMisiewicz Oct 27, 2021

hokru Oct 29, 2021

JonathonMisiewicz left a comment

JonathonMisiewicz left a comment

jturney left a comment

maxscheurer left a comment

	/// set parent atom for the current block
	/// Set parent atom for the current block

	/// parent atom if the current block
	/// Parent atom if the current block

		mp.w *= nuc.computeNuclearWeight(mp, A, stratmannCutoff); // This ain't gonna fly. Must abate this
		// mickey mouse a most rikky tikki tavi.

	stratmannCutoff); // This ain't gonna fly. Must abate this mickey mouse a most rikky tikki tavi.
	stratmannCutoff);

atomic quadrature grid blocking [part 1] #2336

atomic quadrature grid blocking [part 1] #2336

Conversation

hokru commented Oct 26, 2021 • edited Loading

Description

New options:

primary C++ feature:

fixes

note

Questions

Checklist

Status

Usage

JonathonMisiewicz left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

hokru Oct 28, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

hokru commented Oct 26, 2021

susilehtola commented Oct 26, 2021

JonathonMisiewicz commented Oct 27, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

JonathonMisiewicz left a comment

Choose a reason for hiding this comment

JonathonMisiewicz left a comment

Choose a reason for hiding this comment

jturney left a comment

Choose a reason for hiding this comment

maxscheurer left a comment

Choose a reason for hiding this comment

hokru commented Oct 26, 2021 •

edited

Loading

hokru Oct 28, 2021 •

edited

Loading