Problems with Operator::erf_nuclear/Operator::erfc_nuclear #242
Description
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!