Included new HamiltonianGenerator, particularly efficient for impurit…

…y models. Test case added.

include/macis/bitset_operations.hpp

@@ -7,6 +7,9 @@
  */
 
 #pragma once
+#include <string.h>
+#include <strings.h>
+
 #include <bit>
 #include <cassert>
 #include <climits>

include/macis/hamiltonian_generator.hpp

@@ -183,7 +183,7 @@ class HamiltonianGenerator {
   void rotate_hamiltonian_ordm(const double* ordm);
 
   virtual void SetJustSingles(bool /*_js*/) {}
-  virtual bool GetJustSingles() { return false; }
+  virtual bool GetJustSingles() const { return false; }
   virtual size_t GetNimp() const { return N / 2; }
 };
 

include/macis/hamiltonian_generator/sd_build.hpp

@@ -0,0 +1,274 @@
+/*
+ * MACIS Copyright (c) 2023, The Regents of the University of California,
+ * through Lawrence Berkeley National Laboratory (subject to receipt of
+ * any required approvals from the U.S. Dept. of Energy). All rights reserved.
+ *
+ * See LICENSE.txt for details
+ */
+
+#pragma once
+#include <macis/hamiltonian_generator.hpp>
+#include <macis/sd_operations.hpp>
+#include <macis/util/rdms.hpp>
+#include <set>
+
+namespace macis {
+
+template< size_t N >
+struct det_pos
+{
+  public:
+    std::bitset<N> det;
+    uint32_t  id;
+};
+
+template< size_t N >
+bool operator<( const det_pos<N>& a, const det_pos<N>& b )
+{ return bitset_less<N>( a.det, b.det ); }
+
+template <size_t N>
+class SDBuildHamiltonianGenerator : public HamiltonianGenerator<N> {
+ public:
+  using base_type = HamiltonianGenerator<N>;
+  using full_det_t = typename base_type::full_det_t;
+  using spin_det_t = typename base_type::spin_det_t;
+  using full_det_iterator = typename base_type::full_det_iterator;
+  using matrix_span_t = typename base_type::matrix_span_t;
+  using rank4_span_t = typename base_type::rank4_span_t;
+
+  template <typename index_t>
+  using sparse_matrix_type = sparsexx::csr_matrix<double, index_t>;
+
+ protected:
+  size_t nimp, nimp2, nimp3;
+
+  template <typename index_t>
+  sparse_matrix_type<index_t> make_csr_hamiltonian_block_(
+      full_det_iterator bra_begin, full_det_iterator bra_end,
+      full_det_iterator ket_begin, full_det_iterator ket_end, double H_thresh) {
+    const size_t nbra_dets = std::distance(bra_begin, bra_end);
+    const size_t nket_dets = std::distance(ket_begin, ket_end);
+
+    std::vector<index_t> colind, rowptr(nbra_dets + 1);
+    std::vector<double> nzval;
+
+    // List of impurity orbitals, assumed to be the first nimp.
+    std::vector<uint32_t> imp_orbs(nimp, 0);
+    for( int ii = 0; ii < nimp; ii++)
+      imp_orbs[ii] = ii;
+    std::vector<uint32_t> bra_occ_alpha, bra_occ_beta;
+
+    std::set<det_pos<N> > kets;
+    for( full_det_iterator it = ket_begin; it != ket_end; it++ )
+    { 
+      det_pos<N> a;
+      a.det = *it;
+      a.id  = std::distance(ket_begin, it);
+      kets.insert( a );
+    }
+
+    rowptr[0] = 0;
+
+    // Loop over bra determinants
+    for(size_t i = 0; i < nbra_dets; ++i) {
+      //if( (i%1000) == 0 ) std::cout << i << ", " << rowptr[i] << std::endl;
+      const auto bra = *(bra_begin + i);
+
+      size_t nrow = 0;
+      if(bra.count()) {
+        // Separate out into alpha/beta components
+        spin_det_t bra_alpha = bitset_lo_word(bra);
+        spin_det_t bra_beta  = bitset_hi_word(bra);
+
+        // Get occupied indices
+        bits_to_indices(bra_alpha, bra_occ_alpha);
+        bits_to_indices(bra_beta , bra_occ_beta);
+
+        // Get singles and doubles
+        // (Note that doubles only involve impurity orbitals)
+        std::vector<full_det_t> excs, doubles;
+        if( just_singles )
+          generate_singles_spin( this->norb_, bra, excs );
+        else
+        {
+          std::vector<full_det_t> singls;
+          generate_singles_spin( this->norb_, bra, excs );
+          // This will store in singls sinles among impurity orbitals, which we
+          // have already taken into account. 
+          generate_singles_doubles_spin_as( this->norb_, bra, singls, doubles, imp_orbs );
+          excs.insert( excs.end(), doubles.begin(), doubles.end() );
+        }
+
+        // Diagonal term
+        full_det_t ex_diag   = bra ^ bra;
+        spin_det_t exd_alpha = bitset_lo_word( ex_diag );
+        spin_det_t exd_beta  = bitset_hi_word( ex_diag );
+
+        // Compute Matrix Element
+        const auto h_eld = this->matrix_element_diag( bra_occ_alpha, bra_occ_beta );
+
+        if( std::abs(h_eld) > H_thresh ) {
+          nrow++;
+          colind.emplace_back(i);
+          nzval.emplace_back(h_eld);
+        }
+
+        // Loop over ket determinants
+        for( const auto pos_ket : excs ) {
+          det_pos<N> pos_ket2; pos_ket2.det = pos_ket; pos_ket2.id = 0;
+          auto it = kets.find( pos_ket2 );
+          if( it != kets.end() ) {
+            int j = it->id;
+            spin_det_t ket_alpha = bitset_lo_word( pos_ket );
+            spin_det_t ket_beta  = bitset_hi_word( pos_ket );
+
+            full_det_t ex_total = bra ^ pos_ket;
+
+            spin_det_t ex_alpha = bitset_lo_word( ex_total );
+            spin_det_t ex_beta  = bitset_hi_word( ex_total );
+
+            // Compute Matrix Element
+            const auto h_el = this->matrix_element( 
+	        bra_alpha, ket_alpha, ex_alpha, bra_beta, ket_beta, ex_beta, 
+		bra_occ_alpha, bra_occ_beta );
+
+            if( std::abs(h_el) > H_thresh ) {
+              nrow++;
+              colind.emplace_back(j);
+              nzval.emplace_back(h_el);
+            }
+
+
+          } // Non-zero ket determinant
+        } // Loop over ket determinants
+
+      }  // Non-zero bra determinant
+
+      rowptr[i + 1] = rowptr[i] + nrow;  // Update rowptr
+
+    }  // Loop over bra determinants
+
+    colind.shrink_to_fit();
+    nzval.shrink_to_fit();
+
+    return sparse_matrix_type<index_t>(nbra_dets, nket_dets, std::move(rowptr),
+                                       std::move(colind), std::move(nzval));
+  }
+
+  sparse_matrix_type<int32_t> make_csr_hamiltonian_block_32bit_(
+      full_det_iterator bra_begin, full_det_iterator bra_end,
+      full_det_iterator ket_begin, full_det_iterator ket_end,
+      double H_thresh) override {
+    return make_csr_hamiltonian_block_<int32_t>(bra_begin, bra_end, ket_begin,
+                                                ket_end, H_thresh);
+  }
+
+  sparse_matrix_type<int64_t> make_csr_hamiltonian_block_64bit_(
+      full_det_iterator bra_begin, full_det_iterator bra_end,
+      full_det_iterator ket_begin, full_det_iterator ket_end,
+      double H_thresh) override {
+    return make_csr_hamiltonian_block_<int64_t>(bra_begin, bra_end, ket_begin,
+                                                ket_end, H_thresh);
+  }
+
+ public:
+  void form_rdms(full_det_iterator bra_begin, full_det_iterator bra_end,
+                 full_det_iterator ket_begin, full_det_iterator ket_end,
+                 double *C, matrix_span_t ordm, rank4_span_t trdm) override {
+    const size_t nbra_dets = std::distance(bra_begin, bra_end);
+    const size_t nket_dets = std::distance(ket_begin, ket_end);
+
+    std::vector<uint32_t> bra_occ_alpha, bra_occ_beta;
+
+    std::set<det_pos<N> > kets;
+    for( full_det_iterator it = ket_begin; it != ket_end; it++ )
+    { 
+      det_pos<N> a;
+      a.det = *it;
+      a.id  = std::distance(ket_begin, it);
+      kets.insert( a );
+    }
+
+    // Loop over bra determinants
+    for(size_t i = 0; i < nbra_dets; ++i) {
+      const auto bra = *(bra_begin + i);
+      // if( (i%1000) == 0 ) std::cout << i  << std::endl;
+      if(bra.count()) {
+        // Separate out into alpha/beta components
+        spin_det_t bra_alpha = bitset_lo_word(bra);
+        spin_det_t bra_beta  = bitset_hi_word(bra);
+
+        // Get occupied indices
+        bits_to_indices(bra_alpha, bra_occ_alpha);
+        bits_to_indices(bra_beta, bra_occ_beta);
+
+        // Get singles and doubles
+        std::vector<full_det_t> excs;
+	if( trdm.data_handle() ){
+          std::vector<full_det_t> doubles;
+          generate_singles_doubles_spin( this->norb_, bra, excs, doubles );
+          excs.insert( excs.end(), doubles.begin(), doubles.end() );
+	}
+	else{
+          generate_singles_spin( this->norb_, bra, excs );
+	}
+
+        // Diagonal term
+        full_det_t ex_diag   = bra ^ bra;
+        spin_det_t exd_alpha = bitset_lo_word( ex_diag );
+        spin_det_t exd_beta  = bitset_hi_word( ex_diag );
+
+        // Compute Matrix Element
+        rdm_contributions( bra_alpha, bra_alpha,
+          exd_alpha, bra_beta, bra_beta, exd_beta, bra_occ_alpha,
+          bra_occ_beta, C[i] * C[i], ordm, trdm );
+
+        // Loop over excitations
+        for( const auto pos_ket : excs ) {
+          det_pos<N> pos_ket2; pos_ket2.det = pos_ket; pos_ket2.id = 0;
+          auto it = kets.find( pos_ket2 );
+          if( it != kets.end() ) {
+            int j = it->id;
+            spin_det_t ket_alpha = bitset_lo_word( pos_ket );
+            spin_det_t ket_beta  = bitset_hi_word( pos_ket );
+
+            full_det_t ex_total = bra ^ pos_ket;
+	    int ex_lim = 2;
+	    if(trdm.data_handle()) ex_lim = 4;
+            if( ex_total.count() <= ex_lim ) {
+
+              spin_det_t ex_alpha = bitset_lo_word( ex_total );
+              spin_det_t ex_beta  = bitset_hi_word( ex_total );
+
+              const double val = C[i] * C[j];
+
+              // Compute Matrix Element
+              rdm_contributions( bra_alpha, ket_alpha,
+                ex_alpha, bra_beta, ket_beta, ex_beta, bra_occ_alpha,
+                bra_occ_beta, val, ordm, trdm );
+
+            } // Possible non-zero connection (Hamming distance)
+
+          } // Non-zero ket determinant
+        } // Loop over ket determinants
+
+      } // Non-zero bra determinant
+    } // Loop over bra determinants 
+  }
+
+ public:
+
+  bool just_singles;
+
+  template <typename... Args>
+  SDBuildHamiltonianGenerator(Args &&...args)
+      : HamiltonianGenerator<N>(std::forward<Args>(args)...), just_singles(false)
+        { SetNimp(this->norb_); }
+
+  void SetJustSingles( bool _js ) override{ just_singles = _js; }
+  void SetNimp( size_t _n ){ nimp = _n; nimp2 = _n * _n; nimp3 = nimp2 * _n; }
+  size_t GetNimp() const override { return nimp; }
+  bool GetJustSingles( ) const override{ return just_singles; }
+};
+
+}  // namespace macis

include/macis/util/fcidump.hpp

@@ -70,6 +70,17 @@ void read_fcidump_1body(std::string fname, col_major_span<double, 2> T);
  */
 void read_fcidump_2body(std::string fname, col_major_span<double, 4> V);
 
+/**
+ * @brief Check whether the 2-body contribution of the Hamiltonian is exclusively
+ *        diagonal.
+ *
+ * @param[in] fname: Filename of FCIDUMP file
+ *
+ * @returns bool: Is the 2-body contribution to the Hamiltonian exclusively
+ *                diagonal?
+ */
+bool is_2body_diagonal(std::string fname);
+
 /**
  *  @brief Write an FCIDUMP file from a 2-body hamiltonian
  *

src/macis/fcidump.cxx

@@ -181,6 +181,34 @@ void read_fcidump_2body(std::string fname, double* V, size_t LDV) {
       fname, KokkosEx::submdspan(V_map, sl, sl, sl, Kokkos::full_extent));
 }
 
+bool is_2body_diagonal(std::string fname) {
+
+  auto norb     = read_fcidump_norb(fname);
+  bool all_diag = true;
+
+  std::ifstream file(fname);
+  std::string line;
+  while(std::getline(file, line)) {
+    auto tokens = split(line, " ");
+    if(tokens.size() != 5) continue;  // not a valid FCIDUMP line
+
+    auto [p, q, r, s, integral] = fcidump_line(tokens);
+    auto lc = line_classification(p, q, r, s);
+    if(lc == LineClassification::TwoBody) {
+      p--;
+      q--;
+      r--;
+      s--;
+      if(p != q || r != s)
+      {
+	all_diag = false;
+	break;
+      }
+    }
+  }
+  return all_diag;
+}
+
 void write_fcidump(std::string fname, size_t norb, const double* T, size_t LDT,
                    const double* V, size_t LDV, double E_core) {
   auto logger = spdlog::basic_logger_mt("fcidump", fname);

tests/CMakeLists.txt

@@ -13,6 +13,7 @@ add_executable( macis_test
   fcidump.cxx 
   read_wavefunction.cxx
   double_loop.cxx
+  sd_build.cxx
   csr_hamiltonian.cxx
   davidson.cxx
   transform.cxx

tests/ref_data/hubbard10.fci.dat

@@ -0,0 +1,38 @@
+1  1  1  1      4.
+2  2  2  2      4.
+3  3  3  3      4.
+4  4  4  4      4.
+5  5  5  5      4.
+6  6  6  6      4.
+7  7  7  7      4.
+8  8  8  8      4.
+9  9  9  9      4.
+10 10 10 10     4.
+1  1  0  0     -2.
+2  2  0  0     -2.
+3  3  0  0     -2.
+4  4  0  0     -2.
+5  5  0  0     -2.
+6  6  0  0     -2.
+7  7  0  0     -2.
+8  8  0  0     -2.
+9  9  0  0     -2.
+10 10 0  0     -2.
+1  2  0  0     -1.
+2  1  0  0     -1.
+2  3  0  0     -1.
+3  2  0  0     -1.
+3  4  0  0     -1.
+4  3  0  0     -1.
+4  5  0  0     -1.
+5  4  0  0     -1.
+5  6  0  0     -1.
+6  5  0  0     -1.
+6  7  0  0     -1.
+7  6  0  0     -1.
+7  8  0  0     -1.
+8  7  0  0     -1.
+8  9  0  0     -1.
+9  8  0  0     -1.
+9  10 0  0     -1.
+10 9  0  0     -1.