Problems with Operator::erf_nuclear/Operator::erfc_nuclear

[maxscheurer]

While adding Operator::erf_nuclear/Operator::erfc_nuclear integrals to Psi4 (psi4/psi4#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, #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$