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..37c8cf0f --- /dev/null +++ b/include/macis/hamiltonian_generator/sd_build.hpp @@ -0,0 +1,278 @@ +/* + * 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 +struct det_pos { + public: + std::bitset det; + uint32_t id; +}; + +template +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/model/hubbard.hpp b/include/macis/model/hubbard.hpp new file mode 100644 index 00000000..5e80fd0f --- /dev/null +++ b/include/macis/model/hubbard.hpp @@ -0,0 +1,21 @@ +/* + * 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 + +namespace macis { + +/** + * @brief Generate Hamiltonian Data for 1D Hubbard + */ +void hubbard_1d(size_t nsites, double t, double U, std::vector& T, + std::vector& V, bool pbc = false); + +} // namespace macis diff --git a/include/macis/util/fcidump.hpp b/include/macis/util/fcidump.hpp index 41bb2a5c..823ef7c8 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/CMakeLists.txt b/src/macis/CMakeLists.txt index 26155ac6..74c1995a 100644 --- a/src/macis/CMakeLists.txt +++ b/src/macis/CMakeLists.txt @@ -15,6 +15,7 @@ add_library( macis orbital_rotation_utilities.cxx orbital_hessian.cxx orbital_steps.cxx + model/hubbard.cxx ) target_include_directories( macis PUBLIC diff --git a/src/macis/fcidump.cxx b/src/macis/fcidump.cxx index bf1d7762..9870077c 100644 --- a/src/macis/fcidump.cxx +++ b/src/macis/fcidump.cxx @@ -181,6 +181,32 @@ 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/src/macis/model/hubbard.cxx b/src/macis/model/hubbard.cxx new file mode 100644 index 00000000..5e74b0d3 --- /dev/null +++ b/src/macis/model/hubbard.cxx @@ -0,0 +1,39 @@ +/* + * 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 + +namespace macis { + +void hubbard_1d(size_t nsites, double t, double U, std::vector& T, + std::vector& V, bool pbc) { + T.resize(nsites * nsites); + V.resize(nsites * nsites * nsites * nsites); + + for(size_t p = 0; p < nsites; ++p) { + // Half-filling Chemical Potential + T[p * (nsites + 1)] = -U / 2; + + // On-Site Interaction + V[p * (nsites * nsites * nsites + nsites * nsites + nsites + 1)] = U; + + // Hopping + if(p < nsites - 1) { + T[p + (p + 1) * nsites] = -t; + T[(p + 1) + p * nsites] = -t; + } + } + + // PBC for 1-D + if(pbc) { + T[(nsites - 1)] = -t; + T[(nsites - 1) * nsites] = -t; + } +} + +} // namespace macis diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 038f9dfd..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 @@ -20,6 +21,7 @@ add_executable( macis_test mcscf.cxx asci.cxx dist_quickselect.cxx + hubbard.cxx ) target_link_libraries( macis_test PUBLIC macis Catch2::Catch2 ) target_include_directories( macis_test PUBLIC ${PROJECT_BINARY_DIR}/tests ) diff --git a/tests/hubbard.cxx b/tests/hubbard.cxx new file mode 100644 index 00000000..6f3d3553 --- /dev/null +++ b/tests/hubbard.cxx @@ -0,0 +1,64 @@ +/* + * 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 "ut_common.hpp" + +TEST_CASE("Hubbard") { + ROOT_ONLY(MPI_COMM_WORLD); + + const size_t nsites = 4; + const size_t nsites2 = nsites * nsites; + const size_t nsites3 = nsites2 * nsites; + const double t = 1.0, U = 4.0; + std::vector T, V; + + SECTION("1D") { + bool pbc; + + SECTION("No PBC") { pbc = false; } + + SECTION("PBC") { pbc = true; } + + macis::hubbard_1d(nsites, t, U, T, V, pbc); + + // Check two-body term + for(int p = 0; p < nsites; ++p) + for(int q = 0; q < nsites; ++q) + for(int r = 0; r < nsites; ++r) + for(int s = 0; s < nsites; ++s) { + const auto mat_el = V[p + nsites * q + nsites2 * r + nsites3 * s]; + if(p == q and p == r and p == s) + REQUIRE(mat_el == U); + else + REQUIRE(mat_el == 0.0); + } + + // Check 1-body term + for(int p = 0; p < nsites; ++p) + for(int q = 0; q < nsites; ++q) { + const auto mat_el = T[p + q * nsites]; + if(p == q) + REQUIRE(mat_el == Approx(-U / 2)); + else if(std::abs(p - q) == 1) + REQUIRE(mat_el == -t); + else if(pbc) { + if((p == 0 and q == nsites - 1) or (p == nsites - 1 and q == 0)) + REQUIRE(mat_el == -t); + else + REQUIRE(mat_el == 0.0); + } else + REQUIRE(mat_el == 0.0); + } + } +} 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..9c2a0760 --- /dev/null +++ b/tests/sd_build.cxx @@ -0,0 +1,116 @@ +/* + * 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::NumActive; +using macis::NumInactive; +using macis::NumOrbital; + +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";