Merge pull request #236 from ValeevGroup/ajay/fix/projection-with-zer…

…o-rank Fix projection with zero ranks
ValeevGroup · Sep 6, 2024 · 2221798 · 2221798
2 parents 1fd9e6f + 9d86156
commit 2221798
Show file tree

Hide file tree

Showing 2 changed files with 22 additions and 7 deletions.
diff --git a/SeQuant/domain/mbpt/op.cpp b/SeQuant/domain/mbpt/op.cpp
@@ -618,7 +618,7 @@ ExprPtr P(nₚ np, nₕ nh) {
 
 ExprPtr A(nₚ np, nₕ nh) {
   assert(!(np == 0 && nh == 0));
-  // if they are not zero, Kh and Kp should have the same sign
+  // if one of them is not zero, nh and np should have the same sign
   if (np != 0 && nh != 0) {
     assert((np > 0 && nh > 0) || (np < 0 && nh < 0));
   }
@@ -628,20 +628,21 @@ ExprPtr A(nₚ np, nₕ nh) {
   if (nh > 0)  // ex
     for ([[maybe_unused]] auto i : ranges::views::iota(0, nh))
       annihilators.emplace_back(get_hole_space(Spin::any));
-  else  // deex
+  else if (nh < 0)  // deex
     for ([[maybe_unused]] auto i : ranges::views::iota(0, -nh))
       creators.emplace_back(get_hole_space(Spin::any));
   if (np > 0)  // ex
     for ([[maybe_unused]] auto i : ranges::views::iota(0, np))
       creators.emplace_back(get_particle_space(Spin::any));
-  else  // deex
+  else if (np < 0)  // deex
     for ([[maybe_unused]] auto i : ranges::views::iota(0, -np))
       annihilators.emplace_back(get_particle_space(Spin::any));
+  // don't populate if rank is zero
 
   std::optional<OpMaker<Statistics::FermiDirac>::UseDepIdx> dep;
   if (get_default_mbpt_context().csv() == mbpt::CSV::Yes)
-    dep = np > 0 ? OpMaker<Statistics::FermiDirac>::UseDepIdx::Bra
-                 : OpMaker<Statistics::FermiDirac>::UseDepIdx::Ket;
+    dep = (np > 0 || nh > 0) ? OpMaker<Statistics::FermiDirac>::UseDepIdx::Bra
+                             : OpMaker<Statistics::FermiDirac>::UseDepIdx::Ket;
   return OpMaker<Statistics::FermiDirac>(
       OpType::A, cre(creators), ann(annihilators))(dep, {Symmetry::antisymm});
 }
@@ -815,10 +816,12 @@ ExprPtr F(bool use_f_tensor, IndexSpace occupied_density) {
 
 ExprPtr A(nₚ np, nₕ nh) {
   assert(!(nh == 0 && np == 0));
-  // if they are not zero, Kh and Kp should have the same sign
+  // if one of them is not zero, nh and np should have the same sign
   if (nh != 0 && np != 0) {
     assert((nh > 0 && np > 0) || (nh < 0 && np < 0));
   }
+  // if np or nh is negative, it's a deexcitation operator
+  const auto deexcitation = (np < 0 || nh < 0);
 
   auto particle_space = get_particle_space(Spin::any);
   auto hole_space = get_hole_space(Spin::any);
@@ -827,7 +830,7 @@ ExprPtr A(nₚ np, nₕ nh) {
                   [=](qnc_t& qns) {
                     const std::size_t abs_nh = std::abs(nh);
                     const std::size_t abs_np = std::abs(np);
-                    if (np < 0) {
+                    if (deexcitation) {
                       qnc_t op_qnc_t = generic_deexcitation_qns(
                           abs_np, abs_nh, particle_space, hole_space);
                       qns = combine(op_qnc_t, qns);

diff --git a/tests/unit/test_mbpt.cpp b/tests/unit/test_mbpt.cpp
@@ -288,6 +288,18 @@ TEST_CASE("NBodyOp", "[mbpt]") {
     auto lambda2_f = Λ_(2) * H_(1);
     REQUIRE(lowers_rank_to_vacuum(lambda2_f, 2));
 
+    auto expr1 = P(nₚ(0), nₕ(1)) * H() * R(nₚ(0), nₕ(1));
+    auto expr1_tnsr = lower_to_tensor_form(expr1);
+    auto vev1_op = op::vac_av(expr1);
+    auto vev1_t = tensor::vac_av(expr1_tnsr);  // no operator level screening
+    REQUIRE(to_latex(vev1_op) == to_latex(vev1_t));
+
+    auto expr2 = P(nₚ(2), nₕ(1)) * H() * R(nₚ(1), nₕ(0));
+    auto expr2_tnsr = lower_to_tensor_form(expr2);
+    auto vev2_op = op::vac_av(expr2);
+    auto vev2_t = tensor::vac_av(expr2_tnsr);  // no operator level screening
+    REQUIRE(to_latex(vev2_op) == to_latex(vev2_t));
+
   }  // SECTION("screen")
 
   SECTION("operators") {