Merge release (#293)

* Fix EnergyAndForces tests (#277)

* have them work in debug mode too

* Move factor 4pi out og linear solvers (#278)

* Move some code into PoissonSolverFactory (#279)

* Clean up class Potentials (#280)

* Clean up class Ions, add test for it (#281)

* Add test MD_MVP (#290)

* Clean up code related to DM restart data (#292)

* Write dm (#291)

* Update use of DM in restart

* remove a redundant assignment in DensityMatrix::DensityMatrix

---------

Co-authored-by: Jean-Luc Fattebert <[email protected]>

Author: dreamer2368

src/DensityMatrix.cc

@@ -12,6 +12,8 @@
 #include "DistMatrix.h"
 #include "MGmol_MPI.h"
 #include "ReplicatedMatrix.h"
+#include "ReplicatedWorkSpace.h"
+#include "hdf_tools.h"
 
 #include <cmath>
 #include <iomanip>
@@ -22,25 +24,32 @@ const double factor_kernel4dot = 10.;
 
 #define PROCRUSTES 0
 
+#define MGMOL_DENSITYMATRIX_FAIL(X)                                            \
+    {                                                                          \
+        std::cerr << "DensityMatrix failure:" << std::endl;                    \
+        std::cerr << "Error Message: " << X << std::endl;                      \
+    }
+
 // occupations in [0,1]
 // DM eigenvalues in [0,orbital_occupation]
 
 template <class MatrixType>
 DensityMatrix<MatrixType>::DensityMatrix(const int ndim)
+    : dim_(ndim),
+      orbitals_index_(-1),
+      occ_uptodate_(false),
+      uniform_occ_(false),
+      stripped_(false)
 {
     assert(ndim > 0);
 
-    dim_ = ndim;
-
     occ_uptodate_ = false;
     stripped_     = false;
     uniform_occ_  = false;
 
     MGmol_MPI& mmpi     = *(MGmol_MPI::instance());
     orbital_occupation_ = mmpi.nspin() > 1 ? 1. : 2.;
 
-    orbitals_index_ = -1;
-
     dm_         = new MatrixType("DM", ndim, ndim);
     kernel4dot_ = new MatrixType("K4dot", ndim, ndim);
     work_       = new MatrixType("work", ndim, ndim);
@@ -438,6 +447,43 @@ void DensityMatrix<MatrixType>::mix(
     orbitals_index_ = new_orbitals_index;
 }
 
+template <class MatrixType>
+int DensityMatrix<MatrixType>::write(HDFrestart& h5f_file, std::string& name)
+{
+    ReplicatedWorkSpace<double>& wspace(
+        ReplicatedWorkSpace<double>::instance());
+
+    wspace.initSquareMatrix(*dm_);
+
+    DISTMATDTYPE* work_matrix = wspace.square_matrix();
+
+    hid_t file_id = h5f_file.file_id();
+    return mgmol_tools::write_matrix(file_id, name, work_matrix, dim_);
+}
+
+template <class MatrixType>
+int DensityMatrix<MatrixType>::read(HDFrestart& h5f_file, std::string& name)
+{
+    ReplicatedWorkSpace<double>& wspace(
+        ReplicatedWorkSpace<double>::instance());
+    DISTMATDTYPE* work_matrix = wspace.square_matrix();
+
+    MGmol_MPI& mmpi = *(MGmol_MPI::instance());
+
+    int ierr = 0;
+    if (mmpi.instancePE0())
+    {
+        hid_t file_id = h5f_file.file_id();
+        ierr          = mgmol_tools::read_matrix(file_id, name, work_matrix);
+    }
+    mmpi.bcast(&ierr, 1);
+
+    if (ierr >= 0) wspace.mpiBcastSquareMatrix();
+    if (ierr >= 0) initMatrix(work_matrix);
+
+    return ierr;
+}
+
 template class DensityMatrix<dist_matrix::DistMatrix<double>>;
 #ifdef HAVE_MAGMA
 template class DensityMatrix<ReplicatedMatrix>;

src/DensityMatrix.h

@@ -10,6 +10,7 @@
 #ifndef MGMOL_DENSITYMATRIX_H
 #define MGMOL_DENSITYMATRIX_H
 
+#include "HDFrestart.h"
 #include "MGmol_MPI.h"
 #include "global.h"
 
@@ -23,7 +24,7 @@
 template <class MatrixType>
 class DensityMatrix
 {
-    int dim_;
+    const int dim_;
     std::vector<double> occupation_;
 
     MatrixType* dm_;
@@ -35,6 +36,10 @@ class DensityMatrix
     bool occ_uptodate_;
     bool uniform_occ_;
     bool stripped_;
+
+    /*!
+     * Max. occupation of an orbital: 1 with spin, 2 otherwise
+     */
     double orbital_occupation_;
 
     DensityMatrix();
@@ -146,6 +151,9 @@ class DensityMatrix
     double getExpectation(const MatrixType& A);
     void mix(
         const double mix, const MatrixType& matA, const int new_orbitals_index);
+
+    int write(HDFrestart& h5f_file, std::string& name);
+    int read(HDFrestart& h5f_file, std::string& name);
 };
 
 #endif

src/ExtendedGridOrbitals.cc

@@ -584,7 +584,6 @@ int ExtendedGridOrbitals::read_hdf5(HDFrestart& h5f_file)
 
     Control& ct = *(Control::instance());
 
-    hid_t file_id    = h5f_file.file_id();
     std::string name = "Function";
     int ierr         = read_func_hdf5(h5f_file, name);
     if (ierr < 0)
@@ -603,7 +602,7 @@ int ExtendedGridOrbitals::read_hdf5(HDFrestart& h5f_file)
     // Read DM
     if (!ct.fullyOccupied())
     {
-        ierr = proj_matrices_->read_dm_hdf5(file_id);
+        ierr = proj_matrices_->readDM(h5f_file);
         if (ierr < 0)
         {
             (*MPIdata::serr)
@@ -618,28 +617,7 @@ int ExtendedGridOrbitals::read_hdf5(HDFrestart& h5f_file)
     return ierr;
 }
 
-int ExtendedGridOrbitals::write_hdf5(
-    HDFrestart& h5f_file, const std::string& name)
-{
-    assert(proj_matrices_ != nullptr);
-    Control& ct = *(Control::instance());
-
-    if (!ct.fullyOccupied())
-    {
-        MGmol_MPI& mmpi(*(MGmol_MPI::instance()));
-        mmpi.barrier();
-
-        int ierr = proj_matrices_->writeDM_hdf5(h5f_file);
-        if (ierr < 0) return ierr;
-    }
-
-    int ierr = write_func_hdf5(h5f_file, name);
-
-    return ierr;
-}
-
-int ExtendedGridOrbitals::write_func_hdf5(
-    HDFrestart& h5f_file, const std::string& name)
+int ExtendedGridOrbitals::write(HDFrestart& h5f_file, const std::string& name)
 {
     if (onpe0)
         (*MPIdata::sout) << "ExtendedGridOrbitals::write_func_hdf5()...\n";
@@ -829,7 +807,7 @@ int ExtendedGridOrbitals::read_func_hdf5(
         const std::string datasetname(getDatasetName(name, icolor));
 
         // check if dataset exists...
-        int err_id = h5f_file.dset_exists(datasetname);
+        int err_id = h5f_file.checkDataExists(datasetname);
         if (h5f_file.gatherDataX()) mmpi.bcast(&err_id, 1);
         if (err_id == 0) break; // dataset does not exists
 

src/ExtendedGridOrbitals.h

@@ -374,8 +374,7 @@ class ExtendedGridOrbitals : public Orbitals
     void multiplyByMatrix2states(const int st1, const int st2,
         const double* mat, ExtendedGridOrbitals& product);
 
-    int write_hdf5(HDFrestart& h5f_file, const std::string& name = "Function");
-    int write_func_hdf5(HDFrestart&, const std::string& name = "Function");
+    int write(HDFrestart&, const std::string& name = "Function");
     int read_hdf5(HDFrestart& h5f_file);
     int read_func_hdf5(HDFrestart&, const std::string& name = "Function");
 

src/HDFrestart.cc

@@ -933,12 +933,12 @@ int HDFrestart::getLRCenters(std::multimap<std::string, Vector3D>& centers,
 
         std::string datasetname(getDatasetName(name, color));
 
-        int err_id = dset_exists(datasetname);
+        int err_id = checkDataExistsLocal(datasetname);
         if (err_id == 0)
         { // dataset does not exists
             // try older version
             datasetname = getDatasetName_old(name, color);
-            err_id      = dset_exists(datasetname);
+            err_id      = checkDataExistsLocal(datasetname);
         }
 
         if (err_id == 0)
@@ -1052,12 +1052,12 @@ int HDFrestart::getLRs(std::shared_ptr<LocalizationRegions> lrs,
 
         std::string datasetname(getDatasetName(name, color));
 
-        int err_id = dset_exists(datasetname);
+        int err_id = checkDataExistsLocal(datasetname);
         if (err_id == 0)
         { // dataset does not exists
             // try older version
             datasetname = getDatasetName_old(name, color);
-            err_id      = dset_exists(datasetname);
+            err_id      = checkDataExistsLocal(datasetname);
         }
         if (err_id == 0)
         { // dataset does not exists

src/HDFrestart.h

@@ -155,28 +155,30 @@ class HDFrestart
         return 0;
     }
 
-    hid_t dset_exists(const std::string& datasetname) const
+    hid_t checkDataExists(const std::string& datasetname) const
     {
-        if (active_)
-        {
-            return dset_exists(datasetname.c_str());
-        }
-        return 0;
+        return checkDataExists(datasetname.c_str());
     }
-    herr_t dset_exists(const char* const datasetname) const
+    herr_t checkDataExists(const char* const datasetname) const
+    {
+        herr_t err_id = checkDataExistsLocal(datasetname);
+
+        short id = (short)err_id;
+        MPI_Bcast(&id, 1, MPI_SHORT, 0, comm_data_);
+        return (herr_t)id;
+    }
+
+    hid_t checkDataExistsLocal(const std::string& datasetname) const
+    {
+        return checkDataExistsLocal(datasetname.c_str());
+    }
+
+    herr_t checkDataExistsLocal(const char* const datasetname) const
     {
         herr_t err_id = 0;
         if (active_)
         {
             err_id = H5LTfind_dataset(file_id_, datasetname);
-
-            if (err_id < 0)
-            {
-                if (onpe0)
-                    (*MPIdata::sout)
-                        << "HDFrestart::dset_exists() failed for dataset "
-                        << datasetname << std::endl;
-            }
         }
         return err_id;
     }

src/LBFGS_IonicStepper.cc

@@ -328,11 +328,11 @@ int LBFGS_IonicStepper::read_lbfgs(HDFrestart& h5f_file)
 
     short check_data = 0;
     std::vector<double> attr_d(16);
+    // Open an existing dataset.
+    std::string datasetname("/LBFGS");
+    int err_id = h5f_file.checkDataExists(datasetname);
     if (onpe0)
     {
-        // Open an existing dataset.
-        std::string datasetname("/LBFGS");
-        int err_id = h5f_file.dset_exists(datasetname);
         if (err_id < 0)
         { // dataset does not exists
             (*MPIdata::sout) << "Warning: no dataset " << datasetname
@@ -343,10 +343,9 @@ int LBFGS_IonicStepper::read_lbfgs(HDFrestart& h5f_file)
             dataset_id = H5Dopen2(file_id, "/LBFGS", H5P_DEFAULT);
             if (dataset_id < 0)
             {
-                if (onpe0)
-                    (*MPIdata::sout) << "Warning: H5Dopen failed for /LBFGS-> "
-                                        "no restart info for LBFGS"
-                                     << std::endl;
+                (*MPIdata::sout) << "Warning: H5Dopen failed for /LBFGS-> "
+                                    "no restart info for LBFGS"
+                                 << std::endl;
             }
             else
             {

src/LocGridOrbitals.cc

@@ -966,7 +966,6 @@ int LocGridOrbitals::read_hdf5(HDFrestart& h5f_file)
 
     Control& ct = *(Control::instance());
 
-    hid_t file_id    = h5f_file.file_id();
     std::string name = "Function";
     int ierr         = read_func_hdf5(h5f_file, name);
     if (ierr < 0)
@@ -984,7 +983,7 @@ int LocGridOrbitals::read_hdf5(HDFrestart& h5f_file)
     // Read DM
     if (!ct.fullyOccupied())
     {
-        ierr = proj_matrices_->read_dm_hdf5(file_id);
+        ierr = proj_matrices_->readDM(h5f_file);
         if (ierr < 0)
         {
             (*MPIdata::serr)
@@ -997,27 +996,7 @@ int LocGridOrbitals::read_hdf5(HDFrestart& h5f_file)
     return ierr;
 }
 
-int LocGridOrbitals::write_hdf5(HDFrestart& h5f_file, const std::string& name)
-{
-    assert(proj_matrices_ != nullptr);
-    Control& ct = *(Control::instance());
-
-    if (!ct.fullyOccupied())
-    {
-        MGmol_MPI& mmpi(*(MGmol_MPI::instance()));
-        mmpi.barrier();
-
-        int ierr = proj_matrices_->writeDM_hdf5(h5f_file);
-        if (ierr < 0) return ierr;
-    }
-
-    int ierr = write_func_hdf5(h5f_file, name);
-
-    return ierr;
-}
-
-int LocGridOrbitals::write_func_hdf5(
-    HDFrestart& h5f_file, const std::string& name)
+int LocGridOrbitals::write(HDFrestart& h5f_file, const std::string& name)
 {
     Control& ct   = *(Control::instance());
     hid_t file_id = h5f_file.file_id();
@@ -1249,7 +1228,7 @@ int LocGridOrbitals::read_func_hdf5(
         const std::string key(itcenter->first);
 
         // checkif dataset exists...
-        int err_id = h5f_file.dset_exists(key);
+        int err_id = h5f_file.checkDataExistsLocal(key);
         if (h5f_file.gatherDataX()) mmpi.bcast(&err_id, 1);
         if (err_id == 0) break; // dataset does not exists
 

src/LocGridOrbitals.h

@@ -410,8 +410,7 @@ class LocGridOrbitals : public Orbitals
     void multiplyByMatrix2states(const int st1, const int st2,
         const double* mat, LocGridOrbitals& product);
 
-    int write_hdf5(HDFrestart& h5f_file, const std::string& name = "Function");
-    int write_func_hdf5(HDFrestart&, const std::string& name = "Function");
+    int write(HDFrestart&, const std::string& name = "Function");
     int read_hdf5(HDFrestart& h5f_file);
     int read_func_hdf5(HDFrestart&, const std::string& name = "Function");
 

src/MGmol.h

@@ -133,7 +133,7 @@ class MGmol : public MGmolInterface
         KBPsiMatrixSparse* kbpsi, dist_matrix::DistMatrix<DISTMATDTYPE>& hij);
     void computeHnlPhiAndAdd2HPhi(Ions& ions, OrbitalsType& phi,
         OrbitalsType& hphi, const KBPsiMatrixSparse* const kbpsi);
-    int dumpMDrestartFile(OrbitalsType** orbitals, Ions& ions,
+    int dumpMDrestartFile(OrbitalsType& orbitals, Ions& ions,
         Rho<OrbitalsType>& rho, const bool write_extrapolated_wf,
         const short count);