cadfclhf.h

//
// cadfclhf.h
//
// Copyright (C) 2013 David Hollman
//
// Author: David Hollman
// Maintainer: DSH
// Created: Nov 13, 2013
//
// This file is part of the SC Toolkit.
//
// The SC Toolkit is free software; you can redistribute it and/or modify
// it under the terms of the GNU Library General Public License as published by
// the Free Software Foundation; either version 2, or (at your option)
// any later version.
//
// The SC Toolkit is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
//
// The U.S. Government is granted a limited license as per AL 91-7.
//
 
 
#ifndef _chemistry_qc_scf_cadfclhf_h
#define _chemistry_qc_scf_cadfclhf_h
 
#define USE_SPARSE 0
 
// Standard library includes
#include <atomic>
#include <future>
#include <functional>
#include <type_traits>
#include <unordered_set>
#include <unordered_map>
 
// Eigen includes
#include <Eigen/Dense>
#include <Eigen/Sparse>
 
// MPQC includes
#include <chemistry/qc/scf/clhf.h>
#include <chemistry/qc/lcao/wfnworld.h>
#include <util/container/conc_cache_fwd.h>
#include <util/container/thread_wrap.h>
#include <util/misc/hash.h>
 
#include "iters.h"
#include "ordered_shells.h"
#include "treemat_fwd.h"
 
// Allows easy switching between std::shared_ptr and e.g. boost::shared_ptr
template<typename... Types>
using shared_ptr = std::shared_ptr<Types...>;
using std::make_shared;
 
typedef typename boost::integer_range<int> int_range;
typedef unsigned long long ull;
 
#define USE_INTEGRAL_CACHE 0
 
namespace sc {
 
// Forward Declarations
 
class XMLWriter;
class ApproximatePairWriter;
namespace cadf {
  class AssignmentGrid;
  namespace assignments {
    struct Assignments;
  }
}
 
//============================================================================//
 
inline void
do_threaded(int nthread, const std::function<void(int)>& f){
  boost::thread_group compute_threads;
  // Loop over number of threads
  for(int ithr = 0; ithr < nthread; ++ithr) {
    // create each thread that runs f
    compute_threads.create_thread([&, ithr](){
      // run the work
      f(ithr);
    });
  }
  // join the created threads
  compute_threads.join_all();
}
 
namespace cadf {
 
class Histogram2d {
  public:
    typedef Eigen::Matrix<uli, Eigen::Dynamic, Eigen::Dynamic> matrix_t;
 
  private:
 
    matrix_t hist_mat_;
    double interval_h_, interval_v_;
    bool log_h_, log_v_;
    int nbins_h_, nbins_v_;
    bool clip_edges_;
 
  public:
 
    double min_h_, max_h_, min_v_, max_v_;
 
    Histogram2d(
        int nbins_h, double min_h, double max_h,
        int nbins_v, double min_v, double max_v,
        bool log_h=false, bool log_v=false,
        bool clip_edges=false
    );
 
    const matrix_t& matrix() const { return hist_mat_; }
 
    void insert_value(double value_h, double value_v);
 
    void accumulate(const Histogram2d& other) {
      hist_mat_ += other.hist_mat_;
    }
 
};
 
 
}
 
//============================================================================//
//============================================================================//
//============================================================================//
 
 
class CADFCLHF: public CLHF {
 
  protected:
 
    void ao_fock(double accuracy);
    void reset_density();
    double new_density();
    void init_vector();
 
  public:
 
    // typedefs
 
    typedef Eigen::Map<Eigen::VectorXd, Eigen::Unaligned, Eigen::Stride<1, Eigen::Dynamic>> CoefView;
    typedef shared_ptr<CoefView> CoefContainer;
    typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> TwoCenterIntContainer;
    typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> ThreeCenterIntContainer;
    typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> FourCenterIntContainer;
    typedef shared_ptr<TwoCenterIntContainer> TwoCenterIntContainerPtr;
    typedef shared_ptr<ThreeCenterIntContainer> ThreeCenterIntContainerPtr;
    typedef shared_ptr<FourCenterIntContainer> FourCenterIntContainerPtr;
    typedef std::unordered_map<
        std::pair<int, int>,
        std::pair<CoefContainer, CoefContainer>,
        sc::hash<std::pair<int, int>>
    > CoefMap;
    typedef Eigen::HouseholderQR<Eigen::MatrixXd> Decomposition;
    typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> RowMatrix;
    typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> ColMatrix;
    typedef Eigen::Map<RowMatrix, Eigen::Unaligned, Eigen::OuterStride<>> StridedRowMap;
    template<typename T1, typename T2>
    using fast_map = std::unordered_map<T1, T2, sc::hash<T1>>;
 
    // TODO Get rid of this
    typedef ConcurrentCache<
        double,
        int, int, int, int, TwoBodyOper::type
    > FourCenterMaxIntCache;
    typedef ConcurrentCache<shared_ptr<Decomposition>, int, int> DecompositionCache;
 
    CADFCLHF(StateIn&);
 
 
    CADFCLHF(const Ref<KeyVal>&);
 
    ~CADFCLHF();
 
    void save_data_state(StateOut&);
 
    class ScreeningStatistics {
 
      public:
 
        typedef std::atomic_uint_fast64_t accumulate_t;
        typedef decltype(accumulate_t().load()) count_t;
 
        bool histogram_mode = false;
 
        struct Iteration {
 
          public:
 
            Iteration() = default;
 
            Iteration(Iteration&& other)
              : int_screening_values(other.int_screening_values),
                int_actual_values(other.int_actual_values),
                int_distance_factors(other.int_distance_factors),
                int_distances(other.int_distances),
                int_indices(other.int_indices),
                int_ams(other.int_ams)
            {
              K_3c_needed.store(other.K_3c_needed.load());
              K_3c_needed_fxn.store(other.K_3c_needed_fxn.load());
              K_3c_dist_screened.store(other.K_3c_dist_screened.load());
              K_3c_dist_screened_fxn.store(other.K_3c_dist_screened_fxn.load());
              K_3c_underestimated.store(other.K_3c_underestimated.load());
              K_3c_underestimated_fxn.store(other.K_3c_underestimated_fxn.load());
              K_3c_perfect.store(other.K_3c_perfect.load());
              K_3c_perfect_fxn.store(other.K_3c_perfect_fxn.load());
              K_3c_contract_fxn.store(other.K_3c_contract_fxn.load());
              K_2c_contract_fxn.store(other.K_2c_contract_fxn.load());
              Kt_contract1_fxn.store(other.Kt_contract1_fxn.load());
              Kt_contract2_fxn.store(other.Kt_contract2_fxn.load());
              L3_build_compares.store(other.L3_build_compares.load());
 
              parent = other.parent;
              is_first = other.is_first;
            }
 
            accumulate_t K_3c_needed = { 0 };
            accumulate_t K_3c_needed_fxn = { 0 };
            accumulate_t K_3c_dist_screened = { 0 };
            accumulate_t K_3c_dist_screened_fxn = { 0 };
            accumulate_t K_3c_underestimated = { 0 };
            accumulate_t K_3c_underestimated_fxn = { 0 };
            accumulate_t K_3c_perfect = { 0 };
            accumulate_t K_3c_perfect_fxn = { 0 };
            accumulate_t K_3c_contract_fxn = { 0 };
            accumulate_t K_2c_contract_fxn = { 0 };
            accumulate_t Kt_contract1_fxn = { 0 };
            accumulate_t Kt_contract2_fxn = { 0 };
            accumulate_t L3_build_compares = { 0 };
 
            ScreeningStatistics* parent;
 
            ThreadReplicated<std::vector<double>> int_screening_values;
            ThreadReplicated<std::vector<double>> int_actual_values;
            ThreadReplicated<std::vector<double>> int_distance_factors;
            ThreadReplicated<std::vector<double>> int_distances;
            ThreadReplicated<std::vector<std::tuple<int, int, int>>> int_indices;
            ThreadReplicated<std::vector<std::tuple<int, int, int>>> int_ams;
            std::vector<cadf::Histogram2d> values_hists;
            std::vector<cadf::Histogram2d> distance_hists;
            std::vector<cadf::Histogram2d> distance_noschwarz_hists;
            std::vector<cadf::Histogram2d> exponent_ratio_hists;
            std::vector<ull> int_am_counts;
            std::vector<double> int_am_ratio_sums;
            std::vector<double> int_am_ratio_log_sums;
 
            void set_nthread(int nthr) {
              int_screening_values.set_nthread(nthr);
              int_actual_values.set_nthread(nthr);
              int_distance_factors.set_nthread(nthr);
              int_distances.set_nthread(nthr);
              int_indices.set_nthread(nthr);
              int_ams.set_nthread(nthr);
            }
 
            void global_sum(const Ref<MessageGrp>& msg) {
              auto accum_all = [&msg](accumulate_t& val) {
                //decltype(val.load()) tmp = val.load();
                long tmp = (long)val.load();
                msg->sum(&tmp, 1);
                val.store(tmp);
              };
 
              accum_all(K_3c_needed);
              accum_all(K_3c_needed_fxn);
              accum_all(K_3c_dist_screened);
              accum_all(K_3c_dist_screened_fxn);
              accum_all(K_3c_underestimated);
              accum_all(K_3c_underestimated_fxn);
              accum_all(K_3c_perfect);
              accum_all(K_3c_perfect_fxn);
              accum_all(K_3c_contract_fxn);
              accum_all(K_2c_contract_fxn);
              accum_all(Kt_contract1_fxn);
              accum_all(Kt_contract2_fxn);
              accum_all(L3_build_compares);
            }
 
            bool is_first = false;
 
        };
 
        mutable accumulate_t sig_pairs = { 0 };
        mutable accumulate_t sig_pairs_fxn = { 0 };
        int print_level = 0;
        mutable bool xml_stats_saved = false;
 
        void global_sum(const Ref<MessageGrp>& msg) const {
          auto accum_all = [&msg](accumulate_t& val) {
            long tmp = val.load();
            msg->sum(&tmp, 1);
            val.store(tmp);
          };
          accum_all(sig_pairs);
          accum_all(sig_pairs_fxn);
        }
 
        std::vector<Iteration> iterations;
 
        ScreeningStatistics()
          : iterations(),
            print_level(0)
        { }
 
        Iteration& next_iteration() {
          iterations.emplace_back();
          iterations.back().is_first = iterations.size() == 1;
          iterations.back().parent = this;
          return iterations.back();
        }
 
        void print_summary(std::ostream& out,
            const Ref<GaussianBasisSet>& basis,
            const Ref<GaussianBasisSet>& dfbs,
            const CADFCLHF* parent,
            int print_level = -1,
            bool new_exchange = false
        ) const;
    };
 
  private:
 
    typedef enum {
      AllPairs = 0,
      SignificantPairs = 1,
    } PairSet;
 
    int nthread_;
    double screening_thresh_ = 1e-7;
    double pair_screening_thresh_;
    double full_screening_thresh_;
    double coef_screening_thresh_;
    double distance_screening_thresh_;
    bool do_linK_ = true;
    bool linK_block_rho_;
    bool linK_use_distance_ = true;
    int print_screening_stats_;
    double distance_damping_factor_ = 1.0;
    bool print_iteration_timings_ = false;
    bool use_norms_nu_ = true;
    bool use_norms_sigma_ = false;
    bool sigma_norms_chunk_by_atoms_ = false;
    bool xml_debug_ = false;
    bool xml_screening_data_ = false;
    bool use_extents_ = true;
    bool use_max_extents_ = true;
    bool subtract_extents_ = false;
;
    double well_separated_thresh_ = 1e-1;
    bool exact_diagonal_J_ = false;
    bool exact_diagonal_K_ = false;
    bool cadf_K_only_ = false;
    bool B_use_buffer_ = false;
    size_t B_buffer_size_;
    // TODO individual options to scale other screening thresholds
    bool scale_screening_thresh_ = true;
    // TODO individual screening threshold minima
    double full_screening_thresh_min_;
    double full_screening_expon_ = 1.0;
    bool screen_B_ = true;
    bool distribute_coefficients_ = false;
    double screen_B_use_distance_ = true;
    double B_screening_thresh_;
    bool store_coefs_transpose_ = false;
    bool thread_4c_ints_ = false;
    double d_over_screening_thresh_;
    double d_under_screening_thresh_;
    bool use_norms_B_ = false;
    bool shuffle_L_3_keys_ = true;
    bool shuffle_J_assignments_ = true;
    bool match_orbitals_ = true;
    double match_orbitals_max_diff_ = 1e-4;
    double match_orbitals_max_homo_offset_ = 0.05;
    bool match_orbitals_use_svd_ = true;
    bool debug_coulomb_energy_ = false;
    bool debug_exchange_energy_ = false;
    bool new_exchange_algorithm_ = true;
    int min_atoms_per_node_ = 3;
    bool count_ints_only_ = false;
    bool count_ints_use_norms_ = true;
    double count_ints_exclude_thresh_ = 1e-16;
    bool count_ints_histogram_ = false;
    int count_ints_hist_distance_bins_ = 25;
    int count_ints_hist_ratio_bins_ = 300;
    int count_ints_hist_bins_ = 200;
    double count_ints_hist_min_ratio_ = 1e-8;
    double count_ints_hist_max_ratio_ = 1e8;
    double count_ints_hist_max_int_ = 1e3;
    double count_ints_max_distance_ = -1.0;
    bool count_ints_hist_clip_edges_ = false;
    bool count_ints_exclude_two_center_ = false;
    bool qqr_only_ = false;
    bool safe_extents_ = false;
    bool dist_factor_use_overlap_ = false;
    int dist_factor_overlap_max_am_diff_ = 99;
    bool dist_factor_overlap_exclude_contracted_ = false;
    double dist_factor_overlap_schwarz_ratio_cutoff_ = 1e-2;
    int count_ints_n_integral_extrema_ = 1;
    int n_iter_only_ = -1;
 
    std::shared_ptr<ScreeningStatistics> stats_;
    ScreeningStatistics::Iteration* iter_stats_;
 
    double prev_density_frob_;
    double prev_epsilon_;
    double prev_epsilon_dist_;
    double prev_epsilon_B_;
    double prev_epsilon_d_over_;
    double prev_epsilon_d_under_;
 
    RefDiagSCMatrix prev_evals_;
    RefSCMatrix prev_evecs_;
    std::vector<int> prev_occ_;
 
    int max_fxn_obs_ = 0;
    int max_fxn_obs_j_ish_ = 0;
    int max_fxn_obs_j_jsh_ = 0;
    int max_fxn_obs_todo_ = 0;
    int max_fxn_atom_obs_ = 0;
    int max_fxn_dfbs_ = 0;
    int max_fxn_dfbs_todo_ = 0;
    int max_fxn_atom_dfbs_ = 0;
    int max_obs_atom_fxn_on_dfbs_center_todo_ = 0;
    int max_fxn_atom_dfbs_todo_ = 0;
    //int max_fxn_obs_assigned_ = 0;
    //int max_fxn_atom_dfbs_assigned_ = 0;
 
    // for computing J using standard DF
    // only used if cadf_K_only == true
    Ref<WavefunctionWorld> world_;
 
    TwoCenterIntContainerPtr g2_full_ptr_;
 
    RefDiagSCMatrix core_evals_;
    RefSCMatrix core_evecs_;
 
    std::vector<int> orb_to_hcore_orb_;
 
    RefSCMatrix D_;
 
    RowMatrix S_frob_;
    RowMatrix dipole_frob_;
 
    // A (very) rough estimate of the minimum amount of memory that is in use by CADF at the present time
    std::atomic<size_t> memory_used_;
 
    class TrackedMemoryHolder {
        std::atomic<size_t>& used_ref_;
        size_t size_;
      public:
        TrackedMemoryHolder(
            std::atomic<size_t>& accum,
            size_t size
        ) : used_ref_(accum), size_(size)
        {
          used_ref_ += size_;
        }
        ~TrackedMemoryHolder() { used_ref_ -= size_; }
    };
    void consume_memory(size_t size) {
      memory_used_ += size;
    }
    void release_memory(size_t size) {
      memory_used_ -= size;
    }
    TrackedMemoryHolder hold_memory(size_t size) {
      return TrackedMemoryHolder(memory_used_, size);
    }
 
    // Convenience variable for better code readibility
    static const TwoBodyOper::type coulomb_oper_type_ = TwoBodyOper::eri;
    // Currently only metric_oper_type_ == coulomb_oper_type_ is supported
    TwoBodyOper::type metric_oper_type_;
 
    bool get_shell_pair(ShellData& mu, ShellData& nu, PairSet pset = AllPairs);
 
    RefSCMatrix compute_J();
    RefSCMatrix compute_J_cadf();
    RefSCMatrix compute_J_df();
 
    RefSCMatrix compute_K();
 
    void count_ints();
 
    RefSCMatrix new_compute_K();
 
    double get_distance_factor(
        const ShellData& ish,
        const ShellData& jsh,
        const ShellData& Xsh
    ) const;
 
    double get_R(
        const ShellData& ish,
        const ShellData& jsh,
        const ShellData& Xsh,
        bool ignore_extents = false
    ) const;
 
    void compute_coefficients();
 
    template<template<typename...> class container, typename map_type>
    void get_coefs_ish_jsh(
        const ShellData& ish,
        const ShellData& jsh,
        int ithr,
        container<map_type>& coefsA,
        container<map_type>& coefsB
    );
 
    void initialize();
 
    void init_significant_pairs();
 
    void init_threads();
 
    void done_threads();
 
    void print(std::ostream& o) const;
 
    TwoCenterIntContainerPtr ints_to_eigen(
        int ish, int jsh,
        Ref<TwoBodyTwoCenterInt>& ints,
        TwoBodyOper::type ints_type,
        const GaussianBasisSet* bs1=0,
        const GaussianBasisSet* bs2=0
    );
 
    template <typename ShellRange>
    TwoCenterIntContainerPtr ints_to_eigen(
        const ShellBlockData<ShellRange>& ish,
        const ShellBlockData<ShellRange>& jsh,
        Ref<TwoBodyTwoCenterInt>& ints,
        TwoBodyOper::type ints_type
    );
 
    template <typename ShellRange>
    TwoCenterIntContainerPtr ints_to_eigen_threaded(
        const ShellBlockData<ShellRange>& ish,
        const ShellBlockData<ShellRange>& jsh,
        std::vector<Ref<TwoBodyTwoCenterInt>>& ints_for_thread,
        TwoBodyOper::type ints_type
    );
 
    template <typename ShellRange>
    Eigen::Map<TwoCenterIntContainer>
    ints_to_eigen_map_threaded(
        const ShellBlockData<ShellRange>& iblk,
        const ShellBlockData<ShellRange>& jblk,
        std::vector<Ref<TwoBodyTwoCenterInt>>& ints_for_thread,
        TwoBodyOper::type int_type,
        double* buffer
    );
 
    template <typename ShellRange>
    ThreeCenterIntContainerPtr ints_to_eigen(
        const ShellBlockData<ShellRange>& ish, const ShellData& jsh,
        Ref<TwoBodyTwoCenterInt>& ints,
        TwoBodyOper::type ints_type
    );
 
    ThreeCenterIntContainerPtr ints_to_eigen(
        int ish, int jsh, int ksh,
        Ref<TwoBodyThreeCenterInt>& ints,
        TwoBodyOper::type ints_type,
        GaussianBasisSet* bs1 = 0,
        GaussianBasisSet* bs2 = 0,
        GaussianBasisSet* bs3 = 0
    );
 
    template <typename ShellRange>
    ThreeCenterIntContainerPtr ints_to_eigen(
        const ShellData& ish, const ShellData& jsh,
        const ShellBlockData<ShellRange>& kblk,
        Ref<TwoBodyThreeCenterInt>& ints,
        TwoBodyOper::type ints_type
    );
 
    template <typename ShellRange1, typename ShellRange2>
    ThreeCenterIntContainerPtr ints_to_eigen(
        const ShellBlockData<ShellRange1>& ish,
        const ShellData& jsh,
        const ShellBlockData<ShellRange2>& kblk,
        Ref<TwoBodyThreeCenterInt>& ints,
        TwoBodyOper::type ints_type
    );
 
    template <typename ShellRange>
    ThreeCenterIntContainerPtr ints_to_eigen(
        const ShellBlockData<ShellRange>& ish,
        const ShellData& jsh,
        const ShellData& ksh,
        Ref<TwoBodyThreeCenterInt>& ints,
        TwoBodyOper::type ints_type
    );
 
    template <typename ShellRange>
    Eigen::Map<ThreeCenterIntContainer> ints_to_eigen_map(
        const ShellData& ish,
        const ShellData& jsh,
        const ShellBlockData<ShellRange>& Xsh,
        Ref<TwoBodyThreeCenterInt>& ints,
        TwoBodyOper::type ints_type,
        double* __restrict__ buffer
    );
 
    template <typename MapType>
    void ints_to_eigen_map(
        const ShellData& ish,
        const ShellData& jsh,
        const ShellData& Xsh,
        Ref<TwoBodyThreeCenterInt>& ints,
        TwoBodyOper::type int_type,
        MapType& out_map
    );
 
    template <typename ShellRange>
    Eigen::Map<ThreeCenterIntContainer> ints_to_eigen_map(
        const ShellBlockData<ShellRange>& ish,
        const ShellData& jsh,
        const ShellData& ksh,
        Ref<TwoBodyThreeCenterInt>& ints,
        TwoBodyOper::type ints_type,
        double* buffer
    );
 
    template <typename ShellRange1, typename ShellRange2>
    Eigen::Map<ThreeCenterIntContainer> ints_to_eigen_map(
        const ShellBlockData<ShellRange1>& ish,
        const ShellData& jsh,
        const ShellBlockData<ShellRange2>& Xblk,
        Ref<TwoBodyThreeCenterInt>& ints,
        TwoBodyOper::type ints_type,
        double* buffer
    );
 
    void ints_to_buffer(
        int ish, int jsh, int ksh,
        int nbfi, int nbfj, int nbfk,
        Ref<TwoBodyThreeCenterInt>& ints,
        TwoBodyOper::type ints_type,
        double* buffer, int stride=-1
    );
 
    template <typename ShellRange>
    ThreeCenterIntContainerPtr ints_to_eigen(
        const ShellData& ish,
        const ShellBlockData<ShellRange>& jblk,
        const ShellData& ksh,
        Ref<TwoBodyThreeCenterInt>& ints,
        TwoBodyOper::type ints_type
    );
 
    FourCenterIntContainerPtr ints_to_eigen(
        const ShellData& ish, const ShellData& jsh,
        const ShellData& ksh, const ShellData& lsh,
        Ref<TwoBodyInt>& ints, TwoBodyOper::type ints_type
    );
 
    shared_ptr<Decomposition> get_decomposition(int ish, int jsh, Ref<TwoBodyTwoCenterInt> ints);
 
 
    // Non-blocked version of cl_gmat_
    RefSymmSCMatrix gmat_;
 
    bool density_reset_ = true;
 
    bool have_coefficients_;
 
    // The density fitting auxiliary basis
    Ref<GaussianBasisSet> dfbs_ = 0;
 
    // Where are the shell pairs being evaluated?
    std::map<PairSet, std::map<std::pair<int, int>, int>> pair_assignments_;
 
    // Pair assignments for K
    std::map<std::pair<int, ShellBlockSkeleton<>>, int> pair_assignments_k_;
 
    shared_ptr<cadf::AssignmentGrid> atom_pair_assignments_k_;
    shared_ptr<cadf::assignments::Assignments> assignments_new_k_;
 
    // What pairs are being evaluated on the current node?
    std::vector<std::pair<int, int>> local_pairs_all_;
    std::vector<std::pair<int, int>> local_pairs_sig_;
 
    // What pairs are being evaluated on the current node?
    std::vector<std::pair<int, ShellBlockSkeleton<>>> local_pairs_k_;
    std::unordered_set<
      std::pair<int, int>, sc::hash<std::pair<int, int>>
    > local_pairs_linK_;
    std::unordered_map<int, std::vector<int>> linK_local_map_;
    std::unordered_map<int, std::vector<int>> linK_local_map_ish_Xsh_;
 
    // List of the permutationally unique pairs with half-schwarz bounds larger than pair_thresh_
    std::vector<std::pair<int, int>> sig_pairs_;
 
    // Significant partners for a given shell index
    std::vector<std::set<int>> sig_partners_;
 
    // Significant partners in pre-blocked form
    std::vector<ContiguousShellBlockList> sig_partner_blocks_;
 
    Eigen::VectorXd schwarz_df_;
    Eigen::VectorXd schwarz_df_mine_;
 
    // Where are we in the iteration over the local_pairs_?
    std::atomic<int> local_pairs_spot_;
 
    Eigen::MatrixXd schwarz_frob_;
 
    Eigen::MatrixXd C_bar_;
    Eigen::MatrixXd C_bar_mine_;
    Eigen::MatrixXd C_underbar_;
    shared_ptr<cadf::TreeMatrix<>> C_dfsame_;
 
    // TwoBodyThreeCenterInt integral objects for each thread
    std::vector<Ref<TwoBodyThreeCenterInt>> eris_3c_;
    std::vector<Ref<TwoBodyThreeCenterInt>> metric_ints_3c_;
 
    // TwoBodyTwoCenterInt integral objects for each thread
    std::vector<Ref<TwoBodyTwoCenterInt>> eris_2c_;
    std::vector<Ref<TwoBodyTwoCenterInt>> metric_ints_2c_;
 
    // Count of integrals computed locally this iteration
    std::atomic<long> ints_computed_locally_;
 
    std::vector<Eigen::Vector3d> centers_;
 
    fast_map<std::pair<int, int>, Eigen::Vector3d> pair_centers_;
 
    // Extents of the various pairs
    fast_map<std::pair<int, int>, double> pair_extents_;
    fast_map<std::pair<int, int>, double> effective_pair_exponents_;
    std::vector<double> df_extents_;
    std::vector<double> effective_df_exponents_;
 
    double* __restrict__ coefficients_data_ = 0;
    double* __restrict__ dist_coefs_data_ = 0;
    double* __restrict__ dist_coefs_data_df_ = 0;
 
    fast_map<int, Eigen::Map<RowMatrix>> coefs_X_nu;
    fast_map<int, Eigen::Map<RowMatrix>> coefs_X_nu_other;
    fast_map<int, Eigen::Map<RowMatrix>> coefs_mu_X;
    fast_map<int, Eigen::Map<RowMatrix>> coefs_mu_X_other;
 
    CoefMap coefs_;
 
    // for each atom, matrix of mu_a in atom, rho_b, X(ab)
    std::vector<RowMatrix> coefs_blocked_;
    std::vector<std::vector<int>> coef_block_offsets_;
 
    std::vector<Eigen::Map<RowMatrix>> coefs_transpose_blocked_;
    std::vector<Eigen::Map<RowMatrix>> coefs_transpose_blocked_other_;
 
    std::vector<Eigen::Map<RowMatrix>> coefs_transpose_;
 
    shared_ptr<DecompositionCache> decomps_;
 
    bool threads_initialized_ = false;
 
    static ClassDesc cd_;
 
    typedef typename decltype(sig_partners_)::value_type::iterator sig_partners_iter_t;
    typedef range_of<ShellData, sig_partners_iter_t> sig_partners_range_t;
 
    sig_partners_range_t
    iter_significant_partners(const ShellData& ish){
      const auto& sig_parts = sig_partners_[ish];
      return boost::make_iterator_range(
          basis_element_iterator<ShellData, sig_partners_iter_t>(
              ish.basis, ish.dfbasis, sig_parts.begin()
          ),
          basis_element_iterator<ShellData, sig_partners_iter_t>(
              ish.basis, ish.dfbasis, sig_parts.end()
          )
      );
    }
 
    range_of_shell_blocks<sig_partners_iter_t>
    iter_significant_partners_blocked(
        const ShellData& ish,
        int requirements = Contiguous,
        int target_size = DEFAULT_TARGET_BLOCK_SIZE
    ){
      const auto& sig_parts = sig_partners_[ish];
      return boost::make_iterator_range(
          shell_block_iterator<sig_partners_iter_t>(
              sig_parts.begin(), sig_parts.end(),
              ish.basis, ish.dfbasis, requirements, target_size
          ),
          shell_block_iterator<sig_partners_iter_t>(
              sig_parts.end(), sig_parts.end(),
              ish.basis, ish.dfbasis, requirements, target_size
          )
      );
    }
 
    bool is_sig_pair(int ish, int jsh) const {
      const auto& parts = sig_partners_[ish];
      return parts.find(jsh) != parts.end();
    }
 
    // CADF-LinK lists
    typedef madness::ConcurrentHashMap<int, OrderedShellList> IndexListMap;
    typedef madness::ConcurrentHashMap<
        std::pair<int, int>,
        OrderedShellList,
        sc::hash<std::pair<int, int>>
    > IndexListMap2;
    typedef madness::ConcurrentHashMap<
        std::pair<int, int>,
        std::vector<std::pair<uint64_t, uint64_t>>,
        sc::hash<std::pair<int, int>>
    > LinKRangeMap2;
 
    IndexListMap L_schwarz;
    IndexListMap L_DC;
    IndexListMap L_C_under;
    IndexListMap2 L_3;
    IndexListMap2 L_3_star;
    IndexListMap2 L_B;
    IndexListMap2 L_d_over;
    LinKRangeMap2 L_d_under_ranges;
 
    friend class ApproximatePairWriter;
 
    Ref<ApproximatePairWriter> pair_writer_;
 
}; // CADFCLHF
 
// end of addtogroup ChemistryElectronicStructureOneBodyHFCADF
 
void
get_split_range_part(
    const Ref<MessageGrp>& msg,
    int full_begin, int full_end,
    int& out_begin, int& out_size
);
 
void
get_split_range_part(
    int me, int n,
    int full_begin, int full_end,
    int& out_begin, int& out_size
);
 
boost::property_tree::ptree&
write_xml(
    const CADFCLHF::ScreeningStatistics& obj,
    boost::property_tree::ptree& parent,
    const XMLWriter& writer
);
 
boost::property_tree::ptree&
write_xml(
    const CADFCLHF::ScreeningStatistics::Iteration& obj,
    boost::property_tree::ptree& parent,
    const XMLWriter& writer
);
 
// Borrowed from ConsumeableResources
inline std::string data_size_to_string(size_t t) {
  const int prec = 3; // print this many digits
 
  // determine m such that 1024^m <= t <= 1024^(m+1)
  char m = 0;
  double thousand_m = 1;
  double thousand_mp1 = 1024;
  while (t >= thousand_mp1) {
    ++m;
    thousand_m = thousand_mp1;
    thousand_mp1 *= 1024;
  }
 
  // determine units
  std::string unit;
  switch (m) {
    case 0: unit = "B"; break;
    case 1: unit = "kB"; break;
    case 2: unit = "MB"; break;
    case 3: unit = "GB"; break;
    case 4: unit = "TB"; break;
    case 5: unit = "PB"; break;
    case 6: unit = "EB"; break;
    case 7: unit = "ZB"; break;
    case 8: unit = "YB"; break;
    default: MPQC_ASSERT(false); break;
  }
 
  // compute normalized mantissa
  std::ostringstream oss;
  oss.precision(prec);
  oss << (double)t/(double)thousand_m << unit;
  return oss.str();
}
 
} // end namespace sc
 
#include "get_ints_impl.h"
#include "coefs_impl.h"
 
#endif /* _chemistry_qc_scf_cadfclhf_h */