.. _program_listing_file_SeQuant_domain_mbpt_models_cc.cpp:

Program Listing for File cc.cpp
===============================

|exhale_lsh| :ref:`Return to documentation for file <file_SeQuant_domain_mbpt_models_cc.cpp>` (``SeQuant/domain/mbpt/models/cc.cpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   #include <SeQuant/core/expr.hpp>
   #include <SeQuant/core/rational.hpp>
   #include <SeQuant/core/runtime.hpp>
   #include <SeQuant/domain/mbpt/context.hpp>
   #include <SeQuant/domain/mbpt/convention.hpp>
   #include <SeQuant/domain/mbpt/models/cc.hpp>
   #include <SeQuant/domain/mbpt/op.hpp>
   #include <SeQuant/domain/mbpt/spin.hpp>
   
   #include <cassert>
   #include <cstdint>
   #include <memory>
   #include <new>
   #include <stdexcept>
   #include <utility>
   
   namespace sequant::mbpt {
   
   CC::CC(size_t n, Ansatz a) : N(n), ansatz_(a) {}
   
   CC::Ansatz CC::ansatz() const { return ansatz_; }
   
   bool CC::unitary() const {
     return ansatz_ == Ansatz::U || ansatz_ == Ansatz::oU;
   }
   
   ExprPtr CC::sim_tr(ExprPtr expr, size_t commutator_rank) {
     const bool skip_singles = ansatz_ == Ansatz::oT || ansatz_ == Ansatz::oU;
     auto transform_op_op_pdt = [this, &commutator_rank,
                                 skip_singles](const ExprPtr& expr) {
       assert(expr.is<op_t>() || expr.is<Product>());
       auto result = expr;
       auto op_Sk = result;
       for (size_t k = 1; k <= commutator_rank; ++k) {
         ExprPtr op_Sk_comm_w_S;
         op_Sk_comm_w_S =
             op_Sk * T(N, skip_singles);  // traditional SR ansatz: [O,T] = (O T)_c
         if (unitary())  // unitary SR ansatz: [O,T-T^+] = (O T)_c + (T^+ O)_c
           op_Sk_comm_w_S += adjoint(T(N, skip_singles)) * op_Sk;
         op_Sk = ex<Constant>(rational{1, k}) * op_Sk_comm_w_S;
         simplify(op_Sk);
         result += op_Sk;
       }
       return result;
     };
   
     if (expr.is<op_t>()) {
       return transform_op_op_pdt(expr);
     } else if (expr.is<Product>()) {
       auto& product = expr.as<Product>();
       // Expand product as sum
       if (ranges::any_of(product.factors(), [](const auto& factor) {
             return factor.template is<Sum>();
           })) {
         expr = sequant::expand(expr);
         simplify(expr);
         return sim_tr(expr, commutator_rank);
       } else {
         return transform_op_op_pdt(expr);
       }
     } else if (expr.is<Sum>()) {
       auto result = sequant::transform_reduce(
           *expr, ex<Sum>(),
           [](const ExprPtr& running_total, const ExprPtr& summand) {
             return running_total + summand;
           },
           [=](const auto& op_product) {
             return transform_op_op_pdt(op_product);
           });
       return result;
     } else if (expr.is<Constant>() || expr.is<Variable>())
       return expr;
     else
       throw std::invalid_argument(
           "CC::sim_tr(expr): Unsupported expression type");
   }
   
   std::vector<ExprPtr> CC::t(size_t commutator_rank, size_t pmax, size_t pmin) {
     pmax = (pmax == std::numeric_limits<size_t>::max() ? N : pmax);
   
     assert(commutator_rank >= 1 && "commutator rank should be >= 1");
     assert(pmax >= pmin && "pmax should be >= pmin");
   
     // 1. construct hbar(op) in canonical form
     auto hbar = sim_tr(H(), commutator_rank);
   
     // 2. project onto each manifold, screen, lower to tensor form and wick it
     std::vector<ExprPtr> result(pmax + 1);
     for (std::int64_t p = pmax; p >= static_cast<std::int64_t>(pmin); --p) {
       // 2.a. screen out terms that cannot give nonzero after projection onto
       // <p|
       std::shared_ptr<Sum>
           hbar_p;  // products that can produce excitations of rank p
       std::shared_ptr<Sum>
           hbar_le_p;  // keeps products that can produce excitations rank <=p
       for (auto& term : *hbar) {
         assert(term->is<Product>() || term->is<op_t>());
         if (raises_vacuum_up_to_rank(term, p)) {
           if (!hbar_le_p)
             hbar_le_p = std::make_shared<Sum>(ExprPtrList{term});
           else
             hbar_le_p->append(term);
           if (raises_vacuum_to_rank(term, p)) {
             if (!hbar_p)
               hbar_p = std::make_shared<Sum>(ExprPtrList{term});
             else
               hbar_p->append(term);
           }
         }
       }
       hbar = hbar_le_p;
       // 2.b project onto <p| (i.e., multiply by P(p) if p>0) and compute VEV
       result.at(p) = vac_av(p != 0 ? P(nₚ(p)) * hbar_p : hbar_p);
     }
   
     return result;
   }
   
   std::vector<ExprPtr> CC::λ(size_t commutator_rank) {
     assert(commutator_rank >= 1 && "commutator rank should be >= 1");
     assert(!unitary() && "there is no need for CC::λ for unitary ansatz");
   
     // construct hbar
     auto hbar = sim_tr(H(), commutator_rank - 1);
   
     const auto One = ex<Constant>(1);
     auto lhbar = simplify((One + Λ(N)) * hbar);
   
     const auto op_connect =
         concat(default_op_connections(),
                std::vector<std::pair<OpType, OpType>>{{OpType::h, OpType::A},
                                                       {OpType::f, OpType::A},
                                                       {OpType::g, OpType::A},
                                                       {OpType::h, OpType::S},
                                                       {OpType::f, OpType::S},
                                                       {OpType::g, OpType::S}});
   
     // 2. project onto each manifold, screen, lower to tensor form and wick it
     std::vector<ExprPtr> result(N + 1);
     for (auto p = N; p >= 1; --p) {
       // 2.a. screen out terms that cannot give nonzero after projection onto
       // <P|
       std::shared_ptr<Sum>
           hbar_p;  // products that can produce excitations of rank p
       std::shared_ptr<Sum>
           hbar_le_p;  // keeps products that can produce excitations rank <=p
       for (auto& term : *lhbar) {  // pick terms from lhbar
         assert(term->is<Product>() || term->is<op_t>());
   
         if (lowers_rank_or_lower_to_vacuum(term, p)) {
           if (!hbar_le_p)
             hbar_le_p = std::make_shared<Sum>(ExprPtrList{term});
           else
             hbar_le_p->append(term);
           if (lowers_rank_to_vacuum(term, p)) {
             if (!hbar_p)
               hbar_p = std::make_shared<Sum>(ExprPtrList{term});
             else
               hbar_p->append(term);
           }
         }
       }
       lhbar = hbar_le_p;
   
       // 2.b multiply by adjoint of P(p) (i.e., P(-p)) on the right side and
       // compute VEV
       result.at(p) = vac_av(hbar_p * P(nₚ(-p)), op_connect);
     }
     return result;
   }
   
   std::vector<ExprPtr> CC::t_pt(size_t order, size_t rank) {
     assert(order == 1 &&
            "sequant::mbpt::CC::t_pt(): only first-order perturbation is "
            "supported now");
     assert(rank == 1 &&
            "sequant::mbpt::CC::t_pt(): only one-body perturbation "
            "operator is supported now");
     assert(ansatz_ == Ansatz::T && "unitary ansatz is not yet supported");
   
     // construct h1_bar
   
     // truncate h1_bar at rank 2 for one-body perturbation
     // operator and at rank 4 for two-body perturbation operator
     const auto h1_truncate_at = rank == 1 ? 2 : 4;
   
     auto h1_bar = sim_tr(H_pt(1, rank), h1_truncate_at);
   
     // construct [hbar, T(1)]
     auto hbar_pert = sim_tr(H(), 3) * T_pt(order, N);
   
     // [Eq. 34, WIREs Comput Mol Sci. 2019; 9:e1406]
     auto expr = simplify(h1_bar + hbar_pert);
   
     // connectivity:
     // connect t and t1 with {h,f,g}
     // connect h1 with t
     const auto op_connect =
         concat(default_op_connections(),
                std::vector<std::pair<OpType, OpType>>{{OpType::h, OpType::t_1},
                                                       {OpType::f, OpType::t_1},
                                                       {OpType::g, OpType::t_1},
                                                       {OpType::h_1, OpType::t}});
   
     std::vector<ExprPtr> result(N + 1);
     for (auto p = N; p >= 1; --p) {
       auto freq_term = ex<Variable>(L"ω") * P(nₚ(p)) * T_pt_(order, p);
       result.at(p) = vac_av(P(nₚ(p)) * expr, op_connect) - vac_av(freq_term);
     }
     return result;
   }
   
   std::vector<ExprPtr> CC::λ_pt(size_t order, size_t rank) {
     assert(order == 1 &&
            "sequant::mbpt::CC::λ_pt(): only first-order perturbation is "
            "supported now");
     assert(rank == 1 &&
            "sequant::mbpt::CC::λ_pt(): only one-body perturbation "
            "operator is supported now");
     assert(ansatz_ == Ansatz::T && "unitary ansatz is not yet supported");
   
     // construct hbar
     auto hbar = sim_tr(H(), 4);
   
     // construct h1_bar
   
     // truncate h1_bar at rank 2 for one-body perturbation
     // operator and at rank 4 for two-body perturbation operator
     const auto h1_truncate_at = rank == 1 ? 2 : 4;
   
     auto h1_bar = sim_tr(H_pt(1, rank), h1_truncate_at);
     // construct [hbar, T(1)]
     auto hbar_pert = sim_tr(H(), 3) * T_pt(order, N);
   
     // [Eq. 35, WIREs Comput Mol Sci. 2019; 9:e1406]
     const auto One = ex<Constant>(1);
     auto expr =
         simplify((One + Λ(N)) * (h1_bar + hbar_pert) + Λ_pt(order, N) * hbar);
   
     // connectivity:
     // t and t1 with {h,f,g}
     // projectors with {h,f,g}
     // h1 with t
     // h1 with projectors
     const auto op_connect =
         concat(default_op_connections(),
                std::vector<std::pair<OpType, OpType>>{{OpType::h, OpType::t_1},
                                                       {OpType::f, OpType::t_1},
                                                       {OpType::g, OpType::t_1},
                                                       {OpType::h_1, OpType::t},
                                                       {OpType::h, OpType::A},
                                                       {OpType::f, OpType::A},
                                                       {OpType::g, OpType::A},
                                                       {OpType::h, OpType::S},
                                                       {OpType::f, OpType::S},
                                                       {OpType::g, OpType::S},
                                                       {OpType::h_1, OpType::A},
                                                       {OpType::h_1, OpType::S}});
   
     std::vector<ExprPtr> result(N + 1);
     for (auto p = N; p >= 1; --p) {
       auto freq_term = ex<Variable>(L"ω") * Λ_pt_(order, p) * P(nₚ(-p));
       result.at(p) = vac_av(expr * P(nₚ(-p)), op_connect) + vac_av(freq_term);
     }
     return result;
   }
   
   std::vector<ExprPtr> CC::eom_r(nₚ np, nₕ nh) {
     assert(!unitary() && "Unitary ansatz is not yet supported");
     assert(np > 0 || nh > 0 && "Unsupported excitation order");
     assert(np == nh &&
            "Only EE-EOM-CC has been tested ... remove this assert to try "
            "Fock-space EOM-CC");
   
     if (np != nh)
       assert(
           get_default_context().spbasis() != SPBasis::spinfree &&
           "spin-free basis does not yet support non particle-conserving cases");
   
     // construct hbar
     auto hbar = sim_tr(H(), 4);
   
     // hbar * R
     auto hbar_R = hbar * R(np, nh);
   
     // connectivity:
     // default connections + connect R with {h,f,g}
     const auto op_connect =
         concat(default_op_connections(),
                std::vector<std::pair<OpType, OpType>>{{OpType::h, OpType::R},
                                                       {OpType::f, OpType::R},
                                                       {OpType::g, OpType::R}});
   
     // initialize result vector
     std::vector<ExprPtr> result;
     using std::max;
     auto idx = max(np, nh);  // index for populating the result vector
     result.resize(idx + 1);
   
     // start from the highest excitation order, go down to the lowest possible
     for (std::int64_t rp = np, rh = nh; rp > 0 || rh > 0; --rp, --rh) {
       // project with <rp, rh| (i.e., multiply P(rp, rh)) and compute VEV
       result.at(idx) = vac_av(P(nₚ(rp), nₕ(rh)) * hbar_R, op_connect);
       idx--;  // index decrement
     }
   
     return result;
   }
   
   std::vector<ExprPtr> CC::eom_l(nₚ np, nₕ nh) {
     assert(!unitary() && "Unitary ansatz is not yet supported");
     assert(np > 0 || nh > 0 && "Unsupported excitation order");
     assert(np == nh &&
            "Only EE-EOM-CC has been tested ... remove this assert to try "
            "Fock-space EOM-CC");
   
     if (np != nh)
       assert(get_default_context().spbasis() != SPBasis::spinfree &&
              "spin-free basis does not support non particle-conserving cases");
   
     // construct hbar
     auto hbar = sim_tr(H(), 4);
   
     // L * hbar
     auto L_hbar = L(np, nh) * hbar;
   
     // connectivity:
     // default connections + connect H with projectors
     const auto op_connect =
         concat(default_op_connections(),
                std::vector<std::pair<OpType, OpType>>{{OpType::h, OpType::A},
                                                       {OpType::f, OpType::A},
                                                       {OpType::g, OpType::A},
                                                       {OpType::h, OpType::S},
                                                       {OpType::f, OpType::S},
                                                       {OpType::g, OpType::S}});
   
     // initialize result vector
     std::vector<ExprPtr> result;
     using std::max;
     auto idx = max(np, nh);  // index for populating the result vector
     result.resize(idx + 1);
   
     // start from the highest excitation order, go down to the lowest possible
     for (std::int64_t rp = np, rh = nh; rp > 0 || rh > 0; --rp, --rh) {
       // right project with |rp,rh> (i.e., multiply P(-np, -rh)) and compute VEV
       result.at(idx) = vac_av(L_hbar * P(nₚ(-rp), nₕ(-rh)), op_connect);
       idx--;  // index decrement
     }
   
     return result;
   }
   }  // namespace sequant::mbpt