Merge branch 'master' into kshitij/feature/dfbs_generator

evaleev · Nov 14, 2023 · 6f0583a · 6f0583a
2 parents f399f09 + 43bec65
commit 6f0583a
Show file tree

Hide file tree

Showing 3 changed files with 70 additions and 3 deletions.
diff --git a/include/libint2/engine.impl.h b/include/libint2/engine.impl.h
@@ -1054,7 +1054,6 @@ __libint2_engine_inline void Engine::compute_primdata(Libint_t& primdata, const
                      primdata.PC_y[0] * primdata.PC_y[0] +
                      primdata.PC_z[0] * primdata.PC_z[0];
     const scalar_type U = gammap * PC2;
-    const scalar_type rho = rhop;
     const auto mmax = s1.contr[0].l + s2.contr[0].l + deriv_order_;
     auto* fm_ptr = &(primdata.LIBINT_T_S_ELECPOT_S(0)[0]);
     if (oper_ == Operator::nuclear) {
@@ -1069,15 +1068,15 @@ __libint2_engine_inline void Engine::compute_primdata(Libint_t& primdata, const
       const auto& core_ints_params =
           std::get<0>(any_cast<const typename operator_traits<
             Operator::erf_nuclear>::oper_params_type&>(core_ints_params_));
-      core_eval_ptr->eval(fm_ptr, rho, U, mmax, core_ints_params);
+      core_eval_ptr->eval(fm_ptr, gammap, U, mmax, core_ints_params);
     } else if (oper_ == Operator::erfc_nuclear) {
       const auto& core_eval_ptr =
           any_cast<const detail::core_eval_pack_type<Operator::erfc_nuclear>&>(core_eval_pack_)
             .first();
       const auto& core_ints_params =
           std::get<0>(any_cast<const typename operator_traits<
             Operator::erfc_nuclear>::oper_params_type&>(core_ints_params_));
-      core_eval_ptr->eval(fm_ptr, rho, U, mmax, core_ints_params);
+      core_eval_ptr->eval(fm_ptr, gammap, U, mmax, core_ints_params);
     }
 
     decltype(U) two_o_sqrt_PI(1.12837916709551257389615890312);

diff --git a/tests/hartree-fock/hartree-fock.cc b/tests/hartree-fock/hartree-fock.cc
@@ -651,6 +651,9 @@ Matrix compute_2body_fock_simple(const std::vector<libint2::Shell>& shells,
           // exchange contribution to the Fock matrix is from {s1,s3,s2,s4} integrals
           engine.compute(shells[s1], shells[s3], shells[s2], shells[s4]);
           const auto* buf_1324 = buf[0];
+          if (buf_1324 == nullptr)
+            continue; // if all integrals screened out, skip to next quartet
+
 
           for(auto f1=0, f1324=0; f1!=n1; ++f1) {
             const auto bf1 = f1 + bf1_first;

diff --git a/tests/unit/test-1body.cc b/tests/unit/test-1body.cc
@@ -76,3 +76,68 @@ TEST_CASE_METHOD(libint2::unit::DefaultFixture, "electrostatic potential", "[eng
 
 #endif // LIBINT2_SUPPORT_ONEBODY
 }
+
+TEST_CASE_METHOD(libint2::unit::DefaultFixture, "erf correctness", "[engine][1-body]") {
+#if defined(LIBINT2_SUPPORT_ONEBODY)
+  if (LIBINT_SHGSHELL_ORDERING != LIBINT_SHGSHELL_ORDERING_STANDARD)
+    return;
+
+  std::vector<Shell> obs{
+      Shell{{1.0, 3.0}, {{2, true, {1.0, 0.3}}}, {{0.0, 0.0, 0.0}}},
+      Shell{{2.0, 5.0}, {{2, true, {1.0, 0.2}}}, {{1.0, 1.0, 1.0}}}};
+  {
+    const auto lmax = std::min(3, LIBINT2_MAX_AM_elecpot);
+    if (lmax >= 2) {
+      auto engine = Engine(Operator::nuclear, 2, lmax);
+      engine.set_params(make_point_charges(atoms));
+
+      const auto scale = 2.3;
+      engine.prescale_by(scale);
+      engine.compute(obs[0], obs[0]);
+      {
+        // All reference values are pulled from L2.
+        std::vector<double> shellset_ref = {
+            -1.238239259091998e+01, 0.000000000000000e+00,
+            0.000000000000000e+00,  -5.775996163160049e-02,
+            0.000000000000000e+00,  0.000000000000000e+00,
+            -1.301230978657952e+01, -6.796143730068988e-02,
+            0.000000000000000e+00,  1.139389632827834e-01,
+            0.000000000000000e+00,  -6.796143730068988e-02,
+            -1.343732979153083e+01, 0.000000000000000e+00,
+            -1.478824785355970e-02, -5.775996163160049e-02,
+            0.000000000000000e+00,  0.000000000000000e+00,
+            -1.284475452992947e+01, 0.000000000000000e+00,
+            0.000000000000000e+00,  1.139389632827834e-01,
+            -1.478824785355970e-02, 0.000000000000000e+00,
+            -1.241040347301479e+01};
+        for (int i = 0; i != 25; ++i) {
+          REQUIRE(engine.results()[0][i]/scale == Approx(shellset_ref[i]));
+        }
+      }
+
+      engine.prescale_by(1);
+      engine.compute(obs[0], obs[1]);
+      {
+        std::vector<double> shellset_ref = {
+            -4.769186621041819e-01, -9.303619356400431e-01,
+            -1.559058302243514e+00, -9.290824121864600e-01,
+            -5.835786921473129e-04, -1.159266418436018e+00,
+            -3.770080831197964e-01, 9.572841308198474e-01,
+            -8.291498398421207e-01, -1.663667687168316e+00,
+            -2.171951144148577e+00, 1.074249956874296e+00,
+            2.128355904665372e+00,  1.074590109905394e+00,
+            -3.485163651594458e-03, -1.160865205880651e+00,
+            -8.344173649626901e-01, 9.566621490332916e-01,
+            -3.760919234260182e-01, 1.660514988916377e+00,
+            -1.120272634615116e-03, -1.385603731947886e+00,
+            -2.105750177166632e-03, 1.380654897976564e+00,
+            2.115041199099945e+00};
+        for (int i = 0; i != 25; ++i) {
+          REQUIRE(engine.results()[0][i] == Approx(shellset_ref[i]));
+        }
+      }
+    }
+  }
+#endif // LIBINT2_SUPPORT_ONEBODY
+}
+