Use pivoted Cholesky in determining auxiliary basis for density fitting

src/density_fitting.cpp

@@ -52,7 +52,7 @@ bool DensityFit::Bmat_enabled() const {
   return Bmat;
 }
 
-size_t DensityFit::fill(const BasisSet & orbbas, const BasisSet & auxbas, bool dir, double erithr, double linthr, bool bmat) {
+size_t DensityFit::fill(const BasisSet & orbbas, const BasisSet & auxbas, bool dir, double erithr, double linthr, double cholthr, bool bmat) {
   // Construct density fitting basis
 
   // Store amount of functions
@@ -125,28 +125,8 @@ size_t DensityFit::fill(const BasisSet & orbbas, const BasisSet & auxbas, bool d
   }
 
   if(Bmat) {
-    // Form ab^-1 and ab^-1/2
-    arma::mat abvec;
-    arma::vec abval;
-    eig_sym_ordered(abval,abvec,ab);
-
-    // Count linearly independent vectors
-    size_t Nind=0;
-    for(size_t i=0;i<abval.n_elem;i++)
-      if(abval(i)>=linthr)
-	Nind++;
-
-    // and drop the linearly dependent ones
-    abval=abval.subvec(abval.n_elem-Nind,abval.n_elem-1);
-    abvec=abvec.cols(abvec.n_cols-Nind,abvec.n_cols-1);
-
-    // Form matrices
-    ab_inv.zeros(abvec.n_rows,abvec.n_rows);
-    ab_invh.zeros(abvec.n_rows,abvec.n_rows);
-    for(size_t i=0;i<abval.n_elem;i++) {
-      ab_inv+=abvec.col(i)*arma::trans(abvec.col(i))/abval(i);
-      ab_invh+=abvec.col(i)*arma::trans(abvec.col(i))/sqrt(abval(i));
-    }
+    ab_invh = PartialCholeskyOrth(ab, cholthr, linthr);
+    ab_inv = ab_invh * ab_invh.t();
   } 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));

src/density_fitting.h

@@ -135,7 +135,7 @@ class DensityFit {
    * 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, bool bmat=false);
+  size_t fill(const BasisSet & orbbas, const BasisSet & auxbas, bool direct, double erithr, double linthr, double cholthr, bool bmat=false);
 
   /// Compute estimate of necessary memory
   size_t memory_estimate(const BasisSet & orbbas, const BasisSet & auxbas, double erithr, bool direct) const;

src/localization.cpp

@@ -1203,11 +1203,11 @@ void Pipek::cost_func_der(const arma::cx_mat & Wv, double & Dinv, arma::cx_mat &
 Edmiston::Edmiston(const BasisSet & basis, const BasisSet & fitbas, const arma::mat & Cv, bool delocalize) : UnitaryFunction(4,!delocalize) {
   // Store orbitals
   C=Cv;
-  // Initialize fitting integrals. Direct computation, linear dependence threshold 1e-8, use Hartree-Fock routine since it has better tolerance for linear dependencies
+  // Initialize fitting integrals. Direct computation, linear dependence threshold 1e-8, Cholesky threshold 1e-9, use Hartree-Fock routine since it has better tolerance for linear dependencies
   if(!fitbas.get_Nbf())
-    dfit.fill(basis,basis.density_fitting(),true,1e-8,false);
+    dfit.fill(basis,basis.density_fitting(),true,1e-8,1e-9,false);
   else
-    dfit.fill(basis,fitbas,true,1e-8,false);
+    dfit.fill(basis,fitbas,true,1e-8,1e-9,false);
 
   use_chol=false;
 }

src/scf-base.cpp

@@ -144,6 +144,8 @@ SCF::SCF(const BasisSet & basis, Checkpoint & chkpt) {
   fitmem=1000000*settings.get_int("FittingMemory");
   // Linear dependence threshold
   fitthr=settings.get_double("FittingThreshold");
+  // Linear dependence threshold
+  fitcholthr=settings.get_double("FittingCholeskyThreshold");
 
   // Use Cholesky?
   cholesky=settings.get_bool("Cholesky");
@@ -347,7 +349,7 @@ SCF::SCF(const BasisSet & basis, Checkpoint & chkpt) {
     }
 
     // Calculate the fitting integrals, don't run in B-matrix mode
-    size_t Npairs=dfit.fill(*basisp,dfitbas,direct,intthr,fitthr,false);
+    size_t Npairs=dfit.fill(*basisp,dfitbas,direct,intthr,fitthr,fitcholthr,false);
 
     if(verbose) {
       printf("done (%s)\n",t.elapsed().c_str());
@@ -505,7 +507,7 @@ void SCF::fill_rs(double omega) {
       }
 
       t.set();
-      size_t Npairs=dfit_rs.fill(*basisp,dfitbas,direct,intthr,fitthr,dfit.Bmat_enabled());
+      size_t Npairs=dfit_rs.fill(*basisp,dfitbas,direct,intthr,fitthr,fitcholthr,dfit.Bmat_enabled());
       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());

src/scf.h

@@ -287,6 +287,8 @@ class SCF {
   size_t fitmem;
   /// Threshold for density fitting
   double fitthr;
+  /// Pivoted Cholesky threshold for density fitting
+  double fitcholthr;
 
   /// Cholesky calculation?
   bool cholesky;

src/settings.cpp

@@ -142,6 +142,7 @@ void Settings::add_scf_settings() {
   add_int("FittingMemory", "Amount of memory in MB to use for exchange fitting",1000);
   // Threshold for screening eigenvectors
   add_double("FittingThreshold", "Linear dependence threshold for Coulomb integrals in density fitting",1e-8);
+  add_double("FittingCholeskyThreshold", "Linear dependence threshold for pivoted Cholesky of Coulomb integrals in density fitting",1e-9);
 
   // SAP basis
   add_string("SAPBasis", "Tabulated atomic effective potential \"basis set\"","helfem_large.gbs");