Integrals of BF with several primitives seem to evaluate only one

[HalblauterPC]

When using the library as follows:

[...]
 libint2::BasisSet ao_basis = BasisSet(basisname, atoms);
 using libint2::Engine;

 Engine overlap_dipole_engine(Operator::emultipole1, // will compute overlap + electric dipole moments
                               ao_basis.max_nprim(), ao_basis.max_l());

 for (int i = 0; i < Nthreads; i++)
 {
    overlap_dipole_engines[i] = overlap_dipole_engine;
 }

 [...]

 // with a thread id
{
 const auto& buf_vec1 = overlap_dipole_engines[id].results(); // will point to computed shell sets
                                                  // const auto& is very important!

  auto shell2bf = ao_basis.shell2bf();

   // iterate over all shells
   for (int s1 = shells_begin[id]; s1 < shells_end[id]; ++s1) {
    auto bf1 = shell2bf[s1];
      for (auto s2 = 0; s2 < Nshells; ++s2) {
          auto bf2 = shell2bf[s2];

          overlap_dipole_engines[id].compute(ao_basis[s1], ao_basis[s2]);

          auto n1 = ao_basis[s1].size(); // number of basis functions in first shell
          auto n2 = ao_basis[s2].size(); // number of basis functions in second shell

          if (n1 == 1 && n2 == 1) {
              int orbital_index = bf1 * Norbitals + bf2;
              smat[orbital_index] = buf_vec1[0][0];
               Dx[orbital_index] = -buf_vec1[0][1];
               Dy[orbital_index] = -buf_vec1[0][2];
               Dz[orbital_index] = -buf_vec1[0][3];
          } 
          else {
              auto s_shellset = buf_vec1[0]; // skipped nullptr check for simplicity
              auto mu_x_shellset = buf_vec1[1];
              auto mu_y_shellset = buf_vec1[2];
              auto mu_z_shellset = buf_vec1[3];

              // iterate over all basis functions in the shell
              for (auto f1 = 0; f1 < n1; ++f1) {
                  for (auto f2 = 0; f2 < n2; ++f2) {
                      int shell_index = f1 * n2 + f2;
                      int orbital_index = (bf1+f1) * Norbitals + bf2+f2;

                      smat[orbital_index] = s_shellset[shell_index];
                      Dx[orbital_index] = -mu_x_shellset[shell_index];
                      Dy[orbital_index] = -mu_y_shellset[shell_index];
                      Dz[orbital_index] = -mu_z_shellset[shell_index];
                  }
              }
          }
      }
  }
}

After printing out the results (Here the diagonals for H2 in the aug-cc-pVDZ basis), I get:

Shell > Bf (NPrimitives) | Overlap | Dx | Dy | Dz
0 > 0 (3) | 0.345314 | -0 | -0 | 0.24172
1 > 1 (1) | 1 | -0 | -0 | 0.7
2 > 2 (1) | 1 | -0 | -0 | 0.7
2 > 3 (1) | 1 | -0 | -0 | 0.7
2 > 4 (1) | 1 | -0 | -0 | 0.7
3 > 5 (1) | 1 | -0 | -0 | 0.7
4 > 6 (1) | 1 | -0 | -0 | 0.7
4 > 7 (1) | 1 | -0 | -0 | 0.7
4 > 8 (1) | 1 | -0 | -0 | 0.7
5 > 9 (3) | 0.345314 | -0 | -0 | -0.24172
6 > 10 (1) | 1 | -0 | -0 | -0.7
7 > 11 (1) | 1 | -0 | -0 | -0.7
7 > 12 (1) | 1 | -0 | -0 | -0.7
7 > 13 (1) | 1 | -0 | -0 | -0.7
8 > 14 (1) | 1 | -0 | -0 | -0.7
9 > 15 (1) | 1 | -0 | -0 | -0.7
9 > 16 (1) | 1 | -0 | -0 | -0.7
9 > 17 (1) | 1 | -0 | -0 | -0.7

where the BF with several primitives seem only to evaluate one of them. Am I missing something?

[HalblauterPC]

Update:

I get the same integral values, when I sum up the integrals by primitives individually as follows

for (int p1 = 0; p1 < m1; p1++){ 
                for (int p2 = 0; p2 < m2; p2++){

                  auto prim1 = ao_basis[s1].extract_primitive(p1); //normalized primitive
                  auto prim2 = ao_basis[s2].extract_primitive(p2); //normalized primitive

                  overlap_dipole_engines[id].compute(prim1, prim2);
                  kinetic_engines[id].compute(prim1, prim2);
                  potential_engines[id].compute(prim1, prim2);

                  if (n1 == 1 && n2 == 1) {
                      int orbital_index = bf1 * Norbitals + bf2;
                      auto coef = (ao_basis[s1].coeff_normalized(0,p1)*ao_basis[s2].coeff_normalized(0,p2));

                      smat[orbital_index] += coef*buf_vec1[0][0]; // evtl. buf_vec1[0][0]
                      Dx[orbital_index] -= coef*buf_vec1[0][1];
                      Dy[orbital_index] -= coef*buf_vec1[0][2];
                      Dz[orbital_index] -= coef*buf_vec1[0][3];
                  } 
}

I also tried with the direct coefficients (auto c1 = ao_basis[s1].contr.at(0).coeff[p1];) and changing to extract_primitive(p1, false), but I still don't get normalized overlap in any combination of those.

[HalblauterPC]

Also, when printing the basis set, it differs quite from the coefficients saved in lib/basis/ for several basis sets tested. e.g.

Shell:( O={0,0,-0.699999909988376}
   {l=0,sph=0}
  13.01 0.0961066598498265
  1.962 0.163020057680412
  0.4446 0.185545117025848
1 3

Shell:( O={0,0,-0.699999909988376}
   {l=0,sph=0}
  0.122 0.147122794428574
1 1

Shell:( O={0,0,-0.699999909988376}
   {l=1,sph=0}
  0.727 0.956881375059565
3 1
...

vs.

H     0 
S   3   1.00
     13.0100000              0.0196850        
      1.9620000              0.1379770        
      0.4446000              0.4781480        
S   1   1.00
      0.1220000              1.0000000        
P   1   1.00
      0.7270000              1.0000000 

[HalblauterPC]

The problem was that I included the libint2::Shell::do_enforce_unit_normalization(false); line of the Hartree-Fock example. This disables the normalization of convoluted Shells. I don't know why this is disabled in the HF example, but for normalized overlap, which I need here, this is not correct. Removing the line or changing it to true corrects the behavior.