diff --git a/INSTALL.md b/INSTALL.md index 842a7d2c6..ec1c5e26a 100644 --- a/INSTALL.md +++ b/INSTALL.md @@ -5,7 +5,7 @@ prerequisites: * mandatory: * CMake 3.15 or later * a C++17 compiler - * [Boost](https://www.boost.org/), version 1.67 or higher (N.B. critical bugs make the following versions unusable: 1.70, 1.77, 1.78); the following non-header-only Boost libraries are required: + * [Boost](https://www.boost.org/), version 1.81 or higher (N.B. older compilers _may_ work with older Boost releases). *SeQuant can download and build Boost if configured with `Boost_FETCH_IF_MISSING=ON`, but this is not recommended.* The following non-header-only Boost libraries are required, hence Boost must be configured/built: - [Boost.Regex](https://www.boost.org/doc/libs/master/libs/regex/doc/html/index.html) - [Boost.Locale](https://www.boost.org/doc/libs/master/libs/locale/doc/html/index.html) * [Range-V3](https://github.com/ericniebler/range-v3.git), tag 0.12.0, *if not found, SeQuant will download and build Range-V3* @@ -20,8 +20,9 @@ for the impatient (from the top of the SeQuant source directory): * `cmake --build build --target check-sequant` useful CMake variables: - * `BUILD_TESTING` --- enables unit tests targets, e.g. `check-sequant` [default=ON] - * `CMAKE_CXX_COMPILER` --- specifies the C++ compiler to use - * `CMAKE_PREFIX_PATH` --- this semicolon-separated list specifies search paths for dependencies (Boost, Range-V3, etc.) + * [`BUILD_TESTING`](https://cmake.org/cmake/help/latest/module/CTest.html) --- enables unit tests targets, e.g. `check-sequant` [default=ON] + * [`CMAKE_CXX_COMPILER`](https://cmake.org/cmake/help/latest/variable/CMAKE_LANG_COMPILER.html#variable:CMAKE_%3CLANG%3E_COMPILER) --- specifies the C++ compiler to use + * [`CMAKE_PREFIX_PATH`](https://cmake.org/cmake/help/latest/variable/CMAKE_PREFIX_PATH.html) --- this semicolon-separated list specifies search paths for dependencies (Boost, Range-V3, etc.) * `SEQUANT_MIMALLOC` --- use [mimalloc](https://github.com/microsoft/mimalloc) for fast memory allocation * `SEQUANT_EVAL_TRACE` --- enables tracing of expression interpretation; especially useful in combination with TiledArray's memory tracing mechanism (configure TiledArray with `TA_TENSOR_MEM_PROFILE=ON` to enable that) + * `Boost_FETCH_IF_MISSING` --- if set to `ON`, SeQuant will download and build Boost if it is not found by `find_package(Boost ...)`; this is not recommended. [default=OFF] diff --git a/SeQuant/core/hash.hpp b/SeQuant/core/hash.hpp index f71cea284..b7ef212fb 100644 --- a/SeQuant/core/hash.hpp +++ b/SeQuant/core/hash.hpp @@ -14,6 +14,10 @@ namespace sequant_boost = boost; #include #endif +#if SEQUANT_BOOST_VERSION < 108100 +#error "SeQuant requires Boost 1.81 or later for hashing" +#endif + #include "meta.hpp" namespace sequant { @@ -27,17 +31,7 @@ enum class Impl { BoostPre181 = 1, Boost181OrLater = 2 }; /// @return the version of hashing used by SeQuant, depends on the version of /// Boost -constexpr hash::Impl hash_version() { -#ifdef SEQUANT_USE_SYSTEM_BOOST_HASH -#if SEQUANT_BOOST_VERSION < 108100 - return hash::Impl::BoostPre181; -#else - return hash::Impl::Boost181OrLater; -#endif -#else - return hash::Impl::Boost181OrLater; -#endif -} +constexpr hash::Impl hash_version() { return hash::Impl::Boost181OrLater; } namespace detail { template @@ -116,45 +110,7 @@ template inline void combine(std::size_t& seed, T const& v) { _ hasher; -#ifdef SEQUANT_USE_SYSTEM_BOOST_HASH -#if SEQUANT_BOOST_VERSION >= 108100 - boost::hash_combine(seed, hasher(v)); -#else // older boost workarounds - // in boost 1.78 hash_combine_impl implementation changed - // https://github.com/boostorg/container_hash/commit/21f2b5e1db1a118c83a3690055c110d0f5637da3 - // probably no longer need these acrobatics - if constexpr (sizeof(std::size_t) == sizeof(boost::uint32_t) && - sizeof(decltype(hasher(v))) == sizeof(boost::uint32_t)) { - const boost::uint32_t value = hasher(v); -#if SEQUANT_BOOST_VERSION >= 107800 - seed = boost::hash_detail::hash_combine_impl<32>::fn( - static_cast(seed), value); -#else - // N.B. seed passed by reference - boost::hash_detail::hash_combine_impl( - reinterpret_cast(seed), value); -#endif - return; - } else if constexpr (sizeof(std::size_t) == sizeof(boost::uint64_t) && - sizeof(decltype(hasher(v))) == sizeof(boost::uint64_t)) { - const boost::uint64_t value = hasher(v); - -#if SEQUANT_BOOST_VERSION >= 107800 - seed = boost::hash_detail::hash_combine_impl<64>::fn( - static_cast(seed), value); -#else - // N.B. seed passed by reference - boost::hash_detail::hash_combine_impl( - reinterpret_cast(seed), value); -#endif - return; - } else { - seed ^= hasher(v) + 0x9e3779b9 + (seed << 6) + (seed >> 2); - } -#endif // older boost workarounds -#else // !defined(SEQUANT_USE_SYSTEM_BOOST_HASH) sequant_boost::hash_combine(seed, hasher(v)); -#endif // assert(seed == seed_ref); } diff --git a/SeQuant/domain/eval/eval.hpp b/SeQuant/domain/eval/eval.hpp index 008dbf520..9e5358b86 100644 --- a/SeQuant/domain/eval/eval.hpp +++ b/SeQuant/domain/eval/eval.hpp @@ -32,7 +32,7 @@ void log_eval(Args const&... args) noexcept { #endif } -void log_cache_access(size_t key, CacheManager const& cm) { +[[maybe_unused]] void log_cache_access(size_t key, CacheManager const& cm) { #ifdef SEQUANT_EVAL_TRACE auto& l = Logger::get_instance(); if (l.log_level_eval > 0) { @@ -50,7 +50,7 @@ void log_cache_access(size_t key, CacheManager const& cm) { #endif } -void log_cache_store(size_t key, CacheManager const& cm) { +[[maybe_unused]] void log_cache_store(size_t key, CacheManager const& cm) { #ifdef SEQUANT_EVAL_TRACE auto& l = Logger::get_instance(); if (l.log_level_eval > 0) { diff --git a/SeQuant/domain/eval/eval_result.hpp b/SeQuant/domain/eval/eval_result.hpp index a59afac06..368a1d9a4 100644 --- a/SeQuant/domain/eval/eval_result.hpp +++ b/SeQuant/domain/eval/eval_result.hpp @@ -487,7 +487,7 @@ inline void log_constant(Args const&... args) noexcept { void log_ta_tensor_host_memory_use(madness::World& world, std::string_view label = ""); -struct EvalResult; +class EvalResult; using ERPtr = std::shared_ptr; diff --git a/tests/unit/test_mbpt.cpp b/tests/unit/test_mbpt.cpp index da2279fd8..91c6847cc 100644 --- a/tests/unit/test_mbpt.cpp +++ b/tests/unit/test_mbpt.cpp @@ -281,30 +281,18 @@ TEST_CASE("NBodyOp", "[mbpt]") { auto t = t1 + t2; - if constexpr (hash_version() == hash::Impl::BoostPre181) { - REQUIRE( - to_latex(simplify(f * t * t)) == - to_latex(f * t1 * t1 + f * t2 * t2 + ex(2) * f * t1 * t2)); - } else { - // std::wcout << "to_latex(simplify(f * t * t)): " - // << to_latex(simplify(f * t * t)) << std::endl; - REQUIRE( - to_latex(simplify(f * t * t)) == - to_latex(ex(2) * f * t1 * t2 + f * t2 * t2 + f * t1 * t1)); - } - - if constexpr (hash_version() == hash::Impl::BoostPre181) { - REQUIRE(to_latex(simplify(f * t * t * t)) == - to_latex(ex(3) * f * t1 * t2 * t2 + f * t2 * t2 * t2 + - ex(3) * f * t1 * t1 * t2 + f * t1 * t1 * t1)); - } else { - // std::wcout << "to_latex(simplify(f * t * t * t): " - // << to_latex(simplify(f * t * t * t)) << std::endl; - REQUIRE(to_latex(simplify(f * t * t * t)) == - to_latex(f * t2 * t2 * t2 + f * t1 * t1 * t1 + - ex(3) * f * t1 * t1 * t2 + - ex(3) * f * t1 * t2 * t2)); - } + // std::wcout << "to_latex(simplify(f * t * t)): " + // << to_latex(simplify(f * t * t)) << std::endl; + REQUIRE( + to_latex(simplify(f * t * t)) == + to_latex(ex(2) * f * t1 * t2 + f * t2 * t2 + f * t1 * t1)); + + // std::wcout << "to_latex(simplify(f * t * t * t): " + // << to_latex(simplify(f * t * t * t)) << std::endl; + REQUIRE(to_latex(simplify(f * t * t * t)) == + to_latex(f * t2 * t2 * t2 + f * t1 * t1 * t1 + + ex(3) * f * t1 * t1 * t2 + + ex(3) * f * t1 * t2 * t2)); } // SECTION("canonicalize") diff --git a/tests/unit/test_wick.cpp b/tests/unit/test_wick.cpp index caf7241c2..a782ab5a5 100644 --- a/tests/unit/test_wick.cpp +++ b/tests/unit/test_wick.cpp @@ -148,186 +148,196 @@ TEST_CASE("WickTheorem", "[algorithms][wick]") { }); // number operator - {{auto opseq1 = - FNOperatorSeq({FNOperator({L"i_1"}, {}), FNOperator({}, {L"i_2"})}); - auto wick1 = FWickTheorem{opseq1}; - REQUIRE_NOTHROW(wick1.compute()); - // full contractions = null (N is already in normal form) - auto full_contractions = FWickTheorem{opseq1}.compute(); - REQUIRE(full_contractions->is()); - REQUIRE(full_contractions->as().value() == 0); - // partial contractions = N - auto partial_contractions = - FWickTheorem{opseq1}.full_contractions(false).compute(); - // std::wcout << "partial_contractions=" << to_latex(partial_contractions) - // << std::endl; - REQUIRE(partial_contractions->is()); - REQUIRE(partial_contractions->as().size() == 1); - } - { - auto opseq1 = - BNOperatorSeq({BNOperator({L"i_1"}, {}), BNOperator({}, {L"i_2"})}); - auto wick1 = BWickTheorem{opseq1}; - REQUIRE_NOTHROW(wick1.compute()); - // full contractions = null - auto full_contractions = BWickTheorem{opseq1}.compute(); - REQUIRE(full_contractions->is()); - REQUIRE(full_contractions->as().value() == 0); - // partial contractions = N - auto partial_contractions = - BWickTheorem{opseq1}.full_contractions(false).compute(); - // std::wcout << "partial_contractions=" << to_latex(partial_contractions) - // << std::endl; - REQUIRE(partial_contractions->is()); - REQUIRE(partial_contractions->as().size() == 1); - } -} + { + { + auto opseq1 = + FNOperatorSeq({FNOperator({L"i_1"}, {}), FNOperator({}, {L"i_2"})}); + auto wick1 = FWickTheorem{opseq1}; + REQUIRE_NOTHROW(wick1.compute()); + // full contractions = null (N is already in normal form) + auto full_contractions = FWickTheorem{opseq1}.compute(); + REQUIRE(full_contractions->is()); + REQUIRE(full_contractions->as().value() == 0); + // partial contractions = N + auto partial_contractions = + FWickTheorem{opseq1}.full_contractions(false).compute(); + // std::wcout << "partial_contractions=" << + // to_latex(partial_contractions) + // << std::endl; + REQUIRE(partial_contractions->is()); + REQUIRE(partial_contractions->as().size() == 1); + } + { + auto opseq1 = + BNOperatorSeq({BNOperator({L"i_1"}, {}), BNOperator({}, {L"i_2"})}); + auto wick1 = BWickTheorem{opseq1}; + REQUIRE_NOTHROW(wick1.compute()); + // full contractions = null + auto full_contractions = BWickTheorem{opseq1}.compute(); + REQUIRE(full_contractions->is()); + REQUIRE(full_contractions->as().value() == 0); + // partial contractions = N + auto partial_contractions = + BWickTheorem{opseq1}.full_contractions(false).compute(); + // std::wcout << "partial_contractions=" << + // to_latex(partial_contractions) + // << std::endl; + REQUIRE(partial_contractions->is()); + REQUIRE(partial_contractions->as().size() == 1); + } + } -// hole number operator -{{auto opseq1 = - FNOperatorSeq({FNOperator({}, {L"i_1"}), FNOperator({L"i_2"}, {})}); -auto wick1 = FWickTheorem{opseq1}; -REQUIRE_NOTHROW(wick1.compute()); -// full contractions = delta -auto full_contractions = FWickTheorem{opseq1}.compute(); -REQUIRE(full_contractions->is()); -REQUIRE(full_contractions->as().size() == 1); -// partial contractions = delta - N -auto partial_contractions = - FWickTheorem{opseq1}.full_contractions(false).compute(); -// std::wcout << "partial_contractions=" << to_latex(partial_contractions) << -// std::endl; -REQUIRE(partial_contractions->is()); -REQUIRE(partial_contractions->as().size() == 2); -REQUIRE(to_latex(partial_contractions) == - L"{ \\bigl({{s^{{i_2}}_{{i_1}}}} - {{a^{{i_2}}_{{i_1}}}}\\bigr) }"); -} -{ - auto opseq1 = - BNOperatorSeq({BNOperator({}, {L"i_1"}), BNOperator({L"i_2"}, {})}); - auto wick1 = BWickTheorem{opseq1}; - REQUIRE_NOTHROW(wick1.compute()); - // full contractions = delta - auto full_contractions = BWickTheorem{opseq1}.compute(); - REQUIRE(full_contractions->is()); - REQUIRE(full_contractions->as().size() == 1); - // partial contractions = delta + N - auto partial_contractions = - BWickTheorem{opseq1}.full_contractions(false).compute(); - // std::wcout << "partial_contractions=" << to_latex(partial_contractions) << - // std::endl; - REQUIRE(partial_contractions->is()); - REQUIRE(partial_contractions->as().size() == 2); - REQUIRE(to_latex(partial_contractions) == - L"{ \\bigl({{s^{{i_2}}_{{i_1}}}} + {{b^{{i_2}}_{{i_1}}}}\\bigr) }"); -} -} + // hole number operator + { + { + auto opseq1 = + FNOperatorSeq({FNOperator({}, {L"i_1"}), FNOperator({L"i_2"}, {})}); + auto wick1 = FWickTheorem{opseq1}; + REQUIRE_NOTHROW(wick1.compute()); + // full contractions = delta + auto full_contractions = FWickTheorem{opseq1}.compute(); + REQUIRE(full_contractions->is()); + REQUIRE(full_contractions->as().size() == 1); + // partial contractions = delta - N + auto partial_contractions = + FWickTheorem{opseq1}.full_contractions(false).compute(); + // std::wcout << "partial_contractions=" << + // to_latex(partial_contractions) << std::endl; + REQUIRE(partial_contractions->is()); + REQUIRE(partial_contractions->as().size() == 2); + REQUIRE( + to_latex(partial_contractions) == + L"{ \\bigl({{s^{{i_2}}_{{i_1}}}} - {{a^{{i_2}}_{{i_1}}}}\\bigr) }"); + } + { + auto opseq1 = + BNOperatorSeq({BNOperator({}, {L"i_1"}), BNOperator({L"i_2"}, {})}); + auto wick1 = BWickTheorem{opseq1}; + REQUIRE_NOTHROW(wick1.compute()); + // full contractions = delta + auto full_contractions = BWickTheorem{opseq1}.compute(); + REQUIRE(full_contractions->is()); + REQUIRE(full_contractions->as().size() == 1); + // partial contractions = delta + N + auto partial_contractions = + BWickTheorem{opseq1}.full_contractions(false).compute(); + // std::wcout << "partial_contractions=" << + // to_latex(partial_contractions) << std::endl; + REQUIRE(partial_contractions->is()); + REQUIRE(partial_contractions->as().size() == 2); + REQUIRE( + to_latex(partial_contractions) == + L"{ \\bigl({{s^{{i_2}}_{{i_1}}}} + {{b^{{i_2}}_{{i_1}}}}\\bigr) }"); + } + } -// three 1-body operators -{ - auto opseq1 = FNOperatorSeq({FNOperator({L"i_1"}, {L"i_2"}), - FNOperator({L"i_3"}, {L"i_4"}), - FNOperator({L"i_5"}, {L"i_6"})}); - auto wick1 = FWickTheorem{opseq1}; - REQUIRE_NOTHROW(wick1.compute()); - auto result = FWickTheorem{opseq1}.compute(); - REQUIRE(result->is()); - REQUIRE(result->as().value() == 0); -} + // three 1-body operators + { + auto opseq1 = FNOperatorSeq({FNOperator({L"i_1"}, {L"i_2"}), + FNOperator({L"i_3"}, {L"i_4"}), + FNOperator({L"i_5"}, {L"i_6"})}); + auto wick1 = FWickTheorem{opseq1}; + REQUIRE_NOTHROW(wick1.compute()); + auto result = FWickTheorem{opseq1}.compute(); + REQUIRE(result->is()); + REQUIRE(result->as().value() == 0); + } -// two 2-body operators -{ - auto opseq = FNOperatorSeq( - {FNOperator({}, {L"i_1", L"i_2"}), FNOperator({L"i_3", L"i_4"}, {})}); - auto wick = FWickTheorem{opseq}; - REQUIRE_NOTHROW(wick.compute()); - auto result = wick.compute(); - REQUIRE(result->is()); - REQUIRE(result->size() == 2); -} + // two 2-body operators + { + auto opseq = FNOperatorSeq( + {FNOperator({}, {L"i_1", L"i_2"}), FNOperator({L"i_3", L"i_4"}, {})}); + auto wick = FWickTheorem{opseq}; + REQUIRE_NOTHROW(wick.compute()); + auto result = wick.compute(); + REQUIRE(result->is()); + REQUIRE(result->size() == 2); + } -// two 3-body operators -{ - auto opseq = FNOperatorSeq({FNOperator({}, {L"i_1", L"i_2", L"i_3"}), - FNOperator({L"i_4", L"i_5", L"i_6"}, {})}); - auto wick = FWickTheorem{opseq}; - REQUIRE_NOTHROW(wick.compute()); - auto result = wick.compute(); - REQUIRE(result->is()); - REQUIRE(result->size() == 6); -} + // two 3-body operators + { + auto opseq = FNOperatorSeq({FNOperator({}, {L"i_1", L"i_2", L"i_3"}), + FNOperator({L"i_4", L"i_5", L"i_6"}, {})}); + auto wick = FWickTheorem{opseq}; + REQUIRE_NOTHROW(wick.compute()); + auto result = wick.compute(); + REQUIRE(result->is()); + REQUIRE(result->size() == 6); + } -// two 4-body operators -{ - auto opseq = - FNOperatorSeq({FNOperator({}, {L"i_1", L"i_2", L"i_3", L"i_4"}), - FNOperator({L"i_5", L"i_6", L"i_7", L"i_8"}, {})}); - auto wick = FWickTheorem{opseq}; - REQUIRE_NOTHROW(wick.compute()); - auto result = wick.compute(); - REQUIRE(result->is()); - REQUIRE(result->size() == 24); -} + // two 4-body operators + { + auto opseq = + FNOperatorSeq({FNOperator({}, {L"i_1", L"i_2", L"i_3", L"i_4"}), + FNOperator({L"i_5", L"i_6", L"i_7", L"i_8"}, {})}); + auto wick = FWickTheorem{opseq}; + REQUIRE_NOTHROW(wick.compute()); + auto result = wick.compute(); + REQUIRE(result->is()); + REQUIRE(result->size() == 24); + } -// 1/2 * 1 * 1/2 body ops, full contraction -{ - auto opseq = - FNOperatorSeq({FNOperator({}, {L"i_1"}), FNOperator({L"i_2"}, {L"i_3"}), - FNOperator({L"i_4"}, {})}); - auto wick = FWickTheorem{opseq}; - REQUIRE_NOTHROW(wick.compute()); - auto result = wick.compute(); - REQUIRE(result->is()); - REQUIRE(result->size() == 2); // product of 2 terms -} + // 1/2 * 1 * 1/2 body ops, full contraction + { + auto opseq = FNOperatorSeq({FNOperator({}, {L"i_1"}), + FNOperator({L"i_2"}, {L"i_3"}), + FNOperator({L"i_4"}, {})}); + auto wick = FWickTheorem{opseq}; + REQUIRE_NOTHROW(wick.compute()); + auto result = wick.compute(); + REQUIRE(result->is()); + REQUIRE(result->size() == 2); // product of 2 terms + } -// 1/2 * 1 * 1/2 body ops, partial contraction -{ - auto opseq = - FNOperatorSeq({FNOperator({}, {L"i_1"}), FNOperator({L"i_2"}, {L"i_3"}), - FNOperator({L"i_4"}, {})}); - auto wick = FWickTheorem{opseq}; - REQUIRE_NOTHROW(wick.full_contractions(false).compute()); - auto result = wick.full_contractions(false).compute(); - REQUIRE(result->is()); - REQUIRE(result->size() == 5); // sum of 4 terms - REQUIRE(to_latex(result) == - L"{ \\bigl( - {{s^{{i_2}}_{{i_1}}}{a^{{i_4}}_{{i_3}}}} + " - L"{{s^{{i_2}}_{{i_1}}}{s^{{i_4}}_{{i_3}}}} + " - L"{{s^{{i_4}}_{{i_1}}}{a^{{i_2}}_{{i_3}}}} - " - L"{{s^{{i_4}}_{{i_3}}}{a^{{i_2}}_{{i_1}}}} - " - L"{{a^{{i_2}{i_4}}_{{i_3}{i_1}}}}\\bigr) }"); -} + // 1/2 * 1 * 1/2 body ops, partial contraction + { + auto opseq = FNOperatorSeq({FNOperator({}, {L"i_1"}), + FNOperator({L"i_2"}, {L"i_3"}), + FNOperator({L"i_4"}, {})}); + auto wick = FWickTheorem{opseq}; + REQUIRE_NOTHROW(wick.full_contractions(false).compute()); + auto result = wick.full_contractions(false).compute(); + REQUIRE(result->is()); + REQUIRE(result->size() == 5); // sum of 4 terms + REQUIRE(to_latex(result) == + L"{ \\bigl( - {{s^{{i_2}}_{{i_1}}}{a^{{i_4}}_{{i_3}}}} + " + L"{{s^{{i_2}}_{{i_1}}}{s^{{i_4}}_{{i_3}}}} + " + L"{{s^{{i_4}}_{{i_1}}}{a^{{i_2}}_{{i_3}}}} - " + L"{{s^{{i_4}}_{{i_3}}}{a^{{i_2}}_{{i_1}}}} - " + L"{{a^{{i_2}{i_4}}_{{i_3}{i_1}}}}\\bigr) }"); + } -// three 1-body operators, partial contraction -{ - auto opseq1 = FNOperatorSeq({FNOperator({L"i_1"}, {L"i_2"}), - FNOperator({L"i_3"}, {L"i_4"}), - FNOperator({L"i_5"}, {L"i_6"})}); - auto wick1 = FWickTheorem{opseq1}; - REQUIRE_NOTHROW(wick1.full_contractions(false).compute()); - auto result = FWickTheorem{opseq1}.full_contractions(false).compute(); - REQUIRE(result->is()); - REQUIRE(result->size() == 5); - REQUIRE(to_latex(result) == - L"{ \\bigl({{s^{{i_3}}_{{i_2}}}{a^{{i_1}{i_5}}_{{i_4}{i_6}}}} + " - L"{{s^{{i_3}}_{{i_2}}}{s^{{i_5}}_{{i_4}}}{a^{{i_1}}_{{i_6}}}} + " - L"{{s^{{i_5}}_{{i_2}}}{a^{{i_1}{i_3}}_{{i_6}{i_4}}}} + " - L"{{s^{{i_5}}_{{i_4}}}{a^{{i_1}{i_3}}_{{i_2}{i_6}}}} + " - L"{{a^{{i_1}{i_3}{i_5}}_{{i_2}{i_4}{i_6}}}}\\bigr) }"); -} + // three 1-body operators, partial contraction + { + auto opseq1 = FNOperatorSeq({FNOperator({L"i_1"}, {L"i_2"}), + FNOperator({L"i_3"}, {L"i_4"}), + FNOperator({L"i_5"}, {L"i_6"})}); + auto wick1 = FWickTheorem{opseq1}; + REQUIRE_NOTHROW(wick1.full_contractions(false).compute()); + auto result = FWickTheorem{opseq1}.full_contractions(false).compute(); + REQUIRE(result->is()); + REQUIRE(result->size() == 5); + REQUIRE(to_latex(result) == + L"{ \\bigl({{s^{{i_3}}_{{i_2}}}{a^{{i_1}{i_5}}_{{i_4}{i_6}}}} + " + L"{{s^{{i_3}}_{{i_2}}}{s^{{i_5}}_{{i_4}}}{a^{{i_1}}_{{i_6}}}} + " + L"{{s^{{i_5}}_{{i_2}}}{a^{{i_1}{i_3}}_{{i_6}{i_4}}}} + " + L"{{s^{{i_5}}_{{i_4}}}{a^{{i_1}{i_3}}_{{i_2}{i_6}}}} + " + L"{{a^{{i_1}{i_3}{i_5}}_{{i_2}{i_4}{i_6}}}}\\bigr) }"); + } -// two 2-body operators, partial contraction: Eq. 9b of DOI 10.1063/1.474405 -{ - auto opseq = FNOperatorSeq({FNOperator({L"i_1", L"i_2"}, {L"i_3", L"i_4"}), - FNOperator({L"i_5", L"i_6"}, {L"i_7", L"i_8"})}); - auto wick = FWickTheorem{opseq}; - REQUIRE_NOTHROW(wick.full_contractions(false).compute()); - auto result = wick.full_contractions(false).compute(); - auto result_latex = to_latex(result); - REQUIRE(result->is()); - REQUIRE(result->size() == 7); - REQUIRE(result_latex == + // two 2-body operators, partial contraction: Eq. 9b of DOI 10.1063/1.474405 + { + auto opseq = + FNOperatorSeq({FNOperator({L"i_1", L"i_2"}, {L"i_3", L"i_4"}), + FNOperator({L"i_5", L"i_6"}, {L"i_7", L"i_8"})}); + auto wick = FWickTheorem{opseq}; + REQUIRE_NOTHROW(wick.full_contractions(false).compute()); + auto result = wick.full_contractions(false).compute(); + auto result_latex = to_latex(result); + REQUIRE(result->is()); + REQUIRE(result->size() == 7); + REQUIRE( + result_latex == L"{ " L"\\bigl({{s^{{i_5}}_{{i_4}}}{a^{{i_1}{i_2}{i_6}}_{{i_3}{i_7}{i_8}}}}" L" + " @@ -339,256 +349,263 @@ REQUIRE(to_latex(partial_contractions) == L"+ {{s^{{i_6}}_{{i_3}}}{a^{{i_1}{i_2}{i_5}}_{{i_8}{i_4}{i_7}}}} + " L"{{a^{{i_1}{i_2}{i_5}{i_6}}_{{i_3}{i_4}{i_7}{i_8}}}}\\bigr) }"); - // if Wick's theorem's result is in "canonical" (columns-matching-inputs ... - // this is what Kutzelnigg calls generalized Wick's theorem) it works same - // for spinorbital and spinfree basis for physical vacuum - auto raii_tmp = switch_to_spinfree_context(); - REQUIRE_NOTHROW(wick.full_contractions(false).compute()); - auto result_sf = wick.full_contractions(false).compute(); - auto result_sf_latex = to_latex(result_sf); - // std::wcout << "result_sf: " << to_latex(result_sf) << std::endl; - REQUIRE(result->is()); - REQUIRE(result->size() == 7); - auto result_sf_latex_with_E_replaced_by_a = - result_sf_latex | ranges::views::replace(L'E', L'a') | - ranges::to(); - REQUIRE(result_latex == result_sf_latex_with_E_replaced_by_a); -} + // if Wick's theorem's result is in "canonical" (columns-matching-inputs + // ... this is what Kutzelnigg calls generalized Wick's theorem) it works + // same for spinorbital and spinfree basis for physical vacuum + auto raii_tmp = switch_to_spinfree_context(); + REQUIRE_NOTHROW(wick.full_contractions(false).compute()); + auto result_sf = wick.full_contractions(false).compute(); + auto result_sf_latex = to_latex(result_sf); + // std::wcout << "result_sf: " << to_latex(result_sf) << std::endl; + REQUIRE(result->is()); + REQUIRE(result->size() == 7); + auto result_sf_latex_with_E_replaced_by_a = + result_sf_latex | ranges::views::replace(L'E', L'a') | + ranges::to(); + REQUIRE(result_latex == result_sf_latex_with_E_replaced_by_a); + } -} // SECTION("physical vacuum") + } // SECTION("physical vacuum") -SECTION("fermi vacuum") { - // default vacuum is already spin-orbital Fermi vacuum + SECTION("fermi vacuum") { + // default vacuum is already spin-orbital Fermi vacuum - auto switch_to_spinfree_context = detail::NoDiscard([&]() { - auto context_sf = get_default_context(); - context_sf.set(SPBasis::spinfree); - return set_scoped_default_context(context_sf); - }); + auto switch_to_spinfree_context = detail::NoDiscard([&]() { + auto context_sf = get_default_context(); + context_sf.set(SPBasis::spinfree); + return set_scoped_default_context(context_sf); + }); - // two (pure qp) 1-body operators - { - auto opseq = FNOperatorSeq( - {FNOperator({L"i_1"}, {L"a_1"}), FNOperator({L"a_2"}, {L"i_2"})}); - auto wick = FWickTheorem{opseq}; - REQUIRE_NOTHROW(wick.compute()); - auto result = wick.compute(); - REQUIRE(result->is()); - REQUIRE(result->size() == 2); // product of 2 terms - - // spin-free result is simply twice the spin-orbital result - auto raii_tmp = switch_to_spinfree_context(); - auto result_sf = wick.compute(); - REQUIRE(simplify(result_sf - result * ex(2)) == ex(0)); - } + // two (pure qp) 1-body operators + { + auto opseq = FNOperatorSeq( + {FNOperator({L"i_1"}, {L"a_1"}), FNOperator({L"a_2"}, {L"i_2"})}); + auto wick = FWickTheorem{opseq}; + REQUIRE_NOTHROW(wick.compute()); + auto result = wick.compute(); + REQUIRE(result->is()); + REQUIRE(result->size() == 2); // product of 2 terms - // two (pure qp) N-nonconserving 2-body operators - { - auto opseq = FNOperatorSeq({FNOperator({L"i_1", L"i_2"}, {L"a_1"}), - FNOperator({L"a_2"}, {L"i_3", L"i_4"})}); - auto wick = FWickTheorem{opseq}; - REQUIRE_NOTHROW(wick.compute()); - auto result = wick.compute(); - REQUIRE(result->is()); - REQUIRE(result->size() == 2); - } + // spin-free result is simply twice the spin-orbital result + auto raii_tmp = switch_to_spinfree_context(); + auto result_sf = wick.compute(); + REQUIRE(simplify(result_sf - result * ex(2)) == + ex(0)); + } - // two general 1-body operators - { - auto opseq = FNOperatorSeq( - {FNOperator({L"p_1"}, {L"p_2"}), FNOperator({L"p_3"}, {L"p_4"})}); - auto wick = FWickTheorem{opseq}; - REQUIRE_NOTHROW(wick.compute()); - auto result = wick.compute(); - REQUIRE(result->is()); - REQUIRE(result->size() == - 2 * 2); // product of 4 terms (since each contraction of 2 - // *general* indices produces 2 overlaps) - REQUIRE(to_latex(result) == - L"{{s^{{p_1}}_{{m_{102}}}}{s^{{m_{102}}}_{{p_4}}}{s^{{E_{103}}}_{" - L"{p_2}}}{s^{{p_3}}_{{E_{103}}}}}"); - } - // two general 1-body operators, partial contractions: Eq. 21a of - // DOI 10.1063/1.474405 - { - auto opseq = FNOperatorSeq( - {FNOperator({L"p_1"}, {L"p_2"}), FNOperator({L"p_3"}, {L"p_4"})}); - auto wick = FWickTheorem{opseq}; - REQUIRE_NOTHROW(wick.full_contractions(false).compute()); - auto result = wick.full_contractions(false).compute(); - REQUIRE(result->is()); - REQUIRE(result->size() == 4); - REQUIRE( - to_latex(result) == - L"{ \\bigl( - " - L"{{s^{{p_1}}_{{m_{107}}}}{s^{{m_{107}}}_{{p_4}}}{\\tilde{a}^{{p_3}}_" - L"{{p_2}}}} + " - L"{{s^{{p_1}}_{{m_{107}}}}{s^{{m_{107}}}_{{p_4}}}{s^{{E_{108}}}_{{p_" - L"2}}}{s^{{p_3}}_{{E_{108}}}}} + " - L"{{s^{{E_{109}}}_{{p_2}}}{s^{{p_3}}_{{E_{109}}}}{\\tilde{a}^{{p_1}}_" - L"{{p_4}}}} + {{\\tilde{a}^{{p_1}{p_3}}_{{p_2}{p_4}}}}\\bigr) }"); - } + // two (pure qp) N-nonconserving 2-body operators + { + auto opseq = FNOperatorSeq({FNOperator({L"i_1", L"i_2"}, {L"a_1"}), + FNOperator({L"a_2"}, {L"i_3", L"i_4"})}); + auto wick = FWickTheorem{opseq}; + REQUIRE_NOTHROW(wick.compute()); + auto result = wick.compute(); + REQUIRE(result->is()); + REQUIRE(result->size() == 2); + } - // two (pure qp) 2-body operators - { - auto opseq = - FNOperatorSeq({FNOperator({L"i_1", L"i_2"}, {L"a_1", L"a_2"}), - FNOperator({L"a_3", L"a_4"}, {L"i_3", L"i_4"})}); - auto wick = FWickTheorem{opseq}; - REQUIRE_NOTHROW(wick.compute()); - auto result = wick.compute(); - auto result_latex = to_latex(result); - // std::wcout << "<" << to_latex(opseq) << "> = " << result_latex << - // std::endl; - REQUIRE(result->is()); - REQUIRE(result->size() == 4); - REQUIRE(result_latex == - L"{ " - L"\\bigl({{s^{{i_4}}_{{i_1}}}{s^{{i_3}}_{{i_2}}}{s^{{a_3}}_{{a_2}}}" - L"{s^{{a_4}}_{{a_1}}}} - " - L"{{s^{{i_4}}_{{i_1}}}{s^{{i_3}}_{{i_2}}}{s^{{a_4}}_{{a_2}}}{s^{{a_" - L"3}}_{{a_1}}}} - " - L"{{s^{{i_3}}_{{i_1}}}{s^{{i_4}}_{{i_2}}}{s^{{a_3}}_{{a_2}}}{s^{{a_" - L"4}}_{{a_1}}}} + " - L"{{s^{{i_3}}_{{i_1}}}{s^{{i_4}}_{{i_2}}}{s^{{a_4}}_{{a_2}}}{s^{{a_" - L"3}}_{{a_1}}}}\\bigr) }"); - - // in spin-free result first and fourth terms are multiplied by 4, second - // and third terms multiplied by 2 - auto raii_tmp = switch_to_spinfree_context(); - auto result_sf = wick.compute(); - auto result_sf_latex = to_latex(result_sf); - // std::wcout << "<" << to_latex(opseq) << "> = " << to_latex(result_sf) - // << std::endl; - REQUIRE(result_sf->is()); - REQUIRE(result_sf->size() == 4); - REQUIRE(result_sf_latex == - L"{ " - L"\\bigl({{{4}}{s^{{i_4}}_{{i_1}}}{s^{{i_3}}_{{i_2}}}{s^{{a_3}}_{{" - L"a_2}}}{s^{{a_4}}_{{a_1}}}} - " - L"{{{2}}{s^{{i_4}}_{{i_1}}}{s^{{i_3}}_{{i_2}}}{s^{{a_4}}_{{a_2}}}{" - L"s^{{a_3}}_{{a_1}}}} - " - L"{{{2}}{s^{{i_3}}_{{i_1}}}{s^{{i_4}}_{{i_2}}}{s^{{a_3}}_{{a_2}}}{" - L"s^{{a_4}}_{{a_1}}}} + " - L"{{{4}}{s^{{i_3}}_{{i_1}}}{s^{{i_4}}_{{i_2}}}{s^{{a_4}}_{{a_2}}}{" - L"s^{{a_3}}_{{a_1}}}}\\bigr) }"); - } - // two (pure qp) 3-body operators - { - auto opseq = FNOperatorSeq( - {FNOperator({L"i_1", L"i_2", L"i_3"}, {L"a_1", L"a_2", L"a_3"}), - FNOperator({L"a_4", L"a_5", L"a_6"}, {L"i_4", L"i_5", L"i_6"})}); - auto wick = FWickTheorem{opseq}; - REQUIRE_NOTHROW(wick.compute()); - auto result = wick.compute(); - REQUIRE(result->is()); - REQUIRE(result->size() == 36); - } + // two general 1-body operators + { + auto opseq = FNOperatorSeq( + {FNOperator({L"p_1"}, {L"p_2"}), FNOperator({L"p_3"}, {L"p_4"})}); + auto wick = FWickTheorem{opseq}; + REQUIRE_NOTHROW(wick.compute()); + auto result = wick.compute(); + REQUIRE(result->is()); + REQUIRE(result->size() == + 2 * 2); // product of 4 terms (since each contraction of 2 + // *general* indices produces 2 overlaps) + REQUIRE(to_latex(result) == + L"{{s^{{p_1}}_{{m_{102}}}}{s^{{m_{102}}}_{{p_4}}}{s^{{E_{103}}}_{" + L"{p_2}}}{s^{{p_3}}_{{E_{103}}}}}"); + } + // two general 1-body operators, partial contractions: Eq. 21a of + // DOI 10.1063/1.474405 + { + auto opseq = FNOperatorSeq( + {FNOperator({L"p_1"}, {L"p_2"}), FNOperator({L"p_3"}, {L"p_4"})}); + auto wick = FWickTheorem{opseq}; + REQUIRE_NOTHROW(wick.full_contractions(false).compute()); + auto result = wick.full_contractions(false).compute(); + REQUIRE(result->is()); + REQUIRE(result->size() == 4); + REQUIRE( + to_latex(result) == + L"{ \\bigl( - " + L"{{s^{{p_1}}_{{m_{107}}}}{s^{{m_{107}}}_{{p_4}}}{\\tilde{a}^{{p_3}}_" + L"{{p_2}}}} + " + L"{{s^{{p_1}}_{{m_{107}}}}{s^{{m_{107}}}_{{p_4}}}{s^{{E_{108}}}_{{p_" + L"2}}}{s^{{p_3}}_{{E_{108}}}}} + " + L"{{s^{{E_{109}}}_{{p_2}}}{s^{{p_3}}_{{E_{109}}}}{\\tilde{a}^{{p_1}}_" + L"{{p_4}}}} + {{\\tilde{a}^{{p_1}{p_3}}_{{p_2}{p_4}}}}\\bigr) }"); + } - // one general 1-body operator + one general 2-body operator, partial - // contraction: Eq. 9 of DOI 10.1063/1.474405 - { - auto opseq = - FNOperatorSeq({FNOperator({L"p_1"}, {L"p_2"}), - FNOperator({L"p_3", L"p_4"}, {L"p_5", L"p_6"})}); - auto wick = FWickTheorem{opseq}; - REQUIRE_NOTHROW(wick.full_contractions(false).compute()); - auto result = wick.full_contractions(false).compute(); - REQUIRE(result->is()); - REQUIRE(result->size() == 9); - } + // two (pure qp) 2-body operators + { + auto opseq = + FNOperatorSeq({FNOperator({L"i_1", L"i_2"}, {L"a_1", L"a_2"}), + FNOperator({L"a_3", L"a_4"}, {L"i_3", L"i_4"})}); + auto wick = FWickTheorem{opseq}; + REQUIRE_NOTHROW(wick.compute()); + auto result = wick.compute(); + auto result_latex = to_latex(result); + // std::wcout << "<" << to_latex(opseq) << "> = " << result_latex << + // std::endl; + REQUIRE(result->is()); + REQUIRE(result->size() == 4); + REQUIRE( + result_latex == + L"{ " + L"\\bigl({{s^{{i_4}}_{{i_1}}}{s^{{i_3}}_{{i_2}}}{s^{{a_3}}_{{a_2}}}" + L"{s^{{a_4}}_{{a_1}}}} - " + L"{{s^{{i_4}}_{{i_1}}}{s^{{i_3}}_{{i_2}}}{s^{{a_4}}_{{a_2}}}{s^{{a_" + L"3}}_{{a_1}}}} - " + L"{{s^{{i_3}}_{{i_1}}}{s^{{i_4}}_{{i_2}}}{s^{{a_3}}_{{a_2}}}{s^{{a_" + L"4}}_{{a_1}}}} + " + L"{{s^{{i_3}}_{{i_1}}}{s^{{i_4}}_{{i_2}}}{s^{{a_4}}_{{a_2}}}{s^{{a_" + L"3}}_{{a_1}}}}\\bigr) }"); + + // in spin-free result first and fourth terms are multiplied by 4, second + // and third terms multiplied by 2 + auto raii_tmp = switch_to_spinfree_context(); + auto result_sf = wick.compute(); + auto result_sf_latex = to_latex(result_sf); + // std::wcout << "<" << to_latex(opseq) << "> = " << + // to_latex(result_sf) + // << std::endl; + REQUIRE(result_sf->is()); + REQUIRE(result_sf->size() == 4); + REQUIRE( + result_sf_latex == + L"{ " + L"\\bigl({{{4}}{s^{{i_4}}_{{i_1}}}{s^{{i_3}}_{{i_2}}}{s^{{a_3}}_{{" + L"a_2}}}{s^{{a_4}}_{{a_1}}}} - " + L"{{{2}}{s^{{i_4}}_{{i_1}}}{s^{{i_3}}_{{i_2}}}{s^{{a_4}}_{{a_2}}}{" + L"s^{{a_3}}_{{a_1}}}} - " + L"{{{2}}{s^{{i_3}}_{{i_1}}}{s^{{i_4}}_{{i_2}}}{s^{{a_3}}_{{a_2}}}{" + L"s^{{a_4}}_{{a_1}}}} + " + L"{{{4}}{s^{{i_3}}_{{i_1}}}{s^{{i_4}}_{{i_2}}}{s^{{a_4}}_{{a_2}}}{" + L"s^{{a_3}}_{{a_1}}}}\\bigr) }"); + } + // two (pure qp) 3-body operators + { + auto opseq = FNOperatorSeq( + {FNOperator({L"i_1", L"i_2", L"i_3"}, {L"a_1", L"a_2", L"a_3"}), + FNOperator({L"a_4", L"a_5", L"a_6"}, {L"i_4", L"i_5", L"i_6"})}); + auto wick = FWickTheorem{opseq}; + REQUIRE_NOTHROW(wick.compute()); + auto result = wick.compute(); + REQUIRE(result->is()); + REQUIRE(result->size() == 36); + } - // two general 2-body operators - { - auto opseq = - FNOperatorSeq({FNOperator({L"p_1", L"p_2"}, {L"p_3", L"p_4"}), - FNOperator({L"p_5", L"p_6"}, {L"p_7", L"p_8"})}); - auto wick = FWickTheorem{opseq}; - REQUIRE_NOTHROW(wick.compute()); - auto result = wick.compute(); - REQUIRE(result->is()); - REQUIRE(result->size() == 4); - } - // two general 2-body operators, partial contractions: Eqs. 22 of - // DOI 10.1063/1.474405 - { - auto opseq = - FNOperatorSeq({FNOperator({L"p_1", L"p_2"}, {L"p_3", L"p_4"}), - FNOperator({L"p_5", L"p_6"}, {L"p_7", L"p_8"})}); - auto wick = FWickTheorem{opseq}; - REQUIRE_NOTHROW(wick.full_contractions(false).compute()); - auto result = wick.full_contractions(false).compute(); - REQUIRE(result->is()); - REQUIRE(result->size() == 49); // the MK paper only gives 47 terms, - // misses the 2 double-hole contractions - } - // one general 2-body operator and one 2-body excitation operator - { - auto opseq = - FNOperatorSeq({FNOperator({L"p_1", L"p_2"}, {L"p_3", L"p_4"}), - FNOperator({L"a_3", L"a_4"}, {L"i_3", L"i_4"})}); - auto wick = FWickTheorem{opseq}; - REQUIRE_NOTHROW(wick.compute()); - auto result = wick.compute(); - std::wcout << "<" << to_latex(opseq) << "> = " << to_latex(result) - << std::endl; - REQUIRE(result->is()); - REQUIRE(result->size() == 4); - REQUIRE(to_latex(result) == - L"{ " - L"\\bigl({{s^{{p_1}}_{{i_4}}}{s^{{p_2}}_{{i_3}}}{s^{{a_3}}_{{p_4}}}" - L"{s^{{a_4}}_{{p_3}}}} - " - L"{{s^{{p_1}}_{{i_4}}}{s^{{p_2}}_{{i_3}}}{s^{{a_4}}_{{p_4}}}{s^{{a_" - L"3}}_{{p_3}}}} - " - L"{{s^{{p_1}}_{{i_3}}}{s^{{p_2}}_{{i_4}}}{s^{{a_3}}_{{p_4}}}{s^{{a_" - L"4}}_{{p_3}}}} + " - L"{{s^{{p_1}}_{{i_3}}}{s^{{p_2}}_{{i_4}}}{s^{{a_4}}_{{p_4}}}{s^{{a_" - L"3}}_{{p_3}}}}\\bigr) }"); - - // in spin-free result first and fourth terms are multiplied by 4, second - // and third terms multiplied by 2 - auto raii_tmp = switch_to_spinfree_context(); - auto result_sf = wick.compute(); - auto result_sf_latex = to_latex(result_sf); - // std::wcout << "<" << to_latex(opseq) << "> = " << to_latex(result_sf) - // << std::endl; - REQUIRE(result_sf->is()); - REQUIRE(result_sf->size() == 4); - REQUIRE(to_latex(result_sf) == - L"{ " - L"\\bigl({{{4}}{s^{{p_1}}_{{i_4}}}{s^{{p_2}}_{{i_3}}}{s^{{a_3}}_{{" - L"p_4}}}{s^{{a_4}}_{{p_3}}}} - " - L"{{{2}}{s^{{p_1}}_{{i_4}}}{s^{{p_2}}_{{i_3}}}{s^{{a_4}}_{{p_4}}}{" - L"s^{{a_3}}_{{p_3}}}} - " - L"{{{2}}{s^{{p_1}}_{{i_3}}}{s^{{p_2}}_{{i_4}}}{s^{{a_3}}_{{p_4}}}{" - L"s^{{a_4}}_{{p_3}}}} + " - L"{{{4}}{s^{{p_1}}_{{i_3}}}{s^{{p_2}}_{{i_4}}}{s^{{a_4}}_{{p_4}}}{" - L"s^{{a_3}}_{{p_3}}}}\\bigr) }"); - } + // one general 1-body operator + one general 2-body operator, partial + // contraction: Eq. 9 of DOI 10.1063/1.474405 + { + auto opseq = + FNOperatorSeq({FNOperator({L"p_1"}, {L"p_2"}), + FNOperator({L"p_3", L"p_4"}, {L"p_5", L"p_6"})}); + auto wick = FWickTheorem{opseq}; + REQUIRE_NOTHROW(wick.full_contractions(false).compute()); + auto result = wick.full_contractions(false).compute(); + REQUIRE(result->is()); + REQUIRE(result->size() == 9); + } - // two general 3-body operators - { - auto opseq = FNOperatorSeq( - {FNOperator({L"p_1", L"p_2", L"p_3"}, {L"p_4", L"p_5", L"p_6"}), - FNOperator({L"p_7", L"p_8", L"p_9"}, {L"p_10", L"p_11", L"p_12"})}); - auto wick = FWickTheorem{opseq}; - REQUIRE_NOTHROW(wick.compute()); - auto result = wick.compute(); - REQUIRE(result->is()); - REQUIRE(result->size() == 36); - } + // two general 2-body operators + { + auto opseq = + FNOperatorSeq({FNOperator({L"p_1", L"p_2"}, {L"p_3", L"p_4"}), + FNOperator({L"p_5", L"p_6"}, {L"p_7", L"p_8"})}); + auto wick = FWickTheorem{opseq}; + REQUIRE_NOTHROW(wick.compute()); + auto result = wick.compute(); + REQUIRE(result->is()); + REQUIRE(result->size() == 4); + } + // two general 2-body operators, partial contractions: Eqs. 22 of + // DOI 10.1063/1.474405 + { + auto opseq = + FNOperatorSeq({FNOperator({L"p_1", L"p_2"}, {L"p_3", L"p_4"}), + FNOperator({L"p_5", L"p_6"}, {L"p_7", L"p_8"})}); + auto wick = FWickTheorem{opseq}; + REQUIRE_NOTHROW(wick.full_contractions(false).compute()); + auto result = wick.full_contractions(false).compute(); + REQUIRE(result->is()); + REQUIRE(result->size() == 49); // the MK paper only gives 47 terms, + // misses the 2 double-hole contractions + } + // one general 2-body operator and one 2-body excitation operator + { + auto opseq = + FNOperatorSeq({FNOperator({L"p_1", L"p_2"}, {L"p_3", L"p_4"}), + FNOperator({L"a_3", L"a_4"}, {L"i_3", L"i_4"})}); + auto wick = FWickTheorem{opseq}; + REQUIRE_NOTHROW(wick.compute()); + auto result = wick.compute(); + std::wcout << "<" << to_latex(opseq) << "> = " << to_latex(result) + << std::endl; + REQUIRE(result->is()); + REQUIRE(result->size() == 4); + REQUIRE( + to_latex(result) == + L"{ " + L"\\bigl({{s^{{p_1}}_{{i_4}}}{s^{{p_2}}_{{i_3}}}{s^{{a_3}}_{{p_4}}}" + L"{s^{{a_4}}_{{p_3}}}} - " + L"{{s^{{p_1}}_{{i_4}}}{s^{{p_2}}_{{i_3}}}{s^{{a_4}}_{{p_4}}}{s^{{a_" + L"3}}_{{p_3}}}} - " + L"{{s^{{p_1}}_{{i_3}}}{s^{{p_2}}_{{i_4}}}{s^{{a_3}}_{{p_4}}}{s^{{a_" + L"4}}_{{p_3}}}} + " + L"{{s^{{p_1}}_{{i_3}}}{s^{{p_2}}_{{i_4}}}{s^{{a_4}}_{{p_4}}}{s^{{a_" + L"3}}_{{p_3}}}}\\bigr) }"); + + // in spin-free result first and fourth terms are multiplied by 4, second + // and third terms multiplied by 2 + auto raii_tmp = switch_to_spinfree_context(); + auto result_sf = wick.compute(); + auto result_sf_latex = to_latex(result_sf); + // std::wcout << "<" << to_latex(opseq) << "> = " << + // to_latex(result_sf) + // << std::endl; + REQUIRE(result_sf->is()); + REQUIRE(result_sf->size() == 4); + REQUIRE( + to_latex(result_sf) == + L"{ " + L"\\bigl({{{4}}{s^{{p_1}}_{{i_4}}}{s^{{p_2}}_{{i_3}}}{s^{{a_3}}_{{" + L"p_4}}}{s^{{a_4}}_{{p_3}}}} - " + L"{{{2}}{s^{{p_1}}_{{i_4}}}{s^{{p_2}}_{{i_3}}}{s^{{a_4}}_{{p_4}}}{" + L"s^{{a_3}}_{{p_3}}}} - " + L"{{{2}}{s^{{p_1}}_{{i_3}}}{s^{{p_2}}_{{i_4}}}{s^{{a_3}}_{{p_4}}}{" + L"s^{{a_4}}_{{p_3}}}} + " + L"{{{4}}{s^{{p_1}}_{{i_3}}}{s^{{p_2}}_{{i_4}}}{s^{{a_4}}_{{p_4}}}{" + L"s^{{a_3}}_{{p_3}}}}\\bigr) }"); + } - // two N-nonconserving operators - { - auto opseq = FNOperatorSeq( - {FNOperator({L"p_1", L"p_2", L"p_3"}, {L"p_4", L"p_5"}), - FNOperator({L"p_7", L"p_8"}, {L"p_10", L"p_11", L"p_12"})}); - auto wick = FWickTheorem{opseq}; - REQUIRE_NOTHROW(wick.compute()); - auto result = wick.compute(); - REQUIRE(result->is()); - REQUIRE(result->size() == 12); - } + // two general 3-body operators + { + auto opseq = FNOperatorSeq( + {FNOperator({L"p_1", L"p_2", L"p_3"}, {L"p_4", L"p_5", L"p_6"}), + FNOperator({L"p_7", L"p_8", L"p_9"}, {L"p_10", L"p_11", L"p_12"})}); + auto wick = FWickTheorem{opseq}; + REQUIRE_NOTHROW(wick.compute()); + auto result = wick.compute(); + REQUIRE(result->is()); + REQUIRE(result->size() == 36); + } + + // two N-nonconserving operators + { + auto opseq = FNOperatorSeq( + {FNOperator({L"p_1", L"p_2", L"p_3"}, {L"p_4", L"p_5"}), + FNOperator({L"p_7", L"p_8"}, {L"p_10", L"p_11", L"p_12"})}); + auto wick = FWickTheorem{opseq}; + REQUIRE_NOTHROW(wick.compute()); + auto result = wick.compute(); + REQUIRE(result->is()); + REQUIRE(result->size() == 12); + } // more N-nonconserving operators { @@ -609,159 +626,159 @@ SECTION("fermi vacuum") { L"{{{-1}}{\\bar{g}^{{i_1}{a_2}}_{{a_3}{a_4}}}}"); } - // odd number of ops -> full contraction is 0 - { - auto opseq = FNOperatorSeq( - {FNOperator({L"p_1", L"p_2"}, {L"p_4", L"p_5"}), - FNOperator({L"p_7", L"p_8"}, {L"p_10", L"p_11", L"p_12"})}); - auto wick = FWickTheorem{opseq}; - REQUIRE_NOTHROW(wick.compute()); - auto result = wick.compute(); - REQUIRE(result->is()); - REQUIRE(result->as().value() == 0); - } - - // 4-body ^ 4-body - SEQUANT_PROFILE_SINGLE( - "wick(4^4)", - { - auto opseq = - FNOperatorSeq({FNOperator({L"p_1", L"p_2", L"p_3", L"p_4"}, - {L"p_5", L"p_6", L"p_7", L"p_8"}), - FNOperator({L"p_21", L"p_22", L"p_23", L"p_24"}, - {L"p_25", L"p_26", L"p_27", L"p_28"})}); - auto wick = FWickTheorem{opseq}; - auto result = wick.compute(true); - REQUIRE(result->is()); - REQUIRE(result->as().value() == 576); - }) + // odd number of ops -> full contraction is 0 + { + auto opseq = FNOperatorSeq( + {FNOperator({L"p_1", L"p_2"}, {L"p_4", L"p_5"}), + FNOperator({L"p_7", L"p_8"}, {L"p_10", L"p_11", L"p_12"})}); + auto wick = FWickTheorem{opseq}; + REQUIRE_NOTHROW(wick.compute()); + auto result = wick.compute(); + REQUIRE(result->is()); + REQUIRE(result->as().value() == 0); + } - // three general 1-body operators - { - auto opseq = FNOperatorSeq({FNOperator({L"p_1"}, {L"p_2"}), - FNOperator({L"p_3"}, {L"p_4"}), - FNOperator({L"p_5"}, {L"p_6"})}); - auto wick = FWickTheorem{opseq}; - REQUIRE_NOTHROW(wick.compute()); - auto result = wick.compute(); - REQUIRE(result->is()); - REQUIRE(result->size() == 2); - } + // 4-body ^ 4-body + SEQUANT_PROFILE_SINGLE( + "wick(4^4)", + { + auto opseq = + FNOperatorSeq({FNOperator({L"p_1", L"p_2", L"p_3", L"p_4"}, + {L"p_5", L"p_6", L"p_7", L"p_8"}), + FNOperator({L"p_21", L"p_22", L"p_23", L"p_24"}, + {L"p_25", L"p_26", L"p_27", L"p_28"})}); + auto wick = FWickTheorem{opseq}; + auto result = wick.compute(true); + REQUIRE(result->is()); + REQUIRE(result->as().value() == 576); + }) + + // three general 1-body operators + { + auto opseq = FNOperatorSeq({FNOperator({L"p_1"}, {L"p_2"}), + FNOperator({L"p_3"}, {L"p_4"}), + FNOperator({L"p_5"}, {L"p_6"})}); + auto wick = FWickTheorem{opseq}; + REQUIRE_NOTHROW(wick.compute()); + auto result = wick.compute(); + REQUIRE(result->is()); + REQUIRE(result->size() == 2); + } - // 4 general 1-body operators - { - auto opseq = FNOperatorSeq( - {FNOperator({L"p_1"}, {L"p_2"}), FNOperator({L"p_3"}, {L"p_4"}), - FNOperator({L"p_5"}, {L"p_6"}), FNOperator({L"p_7"}, {L"p_8"})}); - auto ext_indices = make_indices>(WstrList{ - L"p_1", L"p_2", L"p_3", L"p_4", L"p_5", L"p_6", L"p_7", L"p_8"}); - auto wick1 = FWickTheorem{opseq}; - auto result1 = wick1.set_external_indices(ext_indices).compute(); - REQUIRE(result1->is()); - REQUIRE(result1->size() == 9); - auto wick2 = FWickTheorem{opseq}; - auto result2 = wick2.set_external_indices(ext_indices) - .set_nop_connections({{1, 2}, {1, 3}}) - .compute(); - REQUIRE(result2->is()); - REQUIRE(result2->size() == 2); - } + // 4 general 1-body operators + { + auto opseq = FNOperatorSeq( + {FNOperator({L"p_1"}, {L"p_2"}), FNOperator({L"p_3"}, {L"p_4"}), + FNOperator({L"p_5"}, {L"p_6"}), FNOperator({L"p_7"}, {L"p_8"})}); + auto ext_indices = make_indices>(WstrList{ + L"p_1", L"p_2", L"p_3", L"p_4", L"p_5", L"p_6", L"p_7", L"p_8"}); + auto wick1 = FWickTheorem{opseq}; + auto result1 = wick1.set_external_indices(ext_indices).compute(); + REQUIRE(result1->is()); + REQUIRE(result1->size() == 9); + auto wick2 = FWickTheorem{opseq}; + auto result2 = wick2.set_external_indices(ext_indices) + .set_nop_connections({{1, 2}, {1, 3}}) + .compute(); + REQUIRE(result2->is()); + REQUIRE(result2->size() == 2); + } - // 4-body ^ 2-body ^ 2-body - { - auto opseq = - FNOperatorSeq({FNOperator({L"p_1", L"p_2", L"p_3", L"p_4"}, - {L"p_5", L"p_6", L"p_7", L"p_8"}), - FNOperator({L"p_9", L"p_10"}, {L"p_11", L"p_12"}), - FNOperator({L"p_13", L"p_14"}, {L"p_15", L"p_16"})}); - auto wick = FWickTheorem{opseq}; - auto result = wick.compute(); - REQUIRE(result->is()); - REQUIRE(result->size() == 576); - } + // 4-body ^ 2-body ^ 2-body + { + auto opseq = + FNOperatorSeq({FNOperator({L"p_1", L"p_2", L"p_3", L"p_4"}, + {L"p_5", L"p_6", L"p_7", L"p_8"}), + FNOperator({L"p_9", L"p_10"}, {L"p_11", L"p_12"}), + FNOperator({L"p_13", L"p_14"}, {L"p_15", L"p_16"})}); + auto wick = FWickTheorem{opseq}; + auto result = wick.compute(); + REQUIRE(result->is()); + REQUIRE(result->size() == 576); + } - // 2-body ^ 2-body ^ 2-body - { - auto opseq = - FNOperatorSeq({FNOperator({L"p_1", L"p_2"}, {L"p_5", L"p_6"}), - FNOperator({L"p_9", L"p_10"}, {L"p_11", L"p_12"}), - FNOperator({L"p_17", L"p_18"}, {L"p_19", L"p_20"})}); - auto wick = FWickTheorem{opseq}; - auto result = wick.compute(); - REQUIRE(result->is()); - REQUIRE(result->size() == 80); - } + // 2-body ^ 2-body ^ 2-body + { + auto opseq = + FNOperatorSeq({FNOperator({L"p_1", L"p_2"}, {L"p_5", L"p_6"}), + FNOperator({L"p_9", L"p_10"}, {L"p_11", L"p_12"}), + FNOperator({L"p_17", L"p_18"}, {L"p_19", L"p_20"})}); + auto wick = FWickTheorem{opseq}; + auto result = wick.compute(); + REQUIRE(result->is()); + REQUIRE(result->size() == 80); + } - // 2-body ^ 2-body ^ 2-body ^ 2-body - SEQUANT_PROFILE_SINGLE("wick(2^2^2^2)", { - auto opseq = - FNOperatorSeq({FNOperator({L"p_1", L"p_2"}, {L"p_5", L"p_6"}), - FNOperator({L"p_9", L"p_10"}, {L"p_11", L"p_12"}), - FNOperator({L"p_13", L"p_14"}, {L"p_15", L"p_16"}), - FNOperator({L"p_17", L"p_18"}, {L"p_19", L"p_20"})}); - auto wick = FWickTheorem{opseq}; - auto result = wick.compute(true); - REQUIRE(result->is()); - REQUIRE(result->as().value() == 4752); - }) + // 2-body ^ 2-body ^ 2-body ^ 2-body + SEQUANT_PROFILE_SINGLE("wick(2^2^2^2)", { + auto opseq = + FNOperatorSeq({FNOperator({L"p_1", L"p_2"}, {L"p_5", L"p_6"}), + FNOperator({L"p_9", L"p_10"}, {L"p_11", L"p_12"}), + FNOperator({L"p_13", L"p_14"}, {L"p_15", L"p_16"}), + FNOperator({L"p_17", L"p_18"}, {L"p_19", L"p_20"})}); + auto wick = FWickTheorem{opseq}; + auto result = wick.compute(true); + REQUIRE(result->is()); + REQUIRE(result->as().value() == 4752); + }) #ifndef SEQUANT_SKIP_LONG_TESTS - // 4-body ^ 2-body ^ 2-body ^ 2-body - SEQUANT_PROFILE_SINGLE("wick(4^2^2^2)", { - auto opseq = - FNOperatorSeq({FNOperator({L"p_1", L"p_2", L"p_3", L"p_4"}, - {L"p_5", L"p_6", L"p_7", L"p_8"}), - FNOperator({L"p_9", L"p_10"}, {L"p_11", L"p_12"}), - FNOperator({L"p_13", L"p_14"}, {L"p_15", L"p_16"}), - FNOperator({L"p_17", L"p_18"}, {L"p_19", L"p_20"})}); - auto wick = FWickTheorem{opseq}; - auto result = wick.use_topology(true).compute(true); - REQUIRE(result->is()); - REQUIRE(result->as().value() == 2088); - }) - - // 3-body ^ 2-body ^ 2-body ^ 3-body - SEQUANT_PROFILE_SINGLE("wick(3^2^2^3)", { - auto opseq = FNOperatorSeq( - {FNOperator({L"p_1", L"p_2", L"p_3"}, {L"p_5", L"p_6", L"p_7"}), - FNOperator({L"p_9", L"p_10"}, {L"p_11", L"p_12"}), - FNOperator({L"p_13", L"p_14"}, {L"p_15", L"p_16"}), - FNOperator({L"p_17", L"p_18", L"p_19"}, {L"p_20", L"p_21", L"p_22"}, - V)}); - auto wick = FWickTheorem{opseq}; - auto result = wick.use_topology(true).compute(true); - REQUIRE(result->is()); - REQUIRE(result->as().value() == 694); - }) - - // 4-body ^ 2-body ^ 4-body - SEQUANT_PROFILE_SINGLE("wick(4^2^4)", { - auto opseq = - FNOperatorSeq({FNOperator({L"p_1", L"p_2", L"p_3", L"p_4"}, - {L"p_5", L"p_6", L"p_7", L"p_8"}), - FNOperator({L"p_9", L"p_10"}, {L"p_11", L"p_12"}), - FNOperator({L"p_21", L"p_22", L"p_23", L"p_24"}, - {L"p_25", L"p_26", L"p_27", L"p_28"})}); - auto wick = FWickTheorem{opseq}; - auto result = wick.use_topology(true).compute(true); - REQUIRE(result->is()); - REQUIRE(result->as().value() == 28); - }) - - // 4-body ^ 4-body ^ 4-body - SEQUANT_PROFILE_SINGLE("wick(4^4^4)", { - auto opseq = - FNOperatorSeq({FNOperator({L"p_1", L"p_2", L"p_3", L"p_4"}, - {L"p_5", L"p_6", L"p_7", L"p_8"}), - FNOperator({L"p_11", L"p_12", L"p_13", L"p_14"}, - {L"p_15", L"p_16", L"p_17", L"p_18"}), - FNOperator({L"p_21", L"p_22", L"p_23", L"p_24"}, - {L"p_25", L"p_26", L"p_27", L"p_28"})}); - auto wick = FWickTheorem{opseq}; - auto result = wick.use_topology(true).compute(true); - REQUIRE(result->is()); - REQUIRE(result->as().value() == 70); - }) + // 4-body ^ 2-body ^ 2-body ^ 2-body + SEQUANT_PROFILE_SINGLE("wick(4^2^2^2)", { + auto opseq = + FNOperatorSeq({FNOperator({L"p_1", L"p_2", L"p_3", L"p_4"}, + {L"p_5", L"p_6", L"p_7", L"p_8"}), + FNOperator({L"p_9", L"p_10"}, {L"p_11", L"p_12"}), + FNOperator({L"p_13", L"p_14"}, {L"p_15", L"p_16"}), + FNOperator({L"p_17", L"p_18"}, {L"p_19", L"p_20"})}); + auto wick = FWickTheorem{opseq}; + auto result = wick.use_topology(true).compute(true); + REQUIRE(result->is()); + REQUIRE(result->as().value() == 2088); + }) + + // 3-body ^ 2-body ^ 2-body ^ 3-body + SEQUANT_PROFILE_SINGLE("wick(3^2^2^3)", { + auto opseq = FNOperatorSeq( + {FNOperator({L"p_1", L"p_2", L"p_3"}, {L"p_5", L"p_6", L"p_7"}), + FNOperator({L"p_9", L"p_10"}, {L"p_11", L"p_12"}), + FNOperator({L"p_13", L"p_14"}, {L"p_15", L"p_16"}), + FNOperator({L"p_17", L"p_18", L"p_19"}, {L"p_20", L"p_21", L"p_22"}, + V)}); + auto wick = FWickTheorem{opseq}; + auto result = wick.use_topology(true).compute(true); + REQUIRE(result->is()); + REQUIRE(result->as().value() == 694); + }) + + // 4-body ^ 2-body ^ 4-body + SEQUANT_PROFILE_SINGLE("wick(4^2^4)", { + auto opseq = + FNOperatorSeq({FNOperator({L"p_1", L"p_2", L"p_3", L"p_4"}, + {L"p_5", L"p_6", L"p_7", L"p_8"}), + FNOperator({L"p_9", L"p_10"}, {L"p_11", L"p_12"}), + FNOperator({L"p_21", L"p_22", L"p_23", L"p_24"}, + {L"p_25", L"p_26", L"p_27", L"p_28"})}); + auto wick = FWickTheorem{opseq}; + auto result = wick.use_topology(true).compute(true); + REQUIRE(result->is()); + REQUIRE(result->as().value() == 28); + }) + + // 4-body ^ 4-body ^ 4-body + SEQUANT_PROFILE_SINGLE("wick(4^4^4)", { + auto opseq = + FNOperatorSeq({FNOperator({L"p_1", L"p_2", L"p_3", L"p_4"}, + {L"p_5", L"p_6", L"p_7", L"p_8"}), + FNOperator({L"p_11", L"p_12", L"p_13", L"p_14"}, + {L"p_15", L"p_16", L"p_17", L"p_18"}), + FNOperator({L"p_21", L"p_22", L"p_23", L"p_24"}, + {L"p_25", L"p_26", L"p_27", L"p_28"})}); + auto wick = FWickTheorem{opseq}; + auto result = wick.use_topology(true).compute(true); + REQUIRE(result->is()); + REQUIRE(result->as().value() == 70); + }) #endif #if 0 @@ -779,72 +796,23 @@ SECTION("fermi vacuum") { auto result = wick.use_topology(true).compute(true); } #endif -} // SECTION("fermi vacuum") - -SECTION("Expression Reduction") { - constexpr Vacuum V = Vacuum::SingleProduct; - // default vacuum is already spin-orbital Fermi vacuum - - auto switch_to_spinfree_context = detail::NoDiscard([&]() { - auto context_sf = get_default_context(); - context_sf.set(SPBasis::spinfree); - return set_scoped_default_context(context_sf); - }); - - // 2-body ^ 2-body - SEQUANT_PROFILE_SINGLE("wick(H2*T2)", { - auto opseq = - FNOperatorSeq({FNOperator({L"p_1", L"p_2"}, {L"p_3", L"p_4"}), - FNOperator({L"a_4", L"a_5"}, {L"i_4", L"i_5"})}); - auto wick = FWickTheorem{opseq}; - auto wick_result = wick.compute(); - REQUIRE(wick_result->is()); - REQUIRE(wick_result->size() == 4); - - // multiply tensor factors and expand - auto wick_result_2 = - ex(L"g", WstrList{L"p_1", L"p_2"}, WstrList{L"p_3", L"p_4"}, - WstrList{}, Symmetry::antisymm) * - ex(L"t", WstrList{L"a_4", L"a_5"}, WstrList{L"i_4", L"i_5"}, - WstrList{}, Symmetry::antisymm) * - wick_result; - expand(wick_result_2); - REQUIRE(to_latex(wick_result_2) == - L"{ " - L"\\bigl({{\\bar{g}^{{p_3}{p_4}}_{{p_1}{p_2}}}{\\bar{t}^{{i_4}{i_" - L"5}}_{{a_4}{a_" - L"5}}}{s^{{p_1}}_{{i_5}}}{s^{{p_2}}_{{i_4}}}{s^{{a_4}}_{{p_4}}}{" - L"s^{{a_5}}_{{p_3}}}} - {" - L"{\\bar{g}^{{p_3}{p_4}}_{{p_1}{p_2}}}{\\bar{t}^{{i_4}{i_5}}_{{a_" - L"4}{a_5}}}{s^{{" - L"p_1}}_{{i_5}}}{s^{{p_2}}_{{i_4}}}{s^{{a_5}}_{{p_4}}}{s^{{a_4}}_" - L"{{p_3}}}} - {" - L"{\\bar{g}^{{p_3}{p_4}}_{{p_1}{p_2}}}{\\bar{t}^{{i_4}{i_5}}_{{a_" - L"4}{a_5}}}{s^{{" - L"p_1}}_{{i_4}}}{s^{{p_2}}_{{i_5}}}{s^{{a_4}}_{{p_4}}}{s^{{a_5}}_" - L"{{p_3}}}} + " - L"{{\\bar{g}^{{p_3}{p_4}}_{{p_1}{p_2}}}{\\bar{t}^{{i_4}{i_5}}_{{" - L"a_4}{a_5}}}{s^{" - L"{p_1}}_{{i_4}}}{s^{{p_2}}_{{i_5}}}{s^{{a_5}}_{{p_4}}}{s^{{a_4}}" - L"_{{p_3}}}}\\bigr) }"); - wick.reduce(wick_result_2); - rapid_simplify(wick_result_2); - TensorCanonicalizer::register_instance( - std::make_shared()); - canonicalize(wick_result_2); - rapid_simplify(wick_result_2); - - std::wcout << L"H2*T2 = " << to_latex(wick_result_2) << std::endl; - std::wcout << L"H2*T2 = " << to_wolfram(wick_result_2) << std::endl; - REQUIRE(to_latex(wick_result_2) == - L"{{{4}}" - L"{\\bar{g}^{{a_1}{a_2}}_{{i_1}{i_2}}}{\\bar{t}^{{i_1}{i_2}}_{{a_" - L"1}{a_2}}}}"); - - // spin-free case will produce 2 terms - { - auto raii_tmp = switch_to_spinfree_context(); + } // SECTION("fermi vacuum") + + SECTION("Expression Reduction") { + constexpr Vacuum V = Vacuum::SingleProduct; + // default vacuum is already spin-orbital Fermi vacuum + + auto switch_to_spinfree_context = detail::NoDiscard([&]() { + auto context_sf = get_default_context(); + context_sf.set(SPBasis::spinfree); + return set_scoped_default_context(context_sf); + }); + // 2-body ^ 2-body + SEQUANT_PROFILE_SINGLE("wick(H2*T2)", { + auto opseq = + FNOperatorSeq({FNOperator({L"p_1", L"p_2"}, {L"p_3", L"p_4"}), + FNOperator({L"a_4", L"a_5"}, {L"i_4", L"i_5"})}); auto wick = FWickTheorem{opseq}; auto wick_result = wick.compute(); REQUIRE(wick_result->is()); @@ -853,258 +821,306 @@ SECTION("Expression Reduction") { // multiply tensor factors and expand auto wick_result_2 = ex(L"g", WstrList{L"p_1", L"p_2"}, WstrList{L"p_3", L"p_4"}, - WstrList{}, Symmetry::nonsymm) * + WstrList{}, Symmetry::antisymm) * ex(L"t", WstrList{L"a_4", L"a_5"}, WstrList{L"i_4", L"i_5"}, - WstrList{}, Symmetry::nonsymm) * + WstrList{}, Symmetry::antisymm) * wick_result; expand(wick_result_2); - REQUIRE(wick_result_2->size() == 4); // still 4 terms - + REQUIRE(to_latex(wick_result_2) == + L"{ " + L"\\bigl({{\\bar{g}^{{p_3}{p_4}}_{{p_1}{p_2}}}{\\bar{t}^{{i_4}{i_" + L"5}}_{{a_4}{a_" + L"5}}}{s^{{p_1}}_{{i_5}}}{s^{{p_2}}_{{i_4}}}{s^{{a_4}}_{{p_4}}}{" + L"s^{{a_5}}_{{p_3}}}} - {" + L"{\\bar{g}^{{p_3}{p_4}}_{{p_1}{p_2}}}{\\bar{t}^{{i_4}{i_5}}_{{a_" + L"4}{a_5}}}{s^{{" + L"p_1}}_{{i_5}}}{s^{{p_2}}_{{i_4}}}{s^{{a_5}}_{{p_4}}}{s^{{a_4}}_" + L"{{p_3}}}} - {" + L"{\\bar{g}^{{p_3}{p_4}}_{{p_1}{p_2}}}{\\bar{t}^{{i_4}{i_5}}_{{a_" + L"4}{a_5}}}{s^{{" + L"p_1}}_{{i_4}}}{s^{{p_2}}_{{i_5}}}{s^{{a_4}}_{{p_4}}}{s^{{a_5}}_" + L"{{p_3}}}} + " + L"{{\\bar{g}^{{p_3}{p_4}}_{{p_1}{p_2}}}{\\bar{t}^{{i_4}{i_5}}_{{" + L"a_4}{a_5}}}{s^{" + L"{p_1}}_{{i_4}}}{s^{{p_2}}_{{i_5}}}{s^{{a_5}}_{{p_4}}}{s^{{a_4}}" + L"_{{p_3}}}}\\bigr) }"); wick.reduce(wick_result_2); rapid_simplify(wick_result_2); TensorCanonicalizer::register_instance( std::make_shared()); canonicalize(wick_result_2); rapid_simplify(wick_result_2); - REQUIRE(wick_result_2->size() == 2); // now 2 terms - std::wcout << L"spinfree H2*T2 = " << to_latex(wick_result_2) - << std::endl; - REQUIRE_SUM_EQUAL( - wick_result_2, - L"{ \\bigl( - " - L"{{{4}}{g^{{a_1}{a_2}}_{{i_1}{i_2}}}{t^{{i_2}{i_1}}_{{a_1}{a_2}}" - L"}} + " - L"{{{8}}{g^{{a_1}{a_2}}_{{i_1}{i_2}}}{t^{{i_1}{i_2}}_{{a_1}{a_2}}" - L"}}\\bigr) }"); - } - }); - - // 2-body ^ 1-body ^ 1-body, with/without using topology - SEQUANT_PROFILE_SINGLE("wick(H2*T1*T1)", { - for (auto&& use_nop_partitions : {false}) { - for (auto&& use_op_partitions : {true, false}) { - std::wostringstream oss; - oss << "use_op_partitions=" << use_op_partitions << "}: H2*T1*T1 = "; - - auto opseq = FNOperatorSeq( - {FNOperator({L"p_1", L"p_2"}, {L"p_3", L"p_4"}), - FNOperator({L"a_4"}, {L"i_4"}), FNOperator({L"a_5"}, {L"i_5"})}); - auto wick = FWickTheorem{opseq}; - wick.use_topology(use_nop_partitions || use_op_partitions); - // if (use_nop_partitions) wick.set_nop_partitions({{1, 2}}); - if (use_op_partitions) wick.set_op_partitions({{0, 1}, {2, 3}}); - auto wick_result = wick.compute(); - // print(oss.str() + L" (nopseq only) ", wick_result); - if (use_op_partitions) { - REQUIRE(wick_result->is()); - REQUIRE(wick_result->size() == 4 /* factors */); - } else { - REQUIRE(wick_result->is()); - REQUIRE(wick_result->size() == 4 /* summands */); - } + std::wcout << L"H2*T2 = " << to_latex(wick_result_2) << std::endl; + std::wcout << L"H2*T2 = " << to_wolfram(wick_result_2) << std::endl; + REQUIRE(to_latex(wick_result_2) == + L"{{{4}}" + L"{\\bar{g}^{{a_1}{a_2}}_{{i_1}{i_2}}}{\\bar{t}^{{i_1}{i_2}}_{{a_" + L"1}{a_2}}}}"); - // multiply tensor factors and expand - auto wick_result_2 = - ex(L"g", WstrList{L"p_1", L"p_2"}, WstrList{L"p_3", L"p_4"}, - WstrList{}, Symmetry::antisymm) * - ex(L"t", WstrList{L"a_4"}, WstrList{L"i_4"}, WstrList{}, - Symmetry::antisymm) * - ex(L"t", WstrList{L"a_5"}, WstrList{L"i_5"}, WstrList{}, - Symmetry::antisymm) * - wick_result; - expand(wick_result_2); - wick.reduce(wick_result_2); - rapid_simplify(wick_result_2); - TensorCanonicalizer::register_instance( - std::make_shared(std::vector{})); - canonicalize(wick_result_2); - rapid_simplify(wick_result_2); + // spin-free case will produce 2 terms + { + auto raii_tmp = switch_to_spinfree_context(); - // print(oss.str(), wick_result_2); - REQUIRE(wick_result_2->size() == 3 /* factors */); - REQUIRE(to_latex(wick_result_2) == - L"{{{4}}" - L"{\\bar{g}^{{a_1}{a_2}}_{{i_1}{i_2}}}{t^{{i_1}}_{{a_1}}}{" - L"t^{{i_" - L"2}}_{{a_" - L"2}}}}"); - } // use_op_partitions - } // use_nop_partitions - }); - - // 2=body ^ 1-body ^ 2-body with dependent (PNO) indices - SEQUANT_PROFILE_SINGLE("wick(P2*H1*T2)", { - auto opseq = FNOperatorSeq({FNOperator(IndexList{L"i_1", L"i_2"}, - {Index(L"a_1", {L"i_1", L"i_2"}), - Index(L"a_2", {L"i_1", L"i_2"})}, - V), - FNOperator({L"p_1"}, {L"p_2"}), - FNOperator({Index(L"a_3", {L"i_3", L"i_4"}), - Index(L"a_4", {L"i_3", L"i_4"})}, - IndexList{L"i_3", L"i_4"})}); - auto wick = FWickTheorem{opseq}; - auto wick_result = wick.compute(); - REQUIRE(wick_result->is()); - REQUIRE(wick_result->size() == 16); - - // multiply tensor factors and expand - auto wick_result_2 = - ex( - L"A", IndexList{L"i_1", L"i_2"}, - IndexList{{L"a_1", {L"i_1", L"i_2"}}, {L"a_2", {L"i_1", L"i_2"}}}, - IndexList{}, Symmetry::antisymm) * - ex(L"f", WstrList{L"p_1"}, WstrList{L"p_2"}, WstrList{}, - Symmetry::antisymm) * - ex( - L"t", - IndexList{{L"a_3", {L"i_3", L"i_4"}}, {L"a_4", {L"i_3", L"i_4"}}}, - IndexList{L"i_3", L"i_4"}, IndexList{}, Symmetry::antisymm) * - wick_result; - expand(wick_result_2); - wick.reduce(wick_result_2); - rapid_simplify(wick_result_2); - TensorCanonicalizer::register_instance( - std::make_shared()); - canonicalize(wick_result_2); - rapid_simplify(wick_result_2); - - // std::wcout << L"P2*H1*T2(PNO) = " << to_latex_align(wick_result_2) - // << std::endl; - // it appears that the two terms are swapped when using gcc 8 on linux - // TODO investigate why sum canonicalization seems to produce - // platform-dependent results. - // REQUIRE(to_latex(wick_result_2) == - // L"{ \\bigl( - {{{8}}" - // L"{A^{{a_1^{{i_1}{i_2}}}{a_2^{{i_1}{i_2}}}}_{{i_1}{i_2}}}{f^{{a_" - // L"3^{{i_1}{i_2}}}}_{{a_1^{{i_1}{i_2}}}}}{t^{{i_1}{i_2}}_{{a_2^{{" - // L"i_1}{i_2}}}{a_3^{{i_1}{i_2}}}}}} + {{{8}}" - // L"{A^{{a_1^{{i_1}{i_2}}}{a_2^{{i_1}{i_2}}}}_{{i_1}{i_2}}}{f^{{i_" - // L"1}}_{{i_3}}}{t^{{i_2}{i_3}}_{{a_3^{{i_2}{i_3}}}{a_4^{{i_2}{i_3}" - // L"}}}}{s^{{a_3^{{i_2}{i_3}}}}_{{a_1^{{i_1}{i_2}}}}}{s^{{a_4^{{i_" - // L"2}{i_3}}}}_{{a_2^{{i_1}{i_2}}}}}}\\bigr) }"); - }); - - // 2=body ^ 2-body ^ 2-body ^ 2-body with dependent (PNO) indices - SEQUANT_PROFILE_SINGLE("wick(P2*H2*T2*T2)", { - for (auto&& use_nop_partitions : {false}) { - for (auto&& use_op_partitions : {true, false}) { - std::wostringstream oss; - oss << "use_{nop,op}_partitions={" << use_nop_partitions << "," - << use_op_partitions << "}: P2*H2*T2*T2(PNO) = "; - - auto opseq = - FNOperatorSeq({FNOperator(IndexList{L"i_1", L"i_2"}, - {Index(L"a_1", {L"i_1", L"i_2"}), - Index(L"a_2", {L"i_1", L"i_2"})}, - V), - FNOperator({L"p_1", L"p_2"}, {L"p_3", L"p_4"}), - FNOperator({Index(L"a_3", {L"i_3", L"i_4"}), - Index(L"a_4", {L"i_3", L"i_4"})}, - IndexList{L"i_3", L"i_4"}), - FNOperator({Index(L"a_5", {L"i_5", L"i_6"}), - Index(L"a_6", {L"i_5", L"i_6"})}, - IndexList{L"i_5", L"i_6"})}); auto wick = FWickTheorem{opseq}; - wick.set_nop_connections({{1, 2}, {1, 3}}).use_topology(true); - - if (use_nop_partitions) wick.set_nop_partitions({{2, 3}}); - if (use_op_partitions) - wick.set_op_partitions({{0, 1}, - {2, 3}, - {4, 5}, - {6, 7}, - {8, 9}, - {10, 11}, - {12, 13}, - {14, 15}}); auto wick_result = wick.compute(); REQUIRE(wick_result->is()); - if (use_op_partitions) { - REQUIRE(wick_result->size() == 7); - } else { - REQUIRE(wick_result->size() == 544); - } + REQUIRE(wick_result->size() == 4); // multiply tensor factors and expand auto wick_result_2 = - ex(rational{1, 256}) * - ex(L"A", IndexList{L"i_1", L"i_2"}, - IndexList{{L"a_1", {L"i_1", L"i_2"}}, - {L"a_2", {L"i_1", L"i_2"}}}, - IndexList{}, Symmetry::antisymm) * ex(L"g", WstrList{L"p_1", L"p_2"}, WstrList{L"p_3", L"p_4"}, - WstrList{}, Symmetry::antisymm) * - ex(L"t", - IndexList{{L"a_3", {L"i_3", L"i_4"}}, - {L"a_4", {L"i_3", L"i_4"}}}, - IndexList{L"i_3", L"i_4"}, IndexList{}, - Symmetry::antisymm) * - ex(L"t", - IndexList{{L"a_5", {L"i_5", L"i_6"}}, - {L"a_6", {L"i_5", L"i_6"}}}, - IndexList{L"i_5", L"i_6"}, IndexList{}, - Symmetry::antisymm) * + WstrList{}, Symmetry::nonsymm) * + ex(L"t", WstrList{L"a_4", L"a_5"}, WstrList{L"i_4", L"i_5"}, + WstrList{}, Symmetry::nonsymm) * wick_result; expand(wick_result_2); + REQUIRE(wick_result_2->size() == 4); // still 4 terms + wick.reduce(wick_result_2); rapid_simplify(wick_result_2); TensorCanonicalizer::register_instance( std::make_shared()); canonicalize(wick_result_2); - canonicalize(wick_result_2); - canonicalize(wick_result_2); rapid_simplify(wick_result_2); + REQUIRE(wick_result_2->size() == 2); // now 2 terms - // std::wcout << oss.str() << to_latex_align(wick_result_2, 20) - // << std::endl; - REQUIRE(wick_result_2->is()); - REQUIRE(wick_result_2->size() == 4); - } // use_op_partitions - } // use_nop_partitions - }); + std::wcout << L"spinfree H2*T2 = " << to_latex(wick_result_2) + << std::endl; + REQUIRE_SUM_EQUAL( + wick_result_2, + L"{ \\bigl( - " + L"{{{4}}{g^{{a_1}{a_2}}_{{i_1}{i_2}}}{t^{{i_2}{i_1}}_{{a_1}{a_2}}" + L"}} + " + L"{{{8}}{g^{{a_1}{a_2}}_{{i_1}{i_2}}}{t^{{i_1}{i_2}}_{{a_1}{a_2}}" + L"}}\\bigr) }"); + } + }); + + // 2-body ^ 1-body ^ 1-body, with/without using topology + SEQUANT_PROFILE_SINGLE("wick(H2*T1*T1)", { + for (auto&& use_nop_partitions : {false}) { + for (auto&& use_op_partitions : {true, false}) { + std::wostringstream oss; + oss << "use_op_partitions=" << use_op_partitions << "}: H2*T1*T1 = "; + + auto opseq = FNOperatorSeq( + {FNOperator({L"p_1", L"p_2"}, {L"p_3", L"p_4"}), + FNOperator({L"a_4"}, {L"i_4"}), FNOperator({L"a_5"}, {L"i_5"})}); + auto wick = FWickTheorem{opseq}; + wick.use_topology(use_nop_partitions || use_op_partitions); + // if (use_nop_partitions) wick.set_nop_partitions({{1, 2}}); + if (use_op_partitions) wick.set_op_partitions({{0, 1}, {2, 3}}); + auto wick_result = wick.compute(); + // print(oss.str() + L" (nopseq only) ", wick_result); + if (use_op_partitions) { + REQUIRE(wick_result->is()); + REQUIRE(wick_result->size() == 4 /* factors */); + } else { + REQUIRE(wick_result->is()); + REQUIRE(wick_result->size() == 4 /* summands */); + } + + // multiply tensor factors and expand + auto wick_result_2 = + ex(L"g", WstrList{L"p_1", L"p_2"}, + WstrList{L"p_3", L"p_4"}, WstrList{}, Symmetry::antisymm) * + ex(L"t", WstrList{L"a_4"}, WstrList{L"i_4"}, + WstrList{}, Symmetry::antisymm) * + ex(L"t", WstrList{L"a_5"}, WstrList{L"i_5"}, + WstrList{}, Symmetry::antisymm) * + wick_result; + expand(wick_result_2); + wick.reduce(wick_result_2); + rapid_simplify(wick_result_2); + TensorCanonicalizer::register_instance( + std::make_shared( + std::vector{})); + canonicalize(wick_result_2); + rapid_simplify(wick_result_2); + + // print(oss.str(), wick_result_2); + REQUIRE(wick_result_2->size() == 3 /* factors */); + REQUIRE(to_latex(wick_result_2) == + L"{{{4}}" + L"{\\bar{g}^{{a_1}{a_2}}_{{i_1}{i_2}}}{t^{{i_1}}_{{a_1}}}{" + L"t^{{i_" + L"2}}_{{a_" + L"2}}}}"); + } // use_op_partitions + } // use_nop_partitions + }); + + // 2=body ^ 1-body ^ 2-body with dependent (PNO) indices + SEQUANT_PROFILE_SINGLE("wick(P2*H1*T2)", { + auto opseq = FNOperatorSeq({FNOperator(IndexList{L"i_1", L"i_2"}, + {Index(L"a_1", {L"i_1", L"i_2"}), + Index(L"a_2", {L"i_1", L"i_2"})}, + V), + FNOperator({L"p_1"}, {L"p_2"}), + FNOperator({Index(L"a_3", {L"i_3", L"i_4"}), + Index(L"a_4", {L"i_3", L"i_4"})}, + IndexList{L"i_3", L"i_4"})}); + auto wick = FWickTheorem{opseq}; + auto wick_result = wick.compute(); + REQUIRE(wick_result->is()); + REQUIRE(wick_result->size() == 16); + + // multiply tensor factors and expand + auto wick_result_2 = + ex( + L"A", IndexList{L"i_1", L"i_2"}, + IndexList{{L"a_1", {L"i_1", L"i_2"}}, {L"a_2", {L"i_1", L"i_2"}}}, + IndexList{}, Symmetry::antisymm) * + ex(L"f", WstrList{L"p_1"}, WstrList{L"p_2"}, + WstrList{},Symmetry::antisymm) * + ex( + L"t", + IndexList{{L"a_3", {L"i_3", L"i_4"}}, {L"a_4", {L"i_3", L"i_4"}}}, + IndexList{L"i_3", L"i_4"}, IndexList{}, Symmetry::antisymm) * + wick_result; + expand(wick_result_2); + wick.reduce(wick_result_2); + rapid_simplify(wick_result_2); + TensorCanonicalizer::register_instance( + std::make_shared()); + canonicalize(wick_result_2); + rapid_simplify(wick_result_2); + + // std::wcout << L"P2*H1*T2(PNO) = " << to_latex_align(wick_result_2) + // << std::endl; + // it appears that the two terms are swapped when using gcc 8 on linux + // TODO investigate why sum canonicalization seems to produce + // platform-dependent results. + // REQUIRE(to_latex(wick_result_2) == + // L"{ \\bigl( - {{{8}}" + // L"{A^{{a_1^{{i_1}{i_2}}}{a_2^{{i_1}{i_2}}}}_{{i_1}{i_2}}}{f^{{a_" + // L"3^{{i_1}{i_2}}}}_{{a_1^{{i_1}{i_2}}}}}{t^{{i_1}{i_2}}_{{a_2^{{" + // L"i_1}{i_2}}}{a_3^{{i_1}{i_2}}}}}} + {{{8}}" + // L"{A^{{a_1^{{i_1}{i_2}}}{a_2^{{i_1}{i_2}}}}_{{i_1}{i_2}}}{f^{{i_" + // L"1}}_{{i_3}}}{t^{{i_2}{i_3}}_{{a_3^{{i_2}{i_3}}}{a_4^{{i_2}{i_3}" + // L"}}}}{s^{{a_3^{{i_2}{i_3}}}}_{{a_1^{{i_1}{i_2}}}}}{s^{{a_4^{{i_" + // L"2}{i_3}}}}_{{a_2^{{i_1}{i_2}}}}}}\\bigr) }"); + }); + + // 2=body ^ 2-body ^ 2-body ^ 2-body with dependent (PNO) indices + SEQUANT_PROFILE_SINGLE("wick(P2*H2*T2*T2)", { + for (auto&& use_nop_partitions : {false}) { + for (auto&& use_op_partitions : {true, false}) { + std::wostringstream oss; + oss << "use_{nop,op}_partitions={" << use_nop_partitions << "," + << use_op_partitions << "}: P2*H2*T2*T2(PNO) = "; + + auto opseq = + FNOperatorSeq({FNOperator(IndexList{L"i_1", L"i_2"}, + {Index(L"a_1", {L"i_1", L"i_2"}), + Index(L"a_2", {L"i_1", L"i_2"})}, + V), + FNOperator({L"p_1", L"p_2"}, {L"p_3", L"p_4"}), + FNOperator({Index(L"a_3", {L"i_3", L"i_4"}), + Index(L"a_4", {L"i_3", L"i_4"})}, + IndexList{L"i_3", L"i_4"}), + FNOperator({Index(L"a_5", {L"i_5", L"i_6"}), + Index(L"a_6", {L"i_5", L"i_6"})}, + IndexList{L"i_5", L"i_6"})}); + auto wick = FWickTheorem{opseq}; + wick.set_nop_connections({{1, 2}, {1, 3}}).use_topology(true); + + if (use_nop_partitions) wick.set_nop_partitions({{2, 3}}); + if (use_op_partitions) + wick.set_op_partitions({{0, 1}, + {2, 3}, + {4, 5}, + {6, 7}, + {8, 9}, + {10, 11}, + {12, 13}, + {14, 15}}); + auto wick_result = wick.compute(); + REQUIRE(wick_result->is()); + if (use_op_partitions) { + REQUIRE(wick_result->size() == 7); + } else { + REQUIRE(wick_result->size() == 544); + } + + // multiply tensor factors and expand + auto wick_result_2 = + ex(rational{1, 256}) * + ex(L"A", IndexList{L"i_1", L"i_2"}, + IndexList{{L"a_1", {L"i_1", L"i_2"}}, + {L"a_2", {L"i_1", L"i_2"}}}, + IndexList{}, Symmetry::antisymm) * + ex(L"g", WstrList{L"p_1", L"p_2"}, + WstrList{L"p_3", L"p_4"}, WstrList{},Symmetry::antisymm) * + ex(L"t", + IndexList{{L"a_3", {L"i_3", L"i_4"}}, + {L"a_4", {L"i_3", L"i_4"}}}, + IndexList{L"i_3", L"i_4"}, IndexList{}, Symmetry::antisymm) * + ex(L"t", + IndexList{{L"a_5", {L"i_5", L"i_6"}}, + {L"a_6", {L"i_5", L"i_6"}}}, + IndexList{L"i_5", L"i_6"}, IndexList{}, Symmetry::antisymm) * + wick_result; + expand(wick_result_2); + wick.reduce(wick_result_2); + rapid_simplify(wick_result_2); + TensorCanonicalizer::register_instance( + std::make_shared()); + canonicalize(wick_result_2); + canonicalize(wick_result_2); + canonicalize(wick_result_2); + rapid_simplify(wick_result_2); + + // std::wcout << oss.str() << to_latex_align(wick_result_2, 20) + // << std::endl; + REQUIRE(wick_result_2->is()); + REQUIRE(wick_result_2->size() == 4); + } // use_op_partitions + } // use_nop_partitions + }); #if 1 - // 3-body ^ 2-body ^ 2-body ^ 3-body - SEQUANT_PROFILE_SINGLE("wick(P3*H2*T2*T3)", { - constexpr bool connected_only = true; - constexpr bool topology = true; - auto P3 = ex(rational{1, 36}) * - ex(L"A", WstrList{L"i_1", L"i_2", L"i_3"}, - WstrList{L"a_1", L"a_2", L"a_3"}, WstrList{}, - Symmetry::antisymm) * - ex(WstrList{L"i_1", L"i_2", L"i_3"}, - WstrList{L"a_1", L"a_2", L"a_3"}); - auto H2 = - ex(rational{1, 4}) * - ex(L"g", WstrList{L"p_1", L"p_2"}, WstrList{L"p_3", L"p_4"}, - WstrList{}, Symmetry::antisymm) * - ex(WstrList{L"p_1", L"p_2"}, WstrList{L"p_3", L"p_4"}); - auto T2 = - ex(rational{1, 4}) * - ex(L"t", WstrList{L"a_4", L"a_5"}, WstrList{L"i_4", L"i_5"}, - WstrList{}, Symmetry::antisymm) * - ex(WstrList{L"a_4", L"a_5"}, WstrList{L"i_4", L"i_5"}); - auto T3 = ex(rational{1, 36}) * - ex(L"t", WstrList{L"a_6", L"a_7", L"a_8"}, - WstrList{L"i_6", L"i_7", L"i_8"}, WstrList{}, - Symmetry::antisymm) * - ex(WstrList{L"a_6", L"a_7", L"a_8"}, - WstrList{L"i_6", L"i_7", L"i_8"}); - FWickTheorem wick{P3 * H2 * T2 * T3}; - wick.use_topology(topology); - if (connected_only) wick.set_nop_connections({{1, 2}, {1, 3}}); - auto wick_result = wick.compute(); - - std::wcout << "P3*H2*T2*T3 = " << to_latex_align(wick_result, 20) - << std::endl; - REQUIRE(wick_result->is()); - REQUIRE( - wick_result->size() == - (connected_only ? 7 : 9)); // 9 = 2 disconnected + 7 connected terms - }); + // 3-body ^ 2-body ^ 2-body ^ 3-body + SEQUANT_PROFILE_SINGLE("wick(P3*H2*T2*T3)", { + constexpr bool connected_only = true; + constexpr bool topology = true; + auto P3 = + ex(rational{1, 36}) * + ex(L"A", WstrList{L"i_1", L"i_2", L"i_3"}, + WstrList{L"a_1", L"a_2", L"a_3"}, WstrList{}, Symmetry::antisymm) * + ex(WstrList{L"i_1", L"i_2", L"i_3"}, + WstrList{L"a_1", L"a_2", L"a_3"}); + auto H2 = + ex(rational{1, 4}) * + ex(L"g", WstrList{L"p_1", L"p_2"}, WstrList{L"p_3", L"p_4"}, + WstrList{}, Symmetry::antisymm) * + ex(WstrList{L"p_1", L"p_2"}, WstrList{L"p_3", L"p_4"}); + auto T2 = + ex(rational{1, 4}) * + ex(L"t", WstrList{L"a_4", L"a_5"}, WstrList{L"i_4", L"i_5"}, + WstrList{}, Symmetry::antisymm) * + ex(WstrList{L"a_4", L"a_5"}, WstrList{L"i_4", L"i_5"}); + auto T3 = + ex(rational{1, 36}) * + ex(L"t", WstrList{L"a_6", L"a_7", L"a_8"}, + WstrList{L"i_6", L"i_7", L"i_8"}, WstrList{}, Symmetry::antisymm) * + ex(WstrList{L"a_6", L"a_7", L"a_8"}, + WstrList{L"i_6", L"i_7", L"i_8"}); + FWickTheorem wick{P3 * H2 * T2 * T3}; + wick.use_topology(topology); + if (connected_only) wick.set_nop_connections({{1, 2}, {1, 3}}); + auto wick_result = wick.compute(); + + std::wcout << "P3*H2*T2*T3 = " << to_latex_align(wick_result, 20) + << std::endl; + REQUIRE(wick_result->is()); + REQUIRE( + wick_result->size() == + (connected_only ? 7 : 9)); // 9 = 2 disconnected + 7 connected terms + }); #endif -} + } } // TEST_CASE("WickTheorem") #endif