Problems with `Operator::erf_nuclear`/`Operator::erfc_nuclear`

[maxscheurer]

While adding Operator::erf_nuclear/Operator::erfc_nuclear integrals to Psi4 (https://github.com/psi4/psi4/pull/2473), I noticed that the results from libint2 don't agree with what WolframAlpha (don't have a Mathematica license, sorry 😄) gives me.

Here's a minimal example (Hydrogen atom with a single primitive) to reproduce this (the Operator::nuclear part is just to verify that my setup is correct):

#include <array>
#include <iomanip>
#include <libint2/engine.h>
#include <memory>
#include <stdio.h>
#include <vector>

using Zxyz_vector = std::vector<std::pair<double, std::array<double, 3>>>;

int main() {
  libint2::initialize();
  auto s = libint2::Shell{{3.42525091},
                          {
                                {0, false, {1.794441832218435}},
                          },
                          {{0.0, 0.0, 0.0}}};

  {
    std::unique_ptr<libint2::Engine> engine0_ =
          std::make_unique<libint2::Engine>(libint2::Operator::nuclear, 1, 0, 0);
    Zxyz_vector pcs;
    pcs.push_back({-1.0, {0.0, 0.0, 0.0}});
    engine0_->set_params(pcs);
    engine0_->compute(s, s);

    auto buf = engine0_->results()[0];
    std::cout << std::setprecision(14) << buf[0] << std::endl;
  }

  {
    std::unique_ptr<libint2::Engine> engine0_ =
          std::make_unique<libint2::Engine>(libint2::Operator::erf_nuclear, 1, 0, 0);
    Zxyz_vector pcs;
    pcs.push_back({-1.0, {0.0, 0.0, 0.0}});
    std::tuple<double, Zxyz_vector> params{1.0, pcs};
    engine0_->set_params(params);
    engine0_->compute(s, s);

    auto buf = engine0_->results()[0];
    std::cout << std::setprecision(14) << buf[0] << std::endl;
  }

  {
    std::unique_ptr<libint2::Engine> engine0_ =
          std::make_unique<libint2::Engine>(libint2::Operator::erfc_nuclear, 1, 0, 0);
    Zxyz_vector pcs;
    pcs.push_back({-1.0, {0.0, 0.0, 0.0}});
    std::tuple<double, Zxyz_vector> params{1.0, pcs};
    engine0_->set_params(params);
    engine0_->compute(s, s);

    auto buf = engine0_->results()[0];
    std::cout << std::setprecision(14) << buf[0] << std::endl;
  }

  libint2::finalize();
  return 0;
}

The output from the above code is:

2.9533590737505  // nuclear
1.7931694693882  // erf
1.1601896043623  //erfc

However, from WolframAlpha, I get 1.05407 (erf) and 1.89929 (erfc), which I could also reproduce using a hand-written McMurchie-Davidson code (thanks to @andysim).

Now, when scaling omega by a factor of 1/2 in libint2, i.e. std::tuple<double, Zxyz_vector> params{0.5, pcs};, I get the "correct" result for the minimal example above. Some more diagnostics:

• ✅ the "scaling-by-1/2-fix" works for increased angular momentum (i.e., single p primitive)
• ✅ the "fix" also works for multiple Hydrogen atoms, each bearing a single primitive
• 🚫 the "fix" does NOT work when I have multiple primitives with different exponents, e.g., H/sto-3g

From the list above, we (@andysim and I) concluded that the discrepancy a) does not depend on the basis set angular momentum nor b) the location of the basis functions. ~We think it's most likely related to the way the Boys function for these integrals is evaluated (because of the primitive exponents...).~

Any hint on how to resolve this would be greatly appreciated. Thanks a lot!

[maxscheurer]

P.S.: Just as "proof" a screenshot of the WolframAlpha result:

[loriab]

Since it's possible that Max and Andy are working from my 2019 fork of L2 for cmake, I just re-ran the test file on #233 that has been re-synced to master this year. Same result:

(basenol2) psilocaluser@bash:psinet:/psi/gits/libint2-efv/sandbox-erfc: (new-cmake-harness-lab-rb6) ./build/minerfc 
2.9533590737505
1.7931694693882
1.1601896043623

[JonathonMisiewicz]

Since you suspect an error in the Boys function, #243 may fix this.

On #243 merge, I suggest we re-sync #233 and have Lori (sorry) re-run the test file.

[maxscheurer]

Since you suspect an error in the Boys function, https://github.com/evaleev/libint/pull/243 may fix this.

Only if FmEval_Taylor is used. I'm not sure but I think the default is FmEval_Chebyshev7...

[maxscheurer]

I'm quite certain now that erf_coulomb_gm_eval/erfc_coulomb_gm_eval are correct, otherwise, the erf_coulomb/erfc_coulomb two-electron integrals would be wrong, too. Any hints on how to get to the bottom of this, @evaleev?

[loriab]

fwiw, no cure on this from #259, which is up-to-date with present master.

2.9533590737505
1.7931694693882
1.1601896043623

[JonathonMisiewicz]

The diff if anybody needs this fixed right now is

diff --git a/include/libint2/engine.impl.h b/include/libint2/engine.impl.h
index a1ee9c1c..353c0cd3 100644
--- 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,7 +1068,7 @@ __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_)
@@ -1077,7 +1076,7 @@ __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::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);

L2 still disagrees with Mathematica in the eighth decimal place for the test case, after this. That tells me that I probably need to learn what knobs to push to get higher precision out of L2.

[evaleev]

@JonathonMisiewicz nice catch! I think default precision in Libint should be sufficient to get 14+ digits unless you are dealing with extreme exponents (or other params). Just double check your coordinates/exponents match exactly.

[JonathonMisiewicz]

I stand corrected. I can get 10 decimal places of agreement, and Mathematica breaks when I tell it to do the numerical integral tightly enough to check further. L2 is victorious.

I can either do a lightning-quick PR tomorrow to close this, or I can add a test to make sure this stays fixed. If you want a test, let me know what template to follow.

[evaleev]

@JonathonMisiewicz if you could add a unit test to https://github.com/evaleev/libint/blob/master/tests/unit/test-1body.cc would be good ... just follow the patterns

[evaleev]

for posterity: 1-body (point charge) integral formulas in DOI 10.1039/b605188j are obtained by substituting $\lim{\eta \to \infty} \rho \equiv \lim{\eta \to \infty} \frac{\zeta \eta}{\zeta + \eta} = \lim_{\eta \to \infty} \frac{\zeta}{\frac{\zeta}{\eta}+1} = \zeta$