From 9ad93435e0421cc7cfd1cbb22d1df629bcf42135 Mon Sep 17 00:00:00 2001
From: "Fattebert J.-L" <fattebertj@ornl.gov>
Date: Wed, 25 Sep 2024 15:08:22 -0400
Subject: [PATCH] Fix EnergyAndForces tests

* have them work in debug mode too
---
 src/DensityMatrix.cc              |  7 ++++++-
 src/DensityMatrix.h               | 13 +++++++++++--
 src/ProjectedMatrices.cc          |  1 +
 tests/WFEnergyAndForces/mgmol.cfg |  2 +-
 4 files changed, 19 insertions(+), 4 deletions(-)

diff --git a/src/DensityMatrix.cc b/src/DensityMatrix.cc
index 6e6ea652..e2c7b473 100644
--- a/src/DensityMatrix.cc
+++ b/src/DensityMatrix.cc
@@ -28,7 +28,7 @@ const double factor_kernel4dot = 10.;
 template <class MatrixType>
 DensityMatrix<MatrixType>::DensityMatrix(const int ndim)
 {
-    assert(ndim >= 0);
+    assert(ndim > 0);
 
     dim_ = ndim;
 
@@ -45,6 +45,7 @@ DensityMatrix<MatrixType>::DensityMatrix(const int ndim)
     kernel4dot_ = new MatrixType("K4dot", ndim, ndim);
     work_       = new MatrixType("work", ndim, ndim);
     occupation_.resize(dim_);
+    setDummyOcc();
 }
 
 template <class MatrixType>
@@ -109,6 +110,7 @@ void DensityMatrix<MatrixType>::build(
     const std::vector<double>& occ, const int new_orbitals_index)
 {
     assert(dm_ != nullptr);
+    assert(!occ.empty());
 
     setOccupations(occ);
 
@@ -149,6 +151,8 @@ template <class MatrixType>
 void DensityMatrix<MatrixType>::setUniform(
     const double nel, const int new_orbitals_index)
 {
+    assert(!occupation_.empty());
+
     MGmol_MPI& mmpi  = *(MGmol_MPI::instance());
     const double occ = (double)((double)nel / (double)dim_);
     if (mmpi.instancePE0())
@@ -314,6 +318,7 @@ void DensityMatrix<MatrixType>::computeOccupations(const MatrixType& ls)
 template <class MatrixType>
 void DensityMatrix<MatrixType>::setOccupations(const std::vector<double>& occ)
 {
+    assert(!occ.empty());
 #ifdef PRINT_OPERATIONS
     MGmol_MPI& mmpi = *(MGmol_MPI::instance());
     if (mmpi.instancePE0())
diff --git a/src/DensityMatrix.h b/src/DensityMatrix.h
index 960ad18e..84804a86 100644
--- a/src/DensityMatrix.h
+++ b/src/DensityMatrix.h
@@ -84,17 +84,24 @@ class DensityMatrix
         *dm_            = mat;
         orbitals_index_ = orbitals_index;
 
-        occupation_.clear();
+        setDummyOcc();
 
         occ_uptodate_ = false;
         uniform_occ_  = false;
         stripped_     = false;
     }
 
+    // set occupations to meaningless values to catch uninitialized use
+    void setDummyOcc()
+    {
+        for (auto& occ : occupation_)
+            occ = -1.;
+    }
+
     void initMatrix(const double* const val)
     {
         dm_->init(val, dim_);
-        occupation_.clear();
+        setDummyOcc();
 
         occ_uptodate_ = false;
         uniform_occ_  = false;
@@ -105,8 +112,10 @@ class DensityMatrix
 
     void getOccupations(std::vector<double>& occ) const
     {
+        assert(!occupation_.empty());
         assert(occ_uptodate_);
         assert((int)occ.size() == dim_);
+
         memcpy(&occ[0], &occupation_[0], dim_ * sizeof(double));
     }
 
diff --git a/src/ProjectedMatrices.cc b/src/ProjectedMatrices.cc
index e0e679e9..3606aa40 100644
--- a/src/ProjectedMatrices.cc
+++ b/src/ProjectedMatrices.cc
@@ -509,6 +509,7 @@ template <class MatrixType>
 void ProjectedMatrices<MatrixType>::setOccupations(
     const std::vector<DISTMATDTYPE>& occ)
 {
+    assert(!occ.empty());
 #ifdef PRINT_OPERATIONS
     if (mmpi.instancePE0())
         (*MPIdata::sout) << "ProjectedMatrices<MatrixType>::setOccupations()"
diff --git a/tests/WFEnergyAndForces/mgmol.cfg b/tests/WFEnergyAndForces/mgmol.cfg
index d543b626..12002703 100644
--- a/tests/WFEnergyAndForces/mgmol.cfg
+++ b/tests/WFEnergyAndForces/mgmol.cfg
@@ -18,7 +18,7 @@ pseudopotential=pseudo.H
 type=QUENCH
 [Quench]
 max_steps=50
-atol=1.e-9
+atol=1.e-8
 num_lin_iterations=2
 [Orbitals]
 initial_type=Gaussian