Assertion failure crash while computing 3-center ERI derivatives

[smparker]

Hi there,

I've run into a crash caused by what looks like an internal assertion failing. I'm still not sure if it's a bug or if I'm just misunderstanding the interface, so I'm hoping you can help figure it out. I'm trying to compute the first derivative of the 3-center ERIs using the c++11 interface. The basic error comes from line 1745 of my engine.impl.h (I realize the line number might change with configure options) and the failing assertion reads buildfnptrs_[buildfnidx] && "null build function ptr".

I have been able to make a minimum example that shows the same behavior that I'll paste below. But first, details.

• I'm using v2.7.1 with a fairly standard configuration (I can share if that'd help)
• It only seems to effect 1st derivatives of 3c integrals.
  • regular 3c integrals (no derivative) seem to work fine with Braket::xs_xx
  • 2c integrals and their derivatives seem to work fine with Braket::xs_xs
• It doesn't seem to matter whether I use the engine.compute() or engine.compute2() interfaces
• It always crashes at the first p function the code encounters. Increasing max_l doesn't seem to change this at all

The test code that triggers this is below.

#include <stdlib.h>

#include <libint2.hpp>

using libint2::BasisSet;
using libint2::Operator;
using libint2::Shell;

int main(int argc, char* argv[])
{

  printf("libint2 version %s\n", LIBINT_VERSION);

  std::vector<libint2::Shell> shells;
  shells.push_back({
      {0.100000},
      {// S shell, spherical=false, contraction coefficients
       {0, false, {1.0}}},
      {{0, 0, 0}} // origin coordinates
  });
  shells.push_back({
      {0.100000},
      {// P shell, spherical=false, contraction coefficients
       {1, false, {1.0}}},
      {{0, 0, 0}} // origin coordinates
  });

#if LIBINT_MAJOR_VERSION == 2 && LIBINT_MINOR_VERSION >= 8
  auto obs = libint2::BasisSet(std::move(shells));
  auto shell2bf = obs.shell2bf();
#else
  libint2::BasisSet obs = libint2::BasisSet();
  for (libint2::Shell shell : shells) {
    obs.push_back(std::move(shell));
  }
  auto shell2bf = BasisSet::compute_shell2bf(obs);
#endif

  libint2::initialize();

  printf("2-center ERI derivatives\n\n\n");
  { // compute 2-center ERI
    libint2::Engine engine = libint2::Engine(libint2::Operator::coulomb, libint2::max_nprim(obs),
                                             libint2::max_l(obs), 1);
    engine.set(libint2::BraKet::xs_xs);

    const auto& buf = engine.results();

    for (auto s1 = 0; s1 != obs.size(); ++s1) {
      const auto bf1_first = shell2bf[s1]; // first basis function in this shell
      auto n1 = obs[s1].size();            // number of basis functions in shell 1
      for (auto s2 = 0; s2 != obs.size(); ++s2) {
        const auto bf2_first = shell2bf[s2]; // first basis function in this shell
        auto n2 = obs[s2].size();            // number of basis functions in shell 2
        engine.compute(obs[s1], obs[s2]);
        const auto* buf_1234 = buf[0];
        for (auto f1 = 0ul, f1234 = 0ul; f1 != n1; ++f1) {
          const auto bf1 = f1 + bf1_first;
          for (auto f2 = 0ul; f2 != n2; ++f2, ++f1234) {
            const auto bf2 = f2 + bf2_first;
            const auto integral = buf_1234[f1234];
            printf("d/dx (%lu, %lu), %lf \n", bf1, bf2, integral);
          }
        }
        printf("\n");
      }
    }
  }

  printf("3-center ERI derivatives\n\n\n");
  { // compute 3-center ERI
    libint2::Engine engine = libint2::Engine(libint2::Operator::coulomb, libint2::max_nprim(obs),
                                             libint2::max_l(obs), 1);
    engine.set(libint2::BraKet::xs_xx);

    const auto& buf = engine.results();

    for (auto s1 = 0; s1 != obs.size(); ++s1) {
      const auto bf1_first = shell2bf[s1]; // first basis function in this shell
      auto n1 = obs[s1].size();            // number of basis functions in shell 1
      for (auto s2 = 0; s2 != obs.size(); ++s2) {
        const auto bf2_first = shell2bf[s2]; // first basis function in this shell
        auto n2 = obs[s2].size();            // number of basis functions in shell 2
        for (auto s3 = 0; s3 != obs.size(); ++s3) {
          const auto bf3_first = shell2bf[s3]; // first basis function in this shell
          auto n3 = obs[s3].size();            // number of basis functions in shell 3
          engine.compute(obs[s1], obs[s2], obs[s3]);
          const auto* buf_1234 = buf[0];
          for (auto f1 = 0ul, f1234 = 0ul; f1 != n1; ++f1) {
            const auto bf1 = f1 + bf1_first;
            for (auto f2 = 0ul; f2 != n2; ++f2) {
              const auto bf2 = f2 + bf2_first;
              for (auto f3 = 0ul; f3 != n3; ++f3, ++f1234) {
                const auto bf3 = f3 + bf3_first;
                const auto integral = buf_1234[f1234];
                printf("d/dx (%lu, %lu, %lu), %lf \n", bf1, bf2, bf3, integral);
              }
            }
          }
          printf("\n");
        }
      }
    }
  }
  libint2::finalize();
  return 0;
}

The full output I see from this is:

libint2 version 2.8.1
2-center ERI derivatives

d/dx (0, 0), 0.000000

d/dx (0, 1), 13.246118
d/dx (0, 2), 0.000000
d/dx (0, 3), 0.000000

d/dx (1, 0), -13.246118
d/dx (2, 0), 0.000000
d/dx (3, 0), 0.000000

d/dx (1, 1), 0.000000
d/dx (1, 2), 0.000000
d/dx (1, 3), 0.000000
d/dx (2, 1), 0.000000
d/dx (2, 2), 0.000000
d/dx (2, 3), 0.000000
d/dx (3, 1), 0.000000
d/dx (3, 2), 0.000000
d/dx (3, 3), 0.000000

3-center ERI derivatives

d/dx (0, 0, 0), 0.000000

mintest: /usr/local/include/libint2/engine.impl.h:1875: const target_ptr_vec& libint2::Engine::compute2(const libint2::Shell&, const libint2::Shell&, const libint2::Shell&, const libint2::Shell&, const libint2::ShellPair*, const libint2::ShellPair*) [with libint2::Operator oper = libint2::Operator::coulomb; libint2::BraKet braket = libint2::BraKet::xs_xx; long unsigned int deriv_order = 1; libint2::Engine::target_ptr_vec = std::vector<const double*, libint2::detail::ext_stack_allocator<const double*, 25> >]: Assertion `buildfnptrs_[buildfnidx] && "null build function ptr"' failed.
Aborted

So could this be a bug? Or maybe it is time for me to learn the C interface...

update: I updated the code snippet and results to correspond to 2.8.1

[evaleev]

Shane, this is likely a logic error in the Engine ... I will try to reproduce. Unlike 1e 2c and 2e 4c derivatives which have been stressed hard 2e 2c and 2e 3c have not been tested.

[evaleev]

So could this be a bug? Or maybe it is time for me to learn the C interface...

NOOOOOOOOOOOOOOOOOOOOOOOO ...

[smparker]

I learned a couple more details about this that I wanted to share just in case they'd be helpful:

• I get basically the same error with 2.8.1 (though I did have to update the example above, which I'll edit shortly)
• It seems that (p|ss) integral derivatives run through, but it crashes at (s|sp) and (s|ps)

[smparker]

Okay, I pretty much know what's going on now. It seems to be caused by using different max angular momentum values for the eri2 and eri3 options. Originally, I used --with-max-am=4 --with-eri-max-am=4 --with-eri2-max-am=6 --with-eri3-max-am=6 when generating the library. When I replace all that with just --with-max-am=6 then the above code seems to work just fine. I can live with that as a solution.

[evaleev]

@smparker this is addressed by #324 , will try to push out 2.8.2 today/tomorrow