Two body Jastrow issue when there are more than 2 particle types but only one particle of each type

[shivupa]

Describe the bug
Thanks to @markdewing, I was able to progress further until I broke things again 😄

When running a calculation with 1 particle in each of >2 particle types, there is a segfault.

To Reproduce
Steps to reproduce the behavior:

1. b1a4692-dirty
2. cmake command : cmake -DCMAKE_BUILD_TYPE=Debug ..
3. qmcpack H_wpos.qmc.in-wfj.xml
   The files needed to replicate this issue are located here https://gist.github.com/shivupa/e390aa40b77d6957d72f59b3c1be7203.

Expected behavior
The calculation wouldn't segfault. I think the two body jastrow builder should bypass building the uu jastrow?

System:

• Workstation Intel i7-2600S

Additional context
Potential Solution?:
src/QMCWaveFunctions/Jastrow/DiffTwoBodyJastrowOrbital.h 94 seems to guard against the 1 elec up/1 elec down case. This might be the source of the issue.

Portion of Output before Seg Fault:

========================================================
--- Memory usage report : Memory Usage after cloning ---
========================================================
Available memory on node 0, free + buffers :    6696 MiB
Memory footprint by rank 0 on node 0       :      25 MiB
========================================================
  Initial partition of walkers 0 1 2 3 4 5 6 7 8 

  Using Particle by Particle moves
  Walker moves without drift
  Total Sample Size   =0
  Walker distribution on root = 0 1 2 3 4 5 6 7 8 
[1]    8353 segmentation fault (core dumped)  ~/Documents/qmcpack/build/bin/qmcpack H_wpos.qmc.in-wfj.xml

gdb trouble shooting

Thread 1 "qmcpack" received signal SIGSEGV, Segmentation fault.
0x00005555558ee984 in qmcplusplus::BsplineFunctor<double>::evaluateVGL (this=0x0, iat=0, iStart=0, 
    iEnd=0, _distArray=0x555556157e20, _valArray=0x55555612bc40, _gradArray=0x55555612bcc0, 
    _laplArray=0x55555612bd40, distArrayCompressed=0x55555602dfe0, distIndices=0x555556128720)
    at /home/shiv/Documents/qmcpack/src/QMCWaveFunctions/Jastrow/BsplineFunctor.h:687
687	  real_type dSquareDeltaRinv = DeltaRInv * DeltaRInv;
(gdb) bt
#0  0x00005555558ee984 in qmcplusplus::BsplineFunctor<double>::evaluateVGL (this=0x0, iat=0, 
    iStart=0, iEnd=0, _distArray=0x555556157e20, _valArray=0x55555612bc40, _gradArray=0x55555612bcc0, 
    _laplArray=0x55555612bd40, distArrayCompressed=0x55555602dfe0, distIndices=0x555556128720)
    at /home/shiv/Documents/qmcpack/src/QMCWaveFunctions/Jastrow/BsplineFunctor.h:687
#1  0x00005555558e5a3d in qmcplusplus::J2OrbitalSoA<qmcplusplus::BsplineFunctor<double> >::computeU3 (
    this=0x55555612b660, P=..., iat=0, dist=warning: RTTI symbol not found for class 'qmcplusplus::Vector<double, qmcplusplus::Mallocator<double, 32ul> >'
..., u=0x55555612bc40, du=0x55555612bcc0, 
    d2u=0x55555612bd40, triangle=true)
    at /home/shiv/Documents/qmcpack/src/QMCWaveFunctions/Jastrow/J2OrbitalSoA.h:485
#2  0x00005555558cd900 in qmcplusplus::J2OrbitalSoA<qmcplusplus::BsplineFunctor<double> >::recompute (
    this=0x55555612b660, P=...)
    at /home/shiv/Documents/qmcpack/src/QMCWaveFunctions/Jastrow/J2OrbitalSoA.h:605
#3  0x00005555558cee34 in qmcplusplus::J2OrbitalSoA<qmcplusplus::BsplineFunctor<double> >::evaluateGL
    (this=0x55555612b660, P=..., G=warning: RTTI symbol not found for class 'qmcplusplus::ParticleAttrib<qmcplusplus::TinyVector<double, 3u>, std::allocator<qmcplusplus::TinyVector<double, 3u> > >'
..., L=..., fromscratch=true)
    at /home/shiv/Documents/qmcpack/src/QMCWaveFunctions/Jastrow/J2OrbitalSoA.h:662
#4  0x00005555558cd798 in qmcplusplus::J2OrbitalSoA<qmcplusplus::BsplineFunctor<double> >::evaluateLog
    (this=0x55555612b660, P=..., G=warning: RTTI symbol not found for class 'qmcplusplus::ParticleAttrib<qmcplusplus::TinyVector<double, 3u>, std::allocator<qmcplusplus::TinyVector<double, 3u> > >'
..., L=...)
    at /home/shiv/Documents/qmcpack/src/QMCWaveFunctions/Jastrow/J2OrbitalSoA.h:652
#5  0x0000555555861c01 in qmcplusplus::TrialWaveFunction::evaluateLog (this=0x5555560274e0, P=...)
    at /home/shiv/Documents/qmcpack/src/QMCWaveFunctions/TrialWaveFunction.cpp:119

I didn't find this incredibly enlightening besides this resulting in a divide by zero.

Let me know what you think.

[markdewing]

Changing J2OrbitalSoA.h, line 409 from

    if (N == 2)

to

    if (N == NumGroups)

fixes the issue for me. I don't know if the same change needs to be done in DiffTwoBodyJastrowOrbital.h or not.

[markdewing]

The problem comes in from J2OrbitalSoA.h::computeU3 where it loops over all the functors for all inter-particle types, including the same particle interactions. These don't exist if there is only one particle of each type, so the functor (f2) is NULL. Access to the member variable DeltaRInv fails. The solutions are either to check for a NULL functor, or ensure all the functors are defined.

The addFunc function attempts to make sure all of the functors are defined. I think the fix I suggested will not work because it sets all the functors to the same value. What might need to be done is a similar check (or modify the existing check) to just set the same-particle interactions.

It looks like the current special case (N == 2) ensures the same-particle functors are set, and sets uu and dd to the same functor. For the fix, we only want the first part of this behavior.

[markdewing]

Here's a fix that should work. Two changes

1. The special case fix only sets the diagonal (same particle) interactions
2. The general case is pulled out of the else so it always runs for dissimilar particles.

I don't know if this will fix the original N==2 situation - that will need to be checked.

    if (N == NumGroups)
    {
      // a very special case, 1 up + 1 down
      // uu/dd was prevented by the builder
      for (int ig = 0; ig < NumGroups; ++ig)
          F[ig * NumGroups + ig] = j;
    }

    // generic case
    F[ia * NumGroups + ib] = j;
    F[ib * NumGroups + ia] = j;

[shivupa]

Ah Thanks for the quick replies. I will give this a try now

[shivupa]

This does fix the issue and it doesn't seem to be breaking the u/d only case. (All the deterministic tests pass and I think the u/d case is covered by https://github.com/QMCPACK/qmcpack/blob/develop/src/QMCWaveFunctions/tests/test_wf.cpp since it was also modified in #364)