From ec73829940e624c790c2a811ec8711337ae44624 Mon Sep 17 00:00:00 2001
From: Susi Lehtola <susi.lehtola@gmail.com>
Date: Sun, 21 Jan 2024 12:59:06 +0200
Subject: [PATCH] Remove method of Eichkorn et al from density fitting code,
 instead always use the numerically stable method of pivoted Cholesky
 orthogonalization.

---
 src/density_fitting.cpp | 76 ++++++++---------------------------------
 src/density_fitting.h   |  7 +---
 src/scf-base.cpp        |  6 ++--
 3 files changed, 19 insertions(+), 70 deletions(-)

diff --git a/src/density_fitting.cpp b/src/density_fitting.cpp
index ee6c00d6..168a0b3b 100644
--- a/src/density_fitting.cpp
+++ b/src/density_fitting.cpp
@@ -48,11 +48,7 @@ void DensityFit::get_range_separation(double & w, double & a, double & b) const
   b=beta;
 }
 
-bool DensityFit::Bmat_enabled() const {
-  return Bmat;
-}
-
-size_t DensityFit::fill(const BasisSet & orbbas, const BasisSet & auxbas, bool dir, double erithr, double linthr, double cholthr, bool bmat) {
+size_t DensityFit::fill(const BasisSet & orbbas, const BasisSet & auxbas, bool dir, double erithr, double linthr, double cholthr) {
   // Construct density fitting basis
 
   // Store amount of functions
@@ -60,7 +56,6 @@ size_t DensityFit::fill(const BasisSet & orbbas, const BasisSet & auxbas, bool d
   Naux=auxbas.get_Nbf();
   Nnuc=orbbas.get_Nnuc();
   direct=dir;
-  Bmat=bmat;
 
   // Fill list of shell pairs
   arma::mat Q, M;
@@ -124,15 +119,8 @@ size_t DensityFit::fill(const BasisSet & orbbas, const BasisSet & auxbas, bool d
     delete eri;
   }
 
-  if(Bmat) {
-    ab_invh = PartialCholeskyOrth(ab, cholthr, linthr);
-    ab_inv = ab_invh * ab_invh.t();
-    //printf("%i auxiliary funtions out of %i are linearly independent\n",ab_invh.n_cols,ab_invh.n_rows);
-  } else {
-    // Just RI-J(K), so use faster method from Eichkorn et al to form ab_inv only
-    ab_inv=arma::inv(ab + DELTA*arma::eye(ab.n_rows,ab.n_cols));
-    //printf("Using method of Eichkorn et al for ab_inv\n");
-  }
+  ab_invh = PartialCholeskyOrth(ab, cholthr, linthr);
+  ab_inv = ab_invh * ab_invh.t();
 
   // Then, compute the three-center integrals
   if(!direct) {
@@ -384,12 +372,8 @@ void DensityFit::digest_K_incore(const arma::mat & C, const arma::vec & occs, ar
     }
 
     // K_uv = (ui|vi) = (a|ui) (a|b)^-1 (b|vi)
-    if(Bmat) {
-      aui = ab_invh.t()*aui;
-      K += occs[io]*arma::trans(aui)*aui;
-    } else {
-      K += occs[io]*arma::trans(aui)*ab_inv*aui;
-    }
+    aui = ab_invh.t()*aui;
+    K += occs[io]*arma::trans(aui)*aui;
   }
 }
 
@@ -465,12 +449,8 @@ void DensityFit::digest_K_incore(const arma::cx_mat & C, const arma::vec & occs,
     }
 
     // K_uv = (ui|vi) = (a|ui) (a|b)^-1 (b|vi)
-    if(Bmat) {
-      aui = ab_invh.t()*aui;
-      K += occs[io]*arma::trans(aui)*aui;
-    } else {
-      K += occs[io]*arma::trans(aui)*ab_inv*aui;
-    }
+    aui = ab_invh.t()*aui;
+    K += occs[io]*arma::trans(aui)*aui;
   }
 }
 
@@ -565,15 +545,9 @@ void DensityFit::digest_K_direct(const arma::mat & C, const arma::vec & occs, ar
   }
 
   // K_uv = (ui|vi) = (a|ui) (a|b)^-1 (b|vi)
-  if(Bmat) {
-    for(size_t io=0;io<C.n_cols;io++) {
-      aui[io] = ab_invh*aui[io];
-      K += occs[io]*arma::trans(aui[io])*aui[io];
-    }
-  } else {
-    for(size_t io=0;io<C.n_cols;io++) {
-      K += occs[io]*arma::trans(aui[io])*ab_inv*aui[io];
-    }
+  for(size_t io=0;io<C.n_cols;io++) {
+    aui[io] = ab_invh*aui[io];
+    K += occs[io]*arma::trans(aui[io])*aui[io];
   }
 }
 
@@ -598,9 +572,8 @@ size_t DensityFit::memory_estimate(const BasisSet & orbbas, const BasisSet & aux
 
   // Memory taken by (\alpha | \beta) and its inverse
   Nmem+=2*Na*Na*sizeof(double);
-  if(Bmat)
-    // We also have (a|b)^(-1/2)
-    Nmem+=Na*Na*sizeof(double);
+  // We also have (a|b)^(-1/2)
+  Nmem+=Na*Na*sizeof(double);
 
   // Memory taken by gamma and expansion coefficients
   Nmem+=2*Na*sizeof(double);
@@ -675,14 +648,7 @@ arma::vec DensityFit::compute_expansion(const arma::mat & P) const {
   }
 
   // Compute and return c
-  if(Bmat) {
-    return ab_inv*gamma;
-  } else {
-    // Compute x0
-    arma::vec x0=ab_inv*gamma;
-    // Compute and return c
-    return x0+ab_inv*(gamma-ab*x0);
-  }
+  return ab_inv*gamma;
 }
 
 std::vector<arma::vec> DensityFit::compute_expansion(const std::vector<arma::mat> & P) const {
@@ -760,18 +726,8 @@ std::vector<arma::vec> DensityFit::compute_expansion(const std::vector<arma::mat
   }
 
   // Compute and return c
-  if(Bmat) {
-    for(size_t ig=0;ig<P.size();ig++)
-      gamma[ig]=ab_inv*gamma[ig];
-
-  } else {
-    for(size_t ig=0;ig<P.size();ig++) {
-      // Compute x0
-      arma::vec x0=ab_inv*gamma[ig];
-      // Compute and return c
-      gamma[ig]=x0+ab_inv*(gamma[ig]-ab*x0);
-    }
-  }
+  for(size_t ig=0;ig<P.size();ig++)
+    gamma[ig]=ab_inv*gamma[ig];
 
   return gamma;
 }
@@ -1238,8 +1194,6 @@ void DensityFit::three_center_integrals(arma::mat & ints) const {
 void DensityFit::B_matrix(arma::mat & B) const {
   if(direct)
     throw std::runtime_error("Must run in tabulated mode!\n");
-  if(!Bmat)
-    throw std::runtime_error("Must be run in B-matrix mode!\n");
   // Compute the integrals
   three_center_integrals(B);
   // Transform into proper B matrix
diff --git a/src/density_fitting.h b/src/density_fitting.h
index 341995f9..76588487 100644
--- a/src/density_fitting.h
+++ b/src/density_fitting.h
@@ -61,8 +61,6 @@ class DensityFit {
   size_t Naux;
   /// Direct calculation? (Compute three-center integrals on-the-fly)
   bool direct;
-  /// B-matrix calculation?
-  bool Bmat;
 
   /// Range separation constants
   double omega, alpha, beta;
@@ -126,16 +124,13 @@ class DensityFit {
   /// Get range separation constants
   void get_range_separation(double & w, double & a, double & b) const;
 
-  /// Running in B-matrix mode?
-  bool Bmat_enabled() const;
-
   /**
    * Compute integrals, use given linear dependency threshold. The HF
    * flag here controls formation of (a|b)^{-1/2} and (a|b)^{-1}; the
    * HF routine should be more tolerant of linear dependencies in the basis.
    * Returns amount of significant orbital shell pairs.
    */
-  size_t fill(const BasisSet & orbbas, const BasisSet & auxbas, bool direct, double erithr, double linthr, double cholthr, bool bmat=false);
+  size_t fill(const BasisSet & orbbas, const BasisSet & auxbas, bool direct, double erithr, double linthr, double cholthr);
 
   /// Compute estimate of necessary memory
   size_t memory_estimate(const BasisSet & orbbas, const BasisSet & auxbas, double erithr, bool direct) const;
diff --git a/src/scf-base.cpp b/src/scf-base.cpp
index f667f542..92bab146 100644
--- a/src/scf-base.cpp
+++ b/src/scf-base.cpp
@@ -348,8 +348,8 @@ SCF::SCF(const BasisSet & basis, Checkpoint & chkpt) {
       t.set();
     }
 
-    // Calculate the fitting integrals, don't run in B-matrix mode
-    size_t Npairs=dfit.fill(*basisp,dfitbas,direct,intthr,fitthr,fitcholthr,false);
+    // Calculate the fitting integrals
+    size_t Npairs=dfit.fill(*basisp,dfitbas,direct,intthr,fitthr,fitcholthr);
 
     if(verbose) {
       printf("done (%s)\n",t.elapsed().c_str());
@@ -507,7 +507,7 @@ void SCF::fill_rs(double omega) {
       }
 
       t.set();
-      size_t Npairs=dfit_rs.fill(*basisp,dfitbas,direct,intthr,fitthr,fitcholthr,dfit.Bmat_enabled());
+      size_t Npairs=dfit_rs.fill(*basisp,dfitbas,direct,intthr,fitthr,fitcholthr);
       if(verbose) {
 	printf("done (%s)\n",t.elapsed().c_str());
 	printf("%i shell pairs out of %i are significant.\n",(int) Npairs, (int) basisp->get_unique_shellpairs().size());