Fix EnergyAndForces tests

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())

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));
     }
 

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()"

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