Fix indexing in diff Jastrow component

Second part for fix of QMCPACK#2814

The myVars in this object contains all the variational parameters
aggregated from all the radial functions (in their myVars), including
the non-active ones.

The OffSet variable maps the indices from the radial function myVars to
myVars in this object.

The output for parameter derivatives should use a contiguous index (i.e.
all the non-active parameters removed).  This is what the Index value in
VariableSet maps (returned by myVars.where).
In the code, index k covers all variables in myVars.  The index kk is
the contiguous index for output.

Some of the arrays are indexed with k, but should be indexed with kk
instead.

src/QMCWaveFunctions/Jastrow/DiffTwoBodyJastrowOrbital.h

@@ -187,7 +187,7 @@ class DiffTwoBodyJastrowOrbital : public DiffWaveFunctionComponent
         continue;
       if (active.recompute(kk))
         recalculate = true;
-      rcsingles[k] = true;
+      rcsingles[kk] = true;
     }
     if (recalculate)
     {
@@ -196,9 +196,9 @@ class DiffTwoBodyJastrowOrbital : public DiffWaveFunctionComponent
         int kk = myVars.where(k);
         if (kk < 0)
           continue;
-        if (rcsingles[k])
+        if (rcsingles[kk])
         {
-          dhpsioverpsi[kk] = -RealType(0.5) * ValueType(Sum(*lapLogPsi[k])) - ValueType(Dot(P.G, *gradLogPsi[k]));
+          dhpsioverpsi[kk] = -RealType(0.5) * ValueType(Sum(*lapLogPsi[kk])) - ValueType(Dot(P.G, *gradLogPsi[kk]));
         }
       }
     }
@@ -217,7 +217,7 @@ class DiffTwoBodyJastrowOrbital : public DiffWaveFunctionComponent
         continue;
       if (active.recompute(kk))
         recalculate = true;
-      rcsingles[k] = true;
+      rcsingles[kk] = true;
     }
     if (recalculate)
     {
@@ -285,9 +285,9 @@ class DiffTwoBodyJastrowOrbital : public DiffWaveFunctionComponent
         int kk = myVars.where(k);
         if (kk < 0)
           continue;
-        if (rcsingles[k])
+        if (rcsingles[kk])
         {
-          dlogpsi[kk] = dLogPsi[k];
+          dlogpsi[kk] = dLogPsi[kk];
         }
         //optVars.setDeriv(p,dLogPsi[ip],-0.5*Sum(*lapLogPsi[ip])-Dot(P.G,*gradLogPsi[ip]));
       }

src/QMCWaveFunctions/tests/test_DiffTwoBodyJastrowOrbital.cpp

@@ -9,14 +9,12 @@
 // File created by: Mark Dewing, mdewing@anl.gov, Argonne National Laboratory
 //////////////////////////////////////////////////////////////////////////////////////
 
-
 #include "catch.hpp"
 
 #include "QMCWaveFunctions/Jastrow/DiffTwoBodyJastrowOrbital.h"
 
 namespace qmcplusplus
 {
-
 // Mock radial function to use in DiffTwoBodyJastrowOrbital
 class FakeJastrow
 {
@@ -27,9 +25,15 @@ class FakeJastrow
 
   using RealType = QMCTraits::RealType;
 
-  bool evaluateDerivatives(RealType r, std::vector<TinyVector<RealType, 3>>& derivs) { return true; }
+  bool evaluateDerivatives(RealType r, std::vector<TinyVector<RealType, 3>>& derivs)
+  {
+    derivs = derivs_;
+    return true;
+  }
 
   void resetParameters(const opt_variables_type& active) {}
+
+  std::vector<TinyVector<RealType, 3>> derivs_;
 };
 
 TEST_CASE("DiffTwoBodyJastrowOrbital simple", "[wavefunction]")
@@ -69,10 +73,20 @@ ParticleSet get_two_species_particleset()
   std::vector<int> ud{2, 2};
   elec.setName("e");
   elec.create(ud);
+
+  elec.R[0] = {1.0, 0.0, 0.0};
+  elec.R[1] = {1.1, 1.0, 0.1};
+  elec.R[2] = {0.9, 0.8, 1.0};
+  elec.R[3] = {0.9, 0.5, 1.1};
+
   SpeciesSet& tspecies = elec.getSpeciesSet();
   int upIdx            = tspecies.addSpecies("u");
   int downIdx          = tspecies.addSpecies("d");
   elec.resetGroups();
+
+  elec.addTable(elec);
+  elec.update();
+
   return elec;
 }
 
@@ -184,4 +198,71 @@ TEST_CASE("DiffTwoBodyJastrowOrbital variables fail", "[wavefunction]")
   CHECK(o4.first == -1);
 }
 
+TEST_CASE("DiffTwoBodyJastrowOrbital three variables", "[wavefunction]")
+{
+  ParticleSet elec = get_two_species_particleset();
+
+  DiffTwoBodyJastrowOrbital<FakeJastrow> jorb(elec);
+
+  FakeJastrow j2a;
+  j2a.myVars.insert("opt1", 1.0);
+  j2a.name = "j2a";
+  // update num_active_vars
+  j2a.myVars.resetIndex();
+  jorb.addFunc(0, 0, &j2a);
+
+  FakeJastrow j2b;
+  j2b.myVars.insert("opt2", 2.0);
+  j2b.name = "j2b";
+  // update num_active_vars
+  j2b.myVars.resetIndex();
+  jorb.addFunc(0, 1, &j2b);
+
+  opt_variables_type global_active;
+  global_active.insertFrom(j2b.myVars);
+  global_active.resetIndex();
+
+
+  jorb.checkOutVariables(global_active);
+
+  //global_active.print(std::cout,0,true);
+  //jorb.getVars().print(std::cout,0,true);
+
+  CHECK(global_active.size_of_active() == 1);
+  // Not optimizing the parameter in this Jastrow factor, indicated by first index is -1
+  auto o1 = jorb.getOffset(0);
+  CHECK(o1.first == -1);
+
+  // Offset into set of active variables (global_active)
+  auto o2 = jorb.getOffset(1);
+  CHECK(o2.first == 0);
+  CHECK(o2.second == 1);
+
+  auto o3 = jorb.getOffset(2);
+  CHECK(o3.first == 0);
+  CHECK(o3.second == 1);
+
+  auto o4 = jorb.getOffset(3);
+  CHECK(o4.first == -1);
+
+
+  using ValueType = QMCTraits::ValueType;
+
+  // Check derivative indexing
+  int num_vars = 1;
+  j2b.derivs_.resize(num_vars);
+  // Elements are d/dp_i u(r), d/dp_i du/dr,  d/dp_i d2u/dr2
+  j2b.derivs_[0] = {0.5, 1.3, 2.4};
+  std::vector<ValueType> dlogpsi(num_vars);
+  jorb.evaluateDerivativesWF(elec, global_active, dlogpsi);
+
+  CHECK(dlogpsi[0] == ValueApprox(-2.0)); // 4 * derivs_[0][0]
+
+  std::vector<ValueType> dlogpsi2(num_vars);
+  std::vector<ValueType> dhpsioverpsi(num_vars);
+
+  jorb.evaluateDerivatives(elec, global_active, dlogpsi2, dhpsioverpsi);
+  CHECK(dlogpsi2[0] == ValueApprox(-2.0));
+}
+
 } // namespace qmcplusplus