[WIP] Energy Density Bug Fixes

src/Containers/OhmmsPETE/tests/CMakeLists.txt

@@ -13,6 +13,6 @@ set(UTEST_EXE test_containers_ohmmspete)
 set(UTEST_NAME deterministic-unit_${UTEST_EXE})
 
 add_executable(${UTEST_EXE} test_Vector.cpp test_Matrix.cpp test_Array.cpp test_TinyVector.cpp)
-target_link_libraries(${UTEST_EXE} catch_main containers)
+target_link_libraries(${UTEST_EXE} catch_main containers utilities_for_test)
 
 add_unit_test(${UTEST_NAME} 1 1 $<TARGET_FILE:${UTEST_EXE}>)

src/Containers/OhmmsPETE/tests/test_Matrix.cpp

@@ -17,6 +17,7 @@
 
 #include "OhmmsPETE/OhmmsMatrix.h"
 #include "config.h"
+#include <checkMatrix.hpp>
 
 using std::string;
 
@@ -53,14 +54,28 @@ TEST_CASE("matrix", "[OhmmsPETE]")
     CHECK(*ia == Approx(3.1));
   }
 
-  REQUIRE( A == A);
-  REQUIRE( A != C);
+  REQUIRE(A == A);
+  REQUIRE(A != C);
 
   // copy constructor and copy method
   Mat D(C);
   CHECK(D.rows() == 3);
   CHECK(D.cols() == 3);
 
+  // Check that the zeroing behavior on construction and resize
+  // follows assumptions about zeroing.
+  Mat E(2, 2);
+  Mat E_expected(2, 2);
+  E_expected        = 0.0;
+  auto check_matrix = checkMatrix(E, E_expected);
+  CHECKED_ELSE(check_matrix.result) { FAIL(check_matrix.result_message); }
+
+  E.resize(5, 5);
+  Mat E2_expected(5, 5);
+  E2_expected  = 0.0;
+  check_matrix = checkMatrix(E, E2_expected);
+  CHECKED_ELSE(check_matrix.result) { FAIL(check_matrix.result_message); }
+
   // swap_rows
   A(0, 0) = 0.0;
   A(0, 1) = 1.0;
@@ -81,49 +96,49 @@ TEST_CASE("matrix", "[OhmmsPETE]")
 
 TEST_CASE("matrix converting assignment", "[OhmmsPETE]")
 {
-  Matrix<double> mat_A(3,3);
-  Matrix<float> mat_B(3,3);
+  Matrix<double> mat_A(3, 3);
+  Matrix<float> mat_B(3, 3);
 
   for (int i = 0; i < 3; i++)
     for (int j = 0; j < 3; j++)
       mat_A(i, j) = (i + j) * 2.1;
 
   mat_B = mat_A;
 
-  CHECK(mat_B(0,0) == Approx(0));
-  CHECK(mat_B(1,1) == Approx(4.2));
+  CHECK(mat_B(0, 0) == Approx(0));
+  CHECK(mat_B(1, 1) == Approx(4.2));
 
-  Matrix<float> mat_C(2,2);
+  Matrix<float> mat_C(2, 2);
   for (int i = 0; i < 2; i++)
     for (int j = 0; j < 2; j++)
       mat_C(i, j) = (i + j) * 2.2;
 
   mat_A.assignUpperLeft(mat_C);
-  CHECK(mat_A(1,0) == Approx(2.2));
-  CHECK(mat_A(1,2) == Approx(6.3));
+  CHECK(mat_A(1, 0) == Approx(2.2));
+  CHECK(mat_A(1, 2) == Approx(6.3));
 
-  Matrix<float> mat_D(4,4);
+  Matrix<float> mat_D(4, 4);
   for (int i = 0; i < 4; i++)
     for (int j = 0; j < 4; j++)
       mat_D(i, j) = (i + j) * 2.3;
 
   mat_A.assignUpperLeft(mat_D);
-  CHECK(mat_A(1,0) == Approx(2.3));
-  CHECK(mat_A(1,2) == Approx(6.9));
+  CHECK(mat_A(1, 0) == Approx(2.3));
+  CHECK(mat_A(1, 2) == Approx(6.9));
 
-  Matrix<float> mat_too_small(2,2);
+  Matrix<float> mat_too_small(2, 2);
   CHECK_THROWS(mat_too_small = mat_A);
-  
-  Matrix<double> mat_too_small_but_larger_value_type(2,2);
+
+  Matrix<double> mat_too_small_but_larger_value_type(2, 2);
   mat_too_small_but_larger_value_type = mat_B;
   CHECK(mat_too_small_but_larger_value_type.rows() == 3);
   CHECK(mat_too_small_but_larger_value_type.cols() == 3);
-  CHECK(mat_too_small_but_larger_value_type(1,1) == Approx(mat_B(1,1)));
-  Matrix<double> too_small_but_same(2,2);
+  CHECK(mat_too_small_but_larger_value_type(1, 1) == Approx(mat_B(1, 1)));
+  Matrix<double> too_small_but_same(2, 2);
   too_small_but_same = mat_A;
   CHECK(too_small_but_same.rows() == 3);
   CHECK(too_small_but_same.cols() == 3);
-  CHECK(too_small_but_same(1,1) == Approx(mat_A(1,1)));  
+  CHECK(too_small_but_same(1, 1) == Approx(mat_A(1, 1)));
 }
 
 } // namespace qmcplusplus

src/Estimators/EnergyDensityEstimator.cpp

@@ -148,8 +148,11 @@ void NEEnergyDensityEstimator::registerListeners(QMCHamiltonian& ham_leader)
   QMCHamiltonian::mw_registerKineticListener(ham_leader, kinetic_listener);
   ListenerVector<Real> potential_listener("potential", getListener(local_pot_values_));
   QMCHamiltonian::mw_registerLocalPotentialListener(ham_leader, potential_listener);
-  ListenerVector<Real> ion_potential_listener("potential", getListener(local_ion_pot_values_));
-  QMCHamiltonian::mw_registerLocalIonPotentialListener(ham_leader, ion_potential_listener);
+  if (pset_static_)
+  {
+    ListenerVector<Real> ion_potential_listener("potential", getListener(local_ion_pot_values_));
+    QMCHamiltonian::mw_registerLocalIonPotentialListener(ham_leader, ion_potential_listener);
+  }
 }
 
 /** This function collects the per particle energies.

src/Estimators/EnergyDensityInput.h

@@ -75,7 +75,9 @@ class EnergyDensityInput
   std::string name_{type_tag};
   std::string type_{type_tag};
   std::string dynamic_{"e"};
-  std::string static_{"ion"};
+  std::string static_{}; // There can't be a default for this since
+                         // some systems, heg namely don't
+                         // have static psets
   bool ion_points_{false};
   EnergyDensityInputSection input_section_;
   ReferencePointsInput ref_points_input_;

src/Estimators/EstimatorManagerCrowd.h

@@ -32,22 +32,22 @@ class QMCHamiltonian;
  *
  *  Stepping away from the CloneManger + clones design which creates EstimatorManagers
  *  Which operate differently based on internal switches.
- *  
+ *
  *  see EstimatorManagerNew.h for full description of the new design.
  */
 class EstimatorManagerCrowd
 {
 public:
-  using MCPWalker = Walker<QMCTraits, PtclOnLatticeTraits>;
-  using RealType  = EstimatorManagerNew::RealType;
+  using MCPWalker        = Walker<QMCTraits, PtclOnLatticeTraits>;
+  using RealType         = EstimatorManagerNew::RealType;
   using FullPrecRealType = EstimatorManagerNew::FullPrecRealType;
 
   /** EstimatorManagerCrowd are always spawn of an EstimatorManagerNew
    */
   EstimatorManagerCrowd(EstimatorManagerNew& em);
 
   ///destructor
-  ~EstimatorManagerCrowd(){};
+  ~EstimatorManagerCrowd() {};
 
   ///return the number of ScalarEstimators
   inline int size() const { return scalar_estimators_.size(); }

src/Estimators/MagnetizationDensity.h

@@ -30,13 +30,13 @@ class MagnetizationDensityTests;
 }
 /** Magnetization density estimator for non-collinear spin calculations.
  *
- * As documented in the manual, the following formula is computed:  
+ * As documented in the manual, the following formula is computed:
  *\mathbf{m}_c = \int d\mathbf{X} \left|{\Psi(\mathbf{X})}\right|^2 \int_{\Omega_c}d\mathbf{r} \sum_i\delta(\mathbf{r}-\hat{\mathbf{r}}_i)\int_0^{2\pi} \frac{ds'_i}{2\pi} \frac{\Psi(\ldots \mathbf{r}_i s'_i \ldots )}{\Psi(\ldots \mathbf{r}_i s_i \ldots)}\langle s_i | \hat{\sigma} | s'_i \rangle
  *
  * The way that the above relates to the underlying data structure data_ is as follows.  We grid space up and assign an index
- * for each of the real space bins (identical to SpinDensityNew).  To account for the fact that magnetization is vectorial, we 
+ * for each of the real space bins (identical to SpinDensityNew).  To account for the fact that magnetization is vectorial, we
  * triple the length of this array.  If grid_i is the index of the real space gridpoint i, then the data is layed out like:
- * [grid_0_x, grid_0_y, grid_0_z, grid_1_x, ..., grid_N_x, grid_N_y, grid_N_z].  This is also the way it is stored in HDF5.  
+ * [grid_0_x, grid_0_y, grid_0_z, grid_1_x, ..., grid_N_x, grid_N_y, grid_N_z].  This is also the way it is stored in HDF5.
  *
  */
 class MagnetizationDensity : public OperatorEstBase
@@ -93,9 +93,9 @@ class MagnetizationDensity : public OperatorEstBase
   /**
   * Generates the spin integrand \Psi(s')/Psi(s)* \langle s | \vec{\sigma} | s'\rangle for a specific
   *  electron iat.  Since this is a vectorial quantity, this function returns sx, sy, and sz in their own
-  *  arrays.  
+  *  arrays.
   *
-  * @param[in] pset ParticleSet  
+  * @param[in] pset ParticleSet
   * @param[in] wfn TrialWaveFunction
   * @param[in] iat electron index
   * @param[out] sx x component of spin integrand
@@ -110,51 +110,51 @@ class MagnetizationDensity : public OperatorEstBase
                              std::vector<Value>& sz);
 
   /**
-  *  Implementation of Simpson's 1/3 rule to integrate a function on a uniform grid. 
+  *  Implementation of Simpson's 1/3 rule to integrate a function on a uniform grid.
   *
-  * @param[in] fgrid f(x), the function to integrate. 
+  * @param[in] fgrid f(x), the function to integrate.
   * @param[in] gridDx, the grid spacing for the uniform grid.  Assumed to be consistent with size of fgrid.
   * @return Value of integral.
   */
   Value integrateBySimpsonsRule(const std::vector<Value>& fgrid, Real gridDx) const;
 
 
   /**
-  * Convenience function to generate a grid between 0 and 2pi, consistent with nsamples_ and integration method.  
-  * Can be uniform or random, depending on choice of integrator.   
+  * Convenience function to generate a grid between 0 and 2pi, consistent with nsamples_ and integration method.
+  * Can be uniform or random, depending on choice of integrator.
   *
   * @param[out] sgrid A grid with nsamples_ points between 0 and 2pi.  Overwritten.
   */
   void generateGrid(std::vector<Real>& sgrid) const;
 
   /**
-  * Generate a uniform grid between [start,stop] for numerical quadrature.  
+  * Generate a uniform grid between [start,stop] for numerical quadrature.
   *
   * @param[out] sgrid Random number grid between "start" and "stop".  Number of points taken from size of sgrid.
-  * @param[in] start start of grid interval 
+  * @param[in] start start of grid interval
   * @param[in] stop end of grid interval
   */
   void generateUniformGrid(std::vector<Real>& sgrid, const Real start, const Real stop) const;
 
   /**
-  * Generate random grid between [start,stop] for MC integration.   
+  * Generate random grid between [start,stop] for MC integration.
   *
   * @tparam RAN_GEN Random number generator type.
   * @param[out] sgrid Random number grid between "start" and "stop".  Number of points taken from size of sgrid.
-  * @param[in] rng Random number generator queried to generate random points.   
-  * @param[in] start start of grid interval 
+  * @param[in] rng Random number generator queried to generate random points.
+  * @param[in] start start of grid interval
   * @param[in] stop end of grid interval
   */
 
   void generateRandomGrid(std::vector<Real>& sgrid, RandomBase<FullPrecReal>& rng, Real start, Real stop) const;
 
   /**
   * For a given spatial position r and spin component s, this returns the bin for accumulating the observable.
-  * 
+  *
   *
   * @param[in]  Position in real space. 3D
   * @param[in]  component, the x=0,y=1,z=2 component of the spin.
-  * @return Index of appropriate bin for this position and spin component 
+  * @return Index of appropriate bin for this position and spin component
   */
   size_t computeBin(const Position& r, const unsigned int component) const;
   MagnetizationDensityInput input_;

src/Estimators/NESpaceGrid.cpp

@@ -731,7 +731,18 @@ void NESpaceGrid<REAL>::collect(NESpaceGrid& reduction_grid, const RefVector<NES
 template<typename REAL>
 void NESpaceGrid<REAL>::normalize(REAL invToWgt)
 {
-  std::transform(data_.begin(), data_.end(), data_.begin(),
+  // I think its better to not transform the weights folling the
+  // convention that value 0 is always the weight.
+
+  // for (int i = 0, n = buffer_offset_; i < ndomains_; i++, n += nvalues_per_domain_)
+  // {
+  //   for (int v = 0; v < nvalues_per_domain_; v++)
+  //   {
+  //     data_[n + v] *= invToWgt;
+  //   }
+  // }
+
+  std::transform(data_.begin() + buffer_offset_, data_.end(), data_.begin() + buffer_offset_,
                  std::bind(std::multiplies<REAL>(), std::placeholders::_1, invToWgt));
 }