From 5175a96907195fac413bc85c8dc5a5aeb38e743e Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Wed, 7 Jun 2023 08:56:32 -0700 Subject: [PATCH 1/7] Added Hubbard Hamiltonian generator + UT --- include/macis/model/hubbard.hpp | 13 +++++++++ src/macis/CMakeLists.txt | 1 + src/macis/model/hubbard.cxx | 29 ++++++++++++++++++++ tests/CMakeLists.txt | 1 + tests/hubbard.cxx | 48 +++++++++++++++++++++++++++++++++ 5 files changed, 92 insertions(+) create mode 100644 include/macis/model/hubbard.hpp create mode 100644 src/macis/model/hubbard.cxx create mode 100644 tests/hubbard.cxx diff --git a/include/macis/model/hubbard.hpp b/include/macis/model/hubbard.hpp new file mode 100644 index 00000000..32ee3d71 --- /dev/null +++ b/include/macis/model/hubbard.hpp @@ -0,0 +1,13 @@ +#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); + +} diff --git a/src/macis/CMakeLists.txt b/src/macis/CMakeLists.txt index 862df2ce..dd65d431 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/model/hubbard.cxx b/src/macis/model/hubbard.cxx new file mode 100644 index 00000000..528349b0 --- /dev/null +++ b/src/macis/model/hubbard.cxx @@ -0,0 +1,29 @@ +#include + +namespace macis { + +void hubbard_1d(size_t nsites, double t, double U, std::vector& T, + std::vector& V) { + + + 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; + } + } + +} + +} diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 9d613244..73694653 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -20,6 +20,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..2f142ab5 --- /dev/null +++ b/tests/hubbard.cxx @@ -0,0 +1,48 @@ + + +#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") { + macis::hubbard_1d( nsites, t, U, T, V ); + + // 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 + REQUIRE(mat_el == 0.0); + } + } + +} From ebb7ef4d4ee3f92b3b0d23260cfd5cef3ee87f12 Mon Sep 17 00:00:00 2001 From: Clang Robot Date: Wed, 7 Jun 2023 15:57:12 +0000 Subject: [PATCH 2/7] Committing clang-format changes --- include/macis/model/hubbard.hpp | 6 ++--- src/macis/model/hubbard.cxx | 18 +++++-------- tests/hubbard.cxx | 47 ++++++++++++++++----------------- 3 files changed, 33 insertions(+), 38 deletions(-) diff --git a/include/macis/model/hubbard.hpp b/include/macis/model/hubbard.hpp index 32ee3d71..4fd849db 100644 --- a/include/macis/model/hubbard.hpp +++ b/include/macis/model/hubbard.hpp @@ -4,10 +4,10 @@ namespace macis { -/** +/** * @brief Generate Hamiltonian Data for 1D Hubbard */ -void hubbard_1d(size_t nsites, double t, double U, std::vector& T, +void hubbard_1d(size_t nsites, double t, double U, std::vector& T, std::vector& V); -} +} // namespace macis diff --git a/src/macis/model/hubbard.cxx b/src/macis/model/hubbard.cxx index 528349b0..cdbfa2d6 100644 --- a/src/macis/model/hubbard.cxx +++ b/src/macis/model/hubbard.cxx @@ -2,28 +2,24 @@ namespace macis { -void hubbard_1d(size_t nsites, double t, double U, std::vector& T, +void hubbard_1d(size_t nsites, double t, double U, std::vector& T, std::vector& V) { - - 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; + T[p * (nsites + 1)] = -U / 2; // On-Site Interaction - V[p * (nsites*nsites*nsites + nsites*nsites + nsites + 1)] = U; + 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; + if(p < nsites - 1) { + T[p + (p + 1) * nsites] = -t; + T[(p + 1) + p * nsites] = -t; } } - } -} +} // namespace macis diff --git a/tests/hubbard.cxx b/tests/hubbard.cxx index 2f142ab5..1f62f1a9 100644 --- a/tests/hubbard.cxx +++ b/tests/hubbard.cxx @@ -1,9 +1,9 @@ -#include #include #include +#include #include #include "ut_common.hpp" @@ -12,37 +12,36 @@ 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 size_t nsites2 = nsites * nsites; + const size_t nsites3 = nsites2 * nsites; const double t = 1.0, U = 4.0; - std::vector T,V; + std::vector T, V; SECTION("1D") { - macis::hubbard_1d( nsites, t, U, T, V ); + macis::hubbard_1d(nsites, t, U, T, V); // 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); - } + 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 - REQUIRE(mat_el == 0.0); - } + 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 + REQUIRE(mat_el == 0.0); + } } - } From 84b219da94e3bd76dafac1199b408682e86b8e15 Mon Sep 17 00:00:00 2001 From: "license[bot]" Date: Wed, 7 Jun 2023 15:57:53 +0000 Subject: [PATCH 3/7] Committing license headers --- include/macis/model/hubbard.hpp | 8 ++++++++ src/macis/model/hubbard.cxx | 8 ++++++++ tests/hubbard.cxx | 8 ++++++++ 3 files changed, 24 insertions(+) diff --git a/include/macis/model/hubbard.hpp b/include/macis/model/hubbard.hpp index 4fd849db..f03f03bb 100644 --- a/include/macis/model/hubbard.hpp +++ b/include/macis/model/hubbard.hpp @@ -1,3 +1,11 @@ +/* + * 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 diff --git a/src/macis/model/hubbard.cxx b/src/macis/model/hubbard.cxx index cdbfa2d6..44f62824 100644 --- a/src/macis/model/hubbard.cxx +++ b/src/macis/model/hubbard.cxx @@ -1,3 +1,11 @@ +/* + * 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 { diff --git a/tests/hubbard.cxx b/tests/hubbard.cxx index 1f62f1a9..de01a4bd 100644 --- a/tests/hubbard.cxx +++ b/tests/hubbard.cxx @@ -1,3 +1,11 @@ +/* + * 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 From 65f56b3e3b4389a0a11c75a5d9d71ac0337d5612 Mon Sep 17 00:00:00 2001 From: David Williams-Young Date: Tue, 13 Jun 2023 10:02:42 -0700 Subject: [PATCH 4/7] Add PBC --- include/macis/model/hubbard.hpp | 5 +++-- src/macis/model/hubbard.cxx | 11 +++++++++-- tests/hubbard.cxx | 20 ++++++++++++++++++-- 3 files changed, 30 insertions(+), 6 deletions(-) diff --git a/include/macis/model/hubbard.hpp b/include/macis/model/hubbard.hpp index f03f03bb..6c723ab3 100644 --- a/include/macis/model/hubbard.hpp +++ b/include/macis/model/hubbard.hpp @@ -15,7 +15,8 @@ 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); +void hubbard_1d(size_t nsites, double t, double U, + std::vector& T, std::vector& V, + bool pbc = false); } // namespace macis diff --git a/src/macis/model/hubbard.cxx b/src/macis/model/hubbard.cxx index 44f62824..0d508af2 100644 --- a/src/macis/model/hubbard.cxx +++ b/src/macis/model/hubbard.cxx @@ -10,8 +10,9 @@ namespace macis { -void hubbard_1d(size_t nsites, double t, double U, std::vector& T, - std::vector& V) { +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); @@ -28,6 +29,12 @@ void hubbard_1d(size_t nsites, double t, double U, std::vector& 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/hubbard.cxx b/tests/hubbard.cxx index de01a4bd..05bf24c5 100644 --- a/tests/hubbard.cxx +++ b/tests/hubbard.cxx @@ -26,7 +26,18 @@ TEST_CASE("Hubbard") { std::vector T, V; SECTION("1D") { - macis::hubbard_1d(nsites, t, U, T, V); + + 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) @@ -48,7 +59,12 @@ TEST_CASE("Hubbard") { REQUIRE(mat_el == Approx(-U / 2)); else if(std::abs(p - q) == 1) REQUIRE(mat_el == -t); - else + 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); } } From e0ef7715cd790f1d67266cf7f47245406b8ccdc9 Mon Sep 17 00:00:00 2001 From: Clang Robot Date: Tue, 13 Jun 2023 17:03:22 +0000 Subject: [PATCH 5/7] Committing clang-format changes --- include/macis/model/hubbard.hpp | 5 ++--- src/macis/model/hubbard.cxx | 9 ++++----- tests/hubbard.cxx | 13 +++---------- 3 files changed, 9 insertions(+), 18 deletions(-) diff --git a/include/macis/model/hubbard.hpp b/include/macis/model/hubbard.hpp index 6c723ab3..5e80fd0f 100644 --- a/include/macis/model/hubbard.hpp +++ b/include/macis/model/hubbard.hpp @@ -15,8 +15,7 @@ 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); +void hubbard_1d(size_t nsites, double t, double U, std::vector& T, + std::vector& V, bool pbc = false); } // namespace macis diff --git a/src/macis/model/hubbard.cxx b/src/macis/model/hubbard.cxx index 0d508af2..5e74b0d3 100644 --- a/src/macis/model/hubbard.cxx +++ b/src/macis/model/hubbard.cxx @@ -10,9 +10,8 @@ namespace macis { -void hubbard_1d(size_t nsites, double t, double U, - std::vector& T, std::vector& V, - bool pbc) { +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); @@ -32,8 +31,8 @@ void hubbard_1d(size_t nsites, double t, double U, // PBC for 1-D if(pbc) { - T[ (nsites-1) ] = -t; - T[ (nsites-1) * nsites] = -t; + T[(nsites - 1)] = -t; + T[(nsites - 1) * nsites] = -t; } } diff --git a/tests/hubbard.cxx b/tests/hubbard.cxx index 05bf24c5..6f3d3553 100644 --- a/tests/hubbard.cxx +++ b/tests/hubbard.cxx @@ -6,8 +6,6 @@ * See LICENSE.txt for details */ - - #include #include @@ -26,16 +24,11 @@ TEST_CASE("Hubbard") { std::vector T, V; SECTION("1D") { - bool pbc; - SECTION("No PBC") { - pbc = false; - } + SECTION("No PBC") { pbc = false; } - SECTION("PBC") { - pbc = true; - } + SECTION("PBC") { pbc = true; } macis::hubbard_1d(nsites, t, U, T, V, pbc); @@ -60,7 +53,7 @@ TEST_CASE("Hubbard") { 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) ) + if((p == 0 and q == nsites - 1) or (p == nsites - 1 and q == 0)) REQUIRE(mat_el == -t); else REQUIRE(mat_el == 0.0); From d56989c0cc112585ee58407b7db656f2348df387 Mon Sep 17 00:00:00 2001 From: Carlos Mejuto Zaera Date: Tue, 19 Sep 2023 16:40:53 +0200 Subject: [PATCH 6/7] Included new HamiltonianGenerator, particularly efficient for impurity models. Test case added. --- include/macis/bitset_operations.hpp | 3 + include/macis/hamiltonian_generator.hpp | 2 +- .../macis/hamiltonian_generator/sd_build.hpp | 274 ++++++++++++++++++ include/macis/util/fcidump.hpp | 11 + src/macis/fcidump.cxx | 28 ++ tests/CMakeLists.txt | 1 + tests/ref_data/hubbard10.fci.dat | 38 +++ tests/sd_build.cxx | 115 ++++++++ tests/ut_common.hpp.in | 2 + 9 files changed, 473 insertions(+), 1 deletion(-) create mode 100644 include/macis/hamiltonian_generator/sd_build.hpp create mode 100644 tests/ref_data/hubbard10.fci.dat create mode 100644 tests/sd_build.cxx 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"; From 614f5ff8e6f7b8e227c30ee87386f8b52d5a176a Mon Sep 17 00:00:00 2001 From: Clang Robot Date: Tue, 19 Sep 2023 14:41:38 +0000 Subject: [PATCH 7/7] Committing clang-format changes --- .../macis/hamiltonian_generator/sd_build.hpp | 200 +++++++++--------- include/macis/util/fcidump.hpp | 4 +- src/macis/fcidump.cxx | 10 +- tests/sd_build.cxx | 31 +-- 4 files changed, 124 insertions(+), 121 deletions(-) diff --git a/include/macis/hamiltonian_generator/sd_build.hpp b/include/macis/hamiltonian_generator/sd_build.hpp index b1170e79..37c8cf0f 100644 --- a/include/macis/hamiltonian_generator/sd_build.hpp +++ b/include/macis/hamiltonian_generator/sd_build.hpp @@ -14,17 +14,17 @@ namespace macis { -template< size_t N > -struct det_pos -{ - public: - std::bitset det; - uint32_t id; +template +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 +bool operator<(const det_pos& a, const det_pos& b) { + return bitset_less(a.det, b.det); +} template class SDBuildHamiltonianGenerator : public HamiltonianGenerator { @@ -54,94 +54,94 @@ class SDBuildHamiltonianGenerator : public HamiltonianGenerator { // 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; + 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++ ) - { + 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 ); + 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; + // 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); + 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); + 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 - { + if(just_singles) + generate_singles_spin(this->norb_, bra, excs); + else { std::vector singls; - generate_singles_spin( this->norb_, bra, excs ); + 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() ); + // 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 ); + 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 ); + const auto h_eld = + this->matrix_element_diag(bra_occ_alpha, bra_occ_beta); - if( std::abs(h_eld) > H_thresh ) { + 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() ) { + 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 ); + 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 ); + + 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 ); + 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 ) { + 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 ket determinant + } // Loop over ket determinants + } // Non-zero bra determinant rowptr[i + 1] = rowptr[i] + nrow; // Update rowptr @@ -174,19 +174,18 @@ class SDBuildHamiltonianGenerator : public HamiltonianGenerator { 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 { + 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++ ) - { + 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 ); + a.id = std::distance(ket_begin, it); + kets.insert(a); } // Loop over bra determinants @@ -196,79 +195,84 @@ class SDBuildHamiltonianGenerator : public HamiltonianGenerator { 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); + 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() ){ + 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 ); - } + 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 ); + 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 ); + 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() ) { + 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 ); + 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 ); + 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 + 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_); } + 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; } + 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; } + bool GetJustSingles() const override { return just_singles; } }; } // namespace macis diff --git a/include/macis/util/fcidump.hpp b/include/macis/util/fcidump.hpp index 6821f434..823ef7c8 100644 --- a/include/macis/util/fcidump.hpp +++ b/include/macis/util/fcidump.hpp @@ -71,8 +71,8 @@ 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. + * @brief Check whether the 2-body contribution of the Hamiltonian is + * exclusively diagonal. * * @param[in] fname: Filename of FCIDUMP file * diff --git a/src/macis/fcidump.cxx b/src/macis/fcidump.cxx index 6af87371..9870077c 100644 --- a/src/macis/fcidump.cxx +++ b/src/macis/fcidump.cxx @@ -182,8 +182,7 @@ void read_fcidump_2body(std::string fname, double* V, size_t LDV) { } bool is_2body_diagonal(std::string fname) { - - auto norb = read_fcidump_norb(fname); + auto norb = read_fcidump_norb(fname); bool all_diag = true; std::ifstream file(fname); @@ -199,10 +198,9 @@ bool is_2body_diagonal(std::string fname) { q--; r--; s--; - if(p != q || r != s) - { - all_diag = false; - break; + if(p != q || r != s) { + all_diag = false; + break; } } } diff --git a/tests/sd_build.cxx b/tests/sd_build.cxx index 0a8c95cb..9c2a0760 100644 --- a/tests/sd_build.cxx +++ b/tests/sd_build.cxx @@ -9,16 +9,16 @@ #include #include #include -#include -#include #include +#include #include +#include #include "ut_common.hpp" -using macis::NumOrbital; using macis::NumActive; using macis::NumInactive; +using macis::NumOrbital; TEST_CASE("Single Double Build") { ROOT_ONLY(MPI_COMM_WORLD); @@ -33,7 +33,7 @@ TEST_CASE("Single Double Build") { 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 ); + bool just_singles = macis::is_2body_diagonal(hubbard10_fcidump); using generator_type = macis::SDBuildHamiltonianGenerator<64>; @@ -78,11 +78,11 @@ TEST_CASE("Single Double Build") { } SECTION("CI") { - size_t nalpha = 5; - size_t nbeta = 5; + size_t nalpha = 5; + size_t nbeta = 5; size_t n_inactive = 0; - size_t n_active = 10; - size_t n_virtual = 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; @@ -101,15 +101,16 @@ TEST_CASE("Single Double Build") { 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); + 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)); + 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)); } }