diff --git a/include/macis/bitset_operations.hpp b/include/macis/bitset_operations.hpp index e0f5ec27..e11b5f69 100644 --- a/include/macis/bitset_operations.hpp +++ b/include/macis/bitset_operations.hpp @@ -7,6 +7,9 @@ */ #pragma once +#include +#include + #include #include #include diff --git a/include/macis/hamiltonian_generator.hpp b/include/macis/hamiltonian_generator.hpp index 9ef9583b..899f2586 100644 --- a/include/macis/hamiltonian_generator.hpp +++ b/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; } }; diff --git a/include/macis/hamiltonian_generator/sd_build.hpp b/include/macis/hamiltonian_generator/sd_build.hpp new file mode 100644 index 00000000..b1170e79 --- /dev/null +++ b/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 +#include +#include +#include + +namespace macis { + +template< size_t N > +struct det_pos +{ + public: + std::bitset det; + uint32_t id; +}; + +template< size_t N > +bool operator<( const det_pos& a, const det_pos& b ) +{ return bitset_less( a.det, b.det ); } + +template +class SDBuildHamiltonianGenerator : public HamiltonianGenerator { + public: + using base_type = HamiltonianGenerator; + 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 + using sparse_matrix_type = sparsexx::csr_matrix; + + protected: + size_t nimp, nimp2, nimp3; + + template + sparse_matrix_type 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 colind, rowptr(nbra_dets + 1); + std::vector nzval; + + // List of impurity orbitals, assumed to be the first nimp. + std::vector imp_orbs(nimp, 0); + for( int ii = 0; ii < nimp; ii++) + imp_orbs[ii] = ii; + std::vector bra_occ_alpha, bra_occ_beta; + + std::set > kets; + for( full_det_iterator it = ket_begin; it != ket_end; it++ ) + { + det_pos 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 excs, doubles; + if( just_singles ) + generate_singles_spin( this->norb_, bra, excs ); + else + { + std::vector 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 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(nbra_dets, nket_dets, std::move(rowptr), + std::move(colind), std::move(nzval)); + } + + sparse_matrix_type 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_(bra_begin, bra_end, ket_begin, + ket_end, H_thresh); + } + + sparse_matrix_type 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_(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 bra_occ_alpha, bra_occ_beta; + + std::set > kets; + for( full_det_iterator it = ket_begin; it != ket_end; it++ ) + { + det_pos 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 excs; + if( trdm.data_handle() ){ + std::vector 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 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 + SDBuildHamiltonianGenerator(Args &&...args) + : HamiltonianGenerator(std::forward(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 diff --git a/include/macis/util/fcidump.hpp b/include/macis/util/fcidump.hpp index 41bb2a5c..6821f434 100644 --- a/include/macis/util/fcidump.hpp +++ b/include/macis/util/fcidump.hpp @@ -70,6 +70,17 @@ void read_fcidump_1body(std::string fname, col_major_span T); */ void read_fcidump_2body(std::string fname, col_major_span 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 * diff --git a/src/macis/fcidump.cxx b/src/macis/fcidump.cxx index bf1d7762..6af87371 100644 --- a/src/macis/fcidump.cxx +++ b/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); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 3fd01966..1e539475 100644 --- a/tests/CMakeLists.txt +++ b/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 diff --git a/tests/ref_data/hubbard10.fci.dat b/tests/ref_data/hubbard10.fci.dat new file mode 100644 index 00000000..aed59ff2 --- /dev/null +++ b/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. diff --git a/tests/sd_build.cxx b/tests/sd_build.cxx new file mode 100644 index 00000000..0a8c95cb --- /dev/null +++ b/tests/sd_build.cxx @@ -0,0 +1,115 @@ +/* + * 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 + */ + +#include +#include +#include +#include +#include +#include +#include + +#include "ut_common.hpp" + +using macis::NumOrbital; +using macis::NumActive; +using macis::NumInactive; + +TEST_CASE("Single Double Build") { + ROOT_ONLY(MPI_COMM_WORLD); + + auto norb = macis::read_fcidump_norb(hubbard10_fcidump); + const auto norb2 = norb * norb; + const auto norb3 = norb2 * norb; + const size_t nocc = 10; + + std::vector T(norb * norb); + std::vector V(norb * norb * norb * norb); + auto E_core = macis::read_fcidump_core(hubbard10_fcidump); + macis::read_fcidump_1body(hubbard10_fcidump, T.data(), norb); + macis::read_fcidump_2body(hubbard10_fcidump, V.data(), norb); + bool just_singles = macis::is_2body_diagonal( hubbard10_fcidump ); + + using generator_type = macis::SDBuildHamiltonianGenerator<64>; + +#if 0 + generator_type ham_gen(norb, V.data(), T.data()); +#else + generator_type ham_gen( + macis::matrix_span(T.data(), norb, norb), + macis::rank4_span(V.data(), norb, norb, norb, norb)); + ham_gen.SetJustSingles(just_singles); +#endif + const auto hf_det = macis::canonical_hf_determinant<64>(nocc, nocc); + + std::vector eps(norb); + for(auto p = 0ul; p < norb; ++p) { + double tmp = 0.; + for(auto i = 0ul; i < nocc; ++i) { + tmp += 2. * V[p * (norb + 1) + i * (norb2 + norb3)] - + V[p * (1 + norb3) + i * (norb + norb2)]; + } + eps[p] = T[p * (norb + 1)] + tmp; + } + const auto EHF = ham_gen.matrix_element(hf_det, hf_det); + + SECTION("HF Energy") { REQUIRE(EHF + E_core == Approx(-0.00)); } + + SECTION("RDM") { + std::vector ordm(norb * norb, 0.0), trdm(norb3 * norb, 0.0); + std::vector> dets = { + macis::canonical_hf_determinant<64>(nocc, nocc)}; + + std::vector C = {1.}; + + ham_gen.form_rdms( + dets.begin(), dets.end(), dets.begin(), dets.end(), C.data(), + macis::matrix_span(ordm.data(), norb, norb), + macis::rank4_span(trdm.data(), norb, norb, norb, norb)); + + auto E_tmp = blas::dot(norb2, ordm.data(), 1, T.data(), 1) + + blas::dot(norb3 * norb, trdm.data(), 1, V.data(), 1); + REQUIRE(E_tmp == Approx(EHF)); + } + + SECTION("CI") { + size_t nalpha = 5; + size_t nbeta = 5; + size_t n_inactive = 0; + size_t n_active = 10; + size_t n_virtual = 0; + std::vector C_local; + std::vector active_ordm(n_active * n_active); + std::vector active_trdm; + macis::MCSCFSettings mcscf_settings; + mcscf_settings.ci_max_subspace = 100; + // Copy integrals into active subsets + std::vector T_active(n_active * n_active); + std::vector V_active(n_active * n_active * n_active * n_active); + + // Compute active-space Hamiltonian and inactive Fock matrix + std::vector F_inactive(norb2); + macis::active_hamiltonian(NumOrbital(norb), NumActive(n_active), + NumInactive(n_inactive), T.data(), norb, V.data(), + norb, F_inactive.data(), norb, T_active.data(), + n_active, V_active.data(), n_active); + auto E_inactive = macis::inactive_energy(NumInactive(n_inactive), T.data(), + norb, F_inactive.data(), norb); + auto dets = macis::generate_hilbert_space<64>(norb, nalpha, nbeta); + double E0 = + macis::selected_ci_diag(dets.begin(), dets.end(), ham_gen, mcscf_settings.ci_matel_tol, + mcscf_settings.ci_max_subspace, mcscf_settings.ci_res_tol, C_local, + MACIS_MPI_CODE(MPI_COMM_WORLD,) true); + E0 += E_inactive + E_core; + // Compute RDMs + ham_gen.form_rdms(dets.begin(), dets.end(), dets.begin(), dets.end(), + C_local.data(), macis::matrix_span(active_ordm.data(), norb, norb), + macis::rank4_span(active_trdm.data(), norb, norb, norb, norb)); + REQUIRE(E0 == Approx(-2.538061882041e+01)); + } +} diff --git a/tests/ut_common.hpp.in b/tests/ut_common.hpp.in index cb9cc085..ba5af02c 100644 --- a/tests/ut_common.hpp.in +++ b/tests/ut_common.hpp.in @@ -37,3 +37,5 @@ const std::string water_ccpvdz_nzval_fname = REF_DATA_PREFIX "/h2o.ccpvdz.cisd.nzval.bin"; const std::string ch4_wfn_fname = REF_DATA_PREFIX "/ch4.wfn.dat"; +const std::string hubbard10_fcidump = + REF_DATA_PREFIX "/hubbard10.fci.dat";