[Draft] DLPNO-CCSD PR #2915

andyj10224 · 2023-03-30T15:49:01Z

Description

HERE IT IS!!! This is a draft of the DLPNO-CCSD PR that will be coming in the next few months. The purpose of this is for the developers and research groups to be able to run and test DLPNO-CCSD before it is officially part of the code.

Credit to @JoseMadriaga for the derivations
LocalCCSD1to10.pdf

Useful References:
Original DLPNO-CCSD Paper
Sparse Maps II Paper

Example Input File

memory 20 GB

molecule mol {
  0 1
  O    0.705    0.744    0.16
  H    -0.071    0.264    0.45
  H    1.356    0.064    -0.014
  symmetry c1
}

set {
  basis cc-pVDZ
  scf_type df
  freeze_core true
  pno_convergence normal
}
energy('dlpno-ccsd')

Results (Waterclusters in TZ)

[Speedups, relative to DF-CCSD]

[Percent Correlation Energy Recovered, relative to DF-CCSD, all >= 99.9%]

User API & Changelog headlines

Implement the DLPNO-CCSD algorithm

Dev notes & details

Feel free to use this code, it is not fully tested yet, but preliminary tests show encouraging results, and is MUCH faster than conventional CCSD
If you benchmark my code, please post results in the thread

Questions

Question1

Checklist

Add documentation
Add references to terms and equations
Tests added for any new features
All or relevant fraction of full tests run

Status

Ready for review
Ready for merge

Co-authored-by: TiborGY <tibor.gyori@chem.u-szeged.hu> Co-authored-by: Lori A. Burns <lori.burns@gmail.com>

andyj10224 · 2023-04-01T23:24:44Z

Just implemented some of Lori and Tibor's suggestions. I have also implemented SC-LMP2 for "weak pairs" to reduce the cost of the LCCSD computation, per the Sparse Maps II paper.

TiborGY

Round 2, mostly more of the same: const correctness, a few unused variables, a docstring suggestion, and a possible suggestion to eliminate a bit of code duplication,

TiborGY · 2023-04-04T12:29:11Z

psi4/src/psi4/dlpno/dlpnobase.cc

+/* Utility function for making C_DGESV calls
+ *
+ * C_DGESV solves AX=B for X, given symmetric NxN matrix A and NxM matrix B
+ * B is expected in fortran layout, which complicates the call when (M > 1)
+ * The workaround used here is to switch the layout of B before and after the call
+ */


Suggested change

/* Utility function for making C_DGESV calls

*

* C_DGESV solves AX=B for X, given symmetric NxN matrix A and NxM matrix B

* B is expected in fortran layout, which complicates the call when (M > 1)

* The workaround used here is to switch the layout of B before and after the call

*/

Moving this to the .h should make it more digestible for IDEs and doc generators.

TiborGY · 2023-04-04T12:32:09Z

psi4/src/psi4/dlpno/dlpno.h

+
+      void common_init();
+
+      // Helper functions


Suggested change

// Helper functions

// Helper functions

/// @brief Utility function for making C_DGESV calls. C_DGESV solves AX=B for X, given symmetric NxN matrix A and NxM

/// matrix B, but B is expected in fortran layout, which complicates the call when (M > 1).

/// The workaround used here is to switch the layout of B before and after the call.

///

/// @param A On entry, the n-by-n coefficient matrix A. On exit, the factors L and U from the factorization A = P*L*U;

/// the unit diagonal elements of L are not stored.

/// @param B On entry, the n-by-m matrix of right hand side matrix B. On exit, the n-by-m solution matrix X.

/// \ingroup DLPNO

Moved here from the .cc, rewritten for machine-readability.

TiborGY · 2023-04-04T12:46:04Z

psi4/src/psi4/dlpno/dlpnobase.cc

+                                     0);
+    auto X = orthog.basis_to_orthog_basis();
+
+    int nmo_initial = X->rowspi(0);


Suggested change

int nmo_initial = X->rowspi(0);

Unused variable

TiborGY · 2023-04-04T12:59:21Z

psi4/src/psi4/dlpno/dlpnobase.cc

+    }
+
+    auto flat = std::make_shared<Vector>("flattened matrix list", total_size);
+    double* flatp = flat->pointer();


Suggested change

double* flatp = flat->pointer();

double* const flatp = flat->pointer();

The pointer itself can be const

TiborGY · 2023-04-04T13:03:59Z

psi4/src/psi4/dlpno/dlpnobase.cc

+void DLPNOBase::copy_flat_mats(SharedVector flat, std::vector<SharedMatrix>& mat_list) {
+    double* flatp = flat->pointer();


Suggested change

void DLPNOBase::copy_flat_mats(SharedVector flat, std::vector<SharedMatrix>& mat_list) {

double* flatp = flat->pointer();

void DLPNOBase::copy_flat_mats(const SharedVector flat, std::vector<SharedMatrix>& mat_list) {

const double* const flatp = flat->pointer();

TiborGY · 2023-04-04T15:11:54Z

psi4/src/psi4/dlpno/dlpnobase.cc

+            int centerU = basisset_->function_to_center(u);
+            double p_uu = P_i->get(u, u);


Suggested change

int centerU = basisset_->function_to_center(u);

double p_uu = P_i->get(u, u);

const int centerU = basisset_->function_to_center(u);

const double p_uu = P_i->get(u, u);

TiborGY · 2023-04-04T15:12:08Z

psi4/src/psi4/dlpno/dlpnobase.cc

+                int centerV = basisset_->function_to_center(v);
+                double p_vv = P_i->get(v, v);


Suggested change

int centerV = basisset_->function_to_center(v);

double p_vv = P_i->get(v, v);

const int centerV = basisset_->function_to_center(v);

const double p_vv = P_i->get(v, v);

TiborGY · 2023-04-04T15:13:01Z

psi4/src/psi4/dlpno/dlpnobase.cc

+                double p_vv = P_i->get(v, v);
+
+                // off-diag pops (p_uv) split between u and v prop to diag pops
+                double p_uv = P_i->get(u, v);


Suggested change

double p_uv = P_i->get(u, v);

const double p_uv = P_i->get(u, v);

TiborGY · 2023-04-04T15:45:58Z

psi4/src/psi4/dlpno/dlpnobase.cc

+    if (options_.get_str("DLPNO_ALGORITHM") == "MP2") {
+        for (size_t i = 0, ij = 0; i < naocc; i++) {
+            for (size_t j = 0; j < naocc; j++) {
+                bool overlap_big = (DOI_ij_->get(i, j) > options_.get_double("T_CUT_DO_ij"));
+                bool energy_big = (fabs(dipole_pair_e_bound_->get(i, j)) > options_.get_double("T_CUT_PRE"));
+
+                if (overlap_big || energy_big) {
+                    i_j_to_ij_[i].push_back(ij);
+                    ij_to_i_j_.push_back(std::make_pair(i, j));
+                    ij++;
+                } else {
+                    de_dipole_ += dipole_pair_e_->get(i, j);
+                    i_j_to_ij_[i].push_back(-1);
+                }
+            }
+        }
+    } else {
+        for (size_t i = 0, ij = 0; i < naocc; i++) {
+            for (size_t j = 0; j < naocc; j++) {
+                bool overlap_big = (DOI_ij_->get(i, j) > options_.get_double("T_CUT_DO_ij"));
+                bool energy_big = (fabs(dipole_pair_e_bound_->get(i, j)) > options_.get_double("T_CUT_PRE"));
+
+                if ((i == j) || (overlap_big && energy_big)) {
+                    i_j_to_ij_[i].push_back(ij);
+                    ij_to_i_j_.push_back(std::make_pair(i, j));
+                    ij++;
+                } else {
+                    if (overlap_big || energy_big)
+                        weak_pairs_.push_back(std::make_pair(i,j));
+                    else
+                        de_dipole_ += dipole_pair_e_->get(i, j);
+
+                    i_j_to_ij_[i].push_back(-1);
+                }
+            }
+        }
+    }


Suggested change

if (options_.get_str("DLPNO_ALGORITHM") == "MP2") {

for (size_t i = 0, ij = 0; i < naocc; i++) {

for (size_t j = 0; j < naocc; j++) {

bool overlap_big = (DOI_ij_->get(i, j) > options_.get_double("T_CUT_DO_ij"));

bool energy_big = (fabs(dipole_pair_e_bound_->get(i, j)) > options_.get_double("T_CUT_PRE"));

if (overlap_big || energy_big) {

i_j_to_ij_[i].push_back(ij);

ij_to_i_j_.push_back(std::make_pair(i, j));

ij++;

} else {

de_dipole_ += dipole_pair_e_->get(i, j);

i_j_to_ij_[i].push_back(-1);

}

}

}

} else {

for (size_t i = 0, ij = 0; i < naocc; i++) {

for (size_t j = 0; j < naocc; j++) {

bool overlap_big = (DOI_ij_->get(i, j) > options_.get_double("T_CUT_DO_ij"));

bool energy_big = (fabs(dipole_pair_e_bound_->get(i, j)) > options_.get_double("T_CUT_PRE"));

if ((i == j) || (overlap_big && energy_big)) {

i_j_to_ij_[i].push_back(ij);

ij_to_i_j_.push_back(std::make_pair(i, j));

ij++;

} else {

if (overlap_big || energy_big)

weak_pairs_.push_back(std::make_pair(i,j));

else

de_dipole_ += dipole_pair_e_->get(i, j);

i_j_to_ij_[i].push_back(-1);

}

}

}

}

const bool DLPNO_MP2 = (options_.get_str("DLPNO_ALGORITHM") == "MP2");

for (size_t i = 0, ij = 0; i < naocc; i++) {

for (size_t j = 0; j < naocc; j++) {

const bool overlap_big = (DOI_ij_->get(i, j) > options_.get_double("T_CUT_DO_ij"));

const bool energy_big = (fabs(dipole_pair_e_bound_->get(i, j)) > options_.get_double("T_CUT_PRE"));

if(DLPNO_MP2){

if (overlap_big || energy_big) {

i_j_to_ij_[i].push_back(ij);

ij_to_i_j_.push_back(std::make_pair(i, j));

ij++;

} else {

de_dipole_ += dipole_pair_e_->get(i, j);

i_j_to_ij_[i].push_back(-1);

}

}else{

if ((i == j) || (overlap_big && energy_big)) {

i_j_to_ij_[i].push_back(ij);

ij_to_i_j_.push_back(std::make_pair(i, j));

ij++;

} else {

if (overlap_big || energy_big)

weak_pairs_.push_back(std::make_pair(i, j));

else

de_dipole_ += dipole_pair_e_->get(i, j);

i_j_to_ij_[i].push_back(-1);

}

}

}

}

This would eliminate some code duplication. Thoughts?

TiborGY · 2023-04-04T15:48:09Z

psi4/src/psi4/dlpno/dlpnobase.cc

+        }
+    }
+
+    int n_lmo_pairs = ij_to_i_j_.size();


Suggested change

int n_lmo_pairs = ij_to_i_j_.size();

const int n_lmo_pairs = ij_to_i_j_.size();

davpoolechem · 2023-04-06T14:35:34Z

One thing we will need to discuss up-front, is what to do with a DLPNO-CC Psi4NumPy implementation. I know @JonathonMisiewicz talked about this in a previous PR, and I see his point. It would be useful to have an easier-to-understand reference for DLPNO-CC somewhere, especially to help others understand the DLPNO formalism.

A lot of the work done in this PR was previously prototyped in Python, so that gives us a start in a potential Psi4NumPy implementation.

andyj10224 added 30 commits December 22, 2022 11:27

Fix BasisFunctions constructor

acacf2c

Expose DLPNO wfn info to Python

c3d690e

add Qij and Qab ints

11b2262

Add getter for qia

5eae125

Implement C_pao slices in DLPNO-MP2

430b480

add K_maef ints C side

c3681d5

Add K_abef integrals

a6eb8a2

Implement K_mbij integrals

8ce3d2d

Implement local K_mnij integrals + update K_mbij

f85c4b5

Implement J_ijab integrals

4fb686f

Implement machinery for DLPNO-(T)

018ae14

Working DLPNO-CCSD in C++ :)

04d97fb

Implemented DIIS

086f282

Remove unnecessary copies from lccsd iterations

5605d5a

Format output

2a18a41

Fix segfault when n_pno == 0

6f065c9

Optimize J_ijab

22136a7

Optimize Fme intermediate

23032bc

Optimize Fbe term

ec31e6e

Optimize memory of R2 update

6da24f3

Get rid of forming K_abef explicitly

c1bc5dc

Remove explicit three virtual integrals

ce08259

Memory optimizations

2c357a2

Remove storage of Qab_ij

7f2c9c6

Make Wmnij more efficient

05474b3

More efficient integrals

2f5448b

Add memory bounds

f7d1a35

Clean up memory printing

7826f12

Clean up four virtual algorithm control

4ff83e8

Add med memory algo

d4f0347

andyj10224 and others added 2 commits April 1, 2023 18:45

Apply suggestions from code review (loriab and TiborGY)

a6592e3

Co-authored-by: TiborGY <tibor.gyori@chem.u-szeged.hu> Co-authored-by: Lori A. Burns <lori.burns@gmail.com>

Incorporate more of TiborGY suggestions

aa8a310

Segfault free SC-LMP2 weak pair algo

8cc1315

andyj10224 force-pushed the dlpno-ccsd branch from a60ac55 to 8cc1315 Compare April 3, 2023 17:49

Make SC-LMP2 more friendly towards icpc

704a420

TiborGY reviewed Apr 4, 2023

View reviewed changes

andyj10224 added 3 commits April 14, 2023 11:21

Refactor Integrals

5ff1cd2

Update integral printouts

1ca5f6f

Fix L_tilde bug

330f277

loriab added this to the Psi4 1.9 milestone Apr 18, 2023

loriab added feature Extends an existing Psi feature or develops a new one. cc For all issues involving the CC module, ground-state energies to response properties. labels Apr 18, 2023

andyj10224 added 14 commits May 16, 2023 10:42

Update definition of strong and weak pairs

3e09b4c

Prettify Print Statements

9a751b3

Update pno_convergence def based on Liakos 2015

8daea85

Smarter integral generation and prescreening

fc66ef4

Add crude pair prescreening

10cf820

Working SVD Qab int factorization

ca9b810

Clear Qab SVD after its use is done

906e6ab

Add T_CUT_EIG option

a54c4e8

Fix race condition on estimating pno overlap mem

a2ba61a

Reset memory of qab

1e4fba3

Optimize CC ints performance

1c39cf2

Cleanups with grid and printing

a59271c

Clean up reshape

bbff8f1

Fix bug in not transposing

1bec86a

loriab modified the milestones: Psi4 1.9, Psi4 1.10 Nov 13, 2023

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[Draft] DLPNO-CCSD PR #2915

[Draft] DLPNO-CCSD PR #2915

andyj10224 commented Mar 30, 2023 •

edited

andyj10224 commented Apr 1, 2023

TiborGY left a comment

TiborGY Apr 4, 2023

TiborGY Apr 4, 2023

TiborGY Apr 4, 2023

TiborGY Apr 4, 2023

TiborGY Apr 4, 2023

TiborGY Apr 4, 2023

TiborGY Apr 4, 2023

TiborGY Apr 4, 2023

TiborGY Apr 4, 2023

TiborGY Apr 4, 2023

davpoolechem commented Apr 6, 2023

-      // Helper functions
+      // Helper functions
+      /// @brief Utility function for making C_DGESV calls. C_DGESV solves AX=B for X, given symmetric NxN matrix A and NxM
+      /// matrix B, but B is expected in fortran layout, which complicates the call when (M > 1).
+      /// The workaround used here is to switch the layout of B before and after the call.
+      ///
+      /// @param A On entry, the n-by-n coefficient matrix A. On exit, the factors L and U from the factorization A = P*L*U;
+      /// the unit diagonal elements of L are not stored.
+      /// @param B On entry, the n-by-m matrix of right hand side matrix B.  On exit, the n-by-m solution matrix X.
+      /// \ingroup DLPNO

	double* flatp = flat->pointer();
	double* const flatp = flat->pointer();

		void DLPNOBase::copy_flat_mats(SharedVector flat, std::vector<SharedMatrix>& mat_list) {
		double* flatp = flat->pointer();

		int centerU = basisset_->function_to_center(u);
		double p_uu = P_i->get(u, u);

		int centerV = basisset_->function_to_center(v);
		double p_vv = P_i->get(v, v);

	double p_uv = P_i->get(u, v);
	const double p_uv = P_i->get(u, v);

	int n_lmo_pairs = ij_to_i_j_.size();
	const int n_lmo_pairs = ij_to_i_j_.size();

[Draft] DLPNO-CCSD PR #2915

Are you sure you want to change the base?

[Draft] DLPNO-CCSD PR #2915

Conversation

andyj10224 commented Mar 30, 2023 • edited

Description

Results (Waterclusters in TZ)

User API & Changelog headlines

Dev notes & details

Questions

Checklist

Status

andyj10224 commented Apr 1, 2023

TiborGY 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

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

davpoolechem commented Apr 6, 2023

andyj10224 commented Mar 30, 2023 •

edited