Merge release

scripts/build_condo-mod.sh

@@ -3,9 +3,6 @@
 #Before compiling, load the following modules:
 source scripts/modules.condo-mod
 
-# We need to define the cmake blas vendor option here to find the right one.
-BLAS_VENDOR=OpenBLAS
-
 MGMOL_ROOT=`pwd`
 
 INSTALL_DIR=${MGMOL_ROOT}/install
@@ -15,16 +12,17 @@ BUILD_DIR=${MGMOL_ROOT}/build
 mkdir -p ${BUILD_DIR}
 cd ${BUILD_DIR}
 
+SCALAPACK_DIR=/home/q8j/Software/ScaLapack/scalapack-2.2.2
+
 # call cmake 
 cmake -DCMAKE_INSTALL_PREFIX=${INSTALL_DIR} \
       -DCMAKE_BUILD_TYPE=Release \
       -DCMAKE_CXX_COMPILER=mpiCC \
       -DCMAKE_Fortran_COMPILER=mpif77 \
-      -DBLA_VENDOR=${BLAS_VENDOR} \
       -DMGMOL_USE_HDF5P=OFF \
       -DMGMOL_WITH_CLANG_FORMAT=ON \
       -DCMAKE_PREFIX_PATH=${HOME}/bin \
-      -DSCALAPACK_LIBRARY="${SCALAPACK_DIR}/lib/libscalapack.a;/lib64/libgfortran.so.3" \
+      -DSCALAPACK_LIBRARY="${SCALAPACK_DIR}/lib/libscalapack.a;/lib64/libgfortran.so.5" \
       -DMPIEXEC_EXECUTABLE=${OPENMPI_DIR}/bin/mpiexec \
       ..
 

scripts/build_condo.sh

@@ -1,41 +1,30 @@
 #/bin/bash
 ## An example script to build on ONRL condo systems (CADES).
-## This script assumes intel/ mkl libraries are being used.
 
 #Before compiling, load the following modules:
 source scripts/modules.condo
 
-# set some environment variables using loaded module path
-export SCALAPACK_ROOT=${MKLROOT}
-
-# We need to define the cmake blas vendor option here to find the right one.
-BLAS_VENDOR=Intel10_64lp
-
-# manually set the location of BLACS libraries for scalapack
-BLACS_LIB=${MKLROOT}/lib/intel64
-
 MGMOL_ROOT=`pwd`
 
-INSTALL_DIR=${MGMOL_ROOT}/mgmol_install
+INSTALL_DIR=${MGMOL_ROOT}/install
 mkdir -p ${INSTALL_DIR}
 
 BUILD_DIR=${MGMOL_ROOT}/build
 mkdir -p ${BUILD_DIR}
 cd ${BUILD_DIR}
 
+SCALAPACK_DIR=/home/q8j/Software/ScaLapack/scalapack-2.2.2
+
 # call cmake 
 cmake -DCMAKE_INSTALL_PREFIX=${INSTALL_DIR} \
       -DCMAKE_BUILD_TYPE=Release \
       -DCMAKE_CXX_COMPILER=mpiCC \
       -DCMAKE_Fortran_COMPILER=mpif77 \
-      -DHDF5_LIBRARIES=${HDF5_DIR}/lib/libhdf5.so \
-      -DHDF5_HL_LIBRARIES=${HDF5_DIR}/lib/libhdf5_hl.so \
-      -DHDF5_INCLUDE_DIRS=${HDF5_DIR}/include \
-      -DBLA_VENDOR=${BLAS_VENDOR} \
+      -DMGMOL_USE_HDF5P=OFF \
       -DMGMOL_WITH_CLANG_FORMAT=ON \
       -DCMAKE_PREFIX_PATH=${HOME}/bin \
+      -DSCALAPACK_LIBRARY="${SCALAPACK_DIR}/lib/libscalapack.a;/lib64/libgfortran.so.5" \
       -DMPIEXEC_EXECUTABLE=${OPENMPI_DIR}/bin/mpiexec \
-      -DSCALAPACK_BLACS_LIBRARY=${BLACS_LIB}/libmkl_blacs_openmpi_lp64.so \
       ..
 
 # call make install

scripts/modules.condo

@@ -1,6 +1,5 @@
-module load PE-intel/3.0
-module load boost/1.67.0-pe3
-module load mkl
-module load hdf5_parallel/1.10.3
-module load cmake/3.18.4
+module load hdf5
+module load boost
+module load cmake
 module load python
+module load openblas

scripts/modules.condo-mod

@@ -1,7 +1,5 @@
-module load PE-gnu/3.0
-module load hdf5-parallel/1.8.20
+module load hdf5
 module load boost
-module load cmake/3.20.3
-module load python/3.6.6
-module load openBLAS/0.2.19
-module load scalapack/2.0.2
+module load cmake
+module load python
+module load openblas

src/Control.cc

@@ -224,7 +224,8 @@ void Control::print(std::ostream& os)
        << conv_tol << std::endl;
     os << std::fixed;
     os << " Density matrix mixing = " << dm_mix << std::endl;
-    os << std::setprecision(4) << std::scientific << " Density matrix tol = " << dm_tol << std::endl;
+    os << std::setprecision(4) << std::scientific
+       << " Density matrix tol = " << dm_tol << std::endl;
     if (DMEigensolver() == DMEigensolverType::Eigensolver)
     {
         os << " Density matrix computation algorithm = "

src/DavidsonSolver.cc

@@ -206,8 +206,8 @@ double DavidsonSolver<OrbitalsType, MatrixType>::evaluateDerivative(
     const double dbeta = 0.0001;
     *work2N_           = dm2Ninit;
     work2N_->axpy(dbeta, delta_dm);
-    // proj_mat2N_->setDM(*work2N_,orbitals.getIterativeIndex());
-    proj_mat2N_->setDM(*work2N_, -1);
+
+    proj_mat2N_->setDM(*work2N_);
     proj_mat2N_->computeOccupationsFromDM();
 
     const double tsd0e = evalEntropy(proj_mat2N_.get(), false, os_);
@@ -250,7 +250,7 @@ void DavidsonSolver<OrbitalsType, MatrixType>::buildTarget2N_MVP(
     // if( onpe0 )os_<<"Compute X2N..."<<endl;
     proj_mat2N_->setHB2H();
 
-    proj_mat2N_->updateDM(orbitals_index);
+    proj_mat2N_->updateDM();
 
     target = proj_mat2N_->dm();
 
@@ -584,7 +584,7 @@ int DavidsonSolver<OrbitalsType, MatrixType>::solve(
                 //
                 if (mmpi.PE0() && ct.verbose > 2)
                     os_ << "Target energy..." << std::endl;
-                proj_mat2N_->setDM(target, orbitals.getIterativeIndex());
+                proj_mat2N_->setDM(target);
                 proj_mat2N_->computeOccupationsFromDM();
                 double nel = proj_mat2N_->getNel();
                 if (mmpi.PE0() && ct.verbose > 2)
@@ -643,7 +643,7 @@ int DavidsonSolver<OrbitalsType, MatrixType>::solve(
             // update DM
             *work2N_ = dm2Ninit;
             work2N_->axpy(beta, delta_dm);
-            proj_mat2N_->setDM(*work2N_, orbitals.getIterativeIndex());
+            proj_mat2N_->setDM(*work2N_);
 
             if (inner_it < ct.dm_inner_steps - 1)
             {
@@ -809,7 +809,7 @@ int DavidsonSolver<OrbitalsType, MatrixType>::solve(
             = dynamic_cast<ProjectedMatrices<MatrixType>*>(
                 orbitals.getProjMatrices());
         assert(pmat);
-        pmat->buildDM(new_occ, orbitals.getIterativeIndex());
+        pmat->buildDM(new_occ);
 
         if (retval == 0) break;
 

src/DensityMatrix.cc

@@ -22,7 +22,11 @@
 
 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;                      \
+    }
 
 #define MGMOL_DENSITYMATRIX_FAIL(X)                                            \
     {                                                                          \
@@ -36,7 +40,7 @@ const double factor_kernel4dot = 10.;
 template <class MatrixType>
 DensityMatrix<MatrixType>::DensityMatrix(const int ndim)
     : dim_(ndim),
-      orbitals_index_(-1),
+      update_index_(-1),
       occ_uptodate_(false),
       uniform_occ_(false),
       stripped_(false)
@@ -66,8 +70,8 @@ DensityMatrix<MatrixType>::~DensityMatrix()
 }
 
 template <class MatrixType>
-void DensityMatrix<MatrixType>::build(const MatrixType& zmat,
-    const std::vector<double>& occ, const int new_orbitals_index)
+void DensityMatrix<MatrixType>::build(
+    const MatrixType& zmat, const std::vector<double>& occ)
 {
 #ifdef PRINT_OPERATIONS
     MGmol_MPI& mmpi = *(MGmol_MPI::instance());
@@ -98,32 +102,30 @@ void DensityMatrix<MatrixType>::build(const MatrixType& zmat,
     work_->symm('r', 'l', 1., gamma, zmat, 0.);
     kernel4dot_->gemm('n', 't', 1., *work_, zmat, 0.);
 
-    stripped_       = false;
-    orbitals_index_ = new_orbitals_index;
+    stripped_ = false;
+    update_index_++;
 }
 
 template <class MatrixType>
-void DensityMatrix<MatrixType>::build(
-    const MatrixType& zmat, const int new_orbitals_index)
+void DensityMatrix<MatrixType>::build(const MatrixType& zmat)
 {
-    build(zmat, occupation_, new_orbitals_index);
+    build(zmat, occupation_);
 }
 
 // build diagonal matrix
 template <class MatrixType>
-void DensityMatrix<MatrixType>::build(
-    const std::vector<double>& occ, const int new_orbitals_index)
+void DensityMatrix<MatrixType>::build(const std::vector<double>& occ)
 {
     assert(dm_ != nullptr);
     assert(!occ.empty());
 
     setOccupations(occ);
 
-    build(new_orbitals_index);
+    build();
 }
 
 template <class MatrixType>
-void DensityMatrix<MatrixType>::build(const int new_orbitals_index)
+void DensityMatrix<MatrixType>::build()
 {
     MGmol_MPI& mmpi = *(MGmol_MPI::instance());
 #ifdef PRINT_OPERATIONS
@@ -148,13 +150,12 @@ void DensityMatrix<MatrixType>::build(const int new_orbitals_index)
                         * std::min(1., factor_kernel4dot * occupation_[i]));
     kernel4dot_->setDiagonal(w);
 
-    stripped_       = false;
-    orbitals_index_ = new_orbitals_index;
+    stripped_ = false;
+    update_index_++;
 }
 
 template <class MatrixType>
-void DensityMatrix<MatrixType>::setUniform(
-    const double nel, const int new_orbitals_index)
+void DensityMatrix<MatrixType>::setUniform(const double nel)
 {
     assert(!occupation_.empty());
 
@@ -171,22 +172,13 @@ void DensityMatrix<MatrixType>::setUniform(
 
     uniform_occ_ = true;
 
-    build(occupation_, new_orbitals_index);
-}
-
-template <class MatrixType>
-void DensityMatrix<MatrixType>::buildFromBlock(const MatrixType& block00)
-{
-    dm_->clear();
-    dm_->assign(block00, 0, 0);
-    dm_->print(std::cout, 0, 0, 25, 25);
+    build(occupation_);
 }
 
 template <class MatrixType>
 void DensityMatrix<MatrixType>::rotate(
     const MatrixType& rotation_matrix, const bool flag_eigen)
 {
-
     if (!flag_eigen)
     {
         MatrixType invU(rotation_matrix);
@@ -204,6 +196,8 @@ void DensityMatrix<MatrixType>::rotate(
         invU.getrs('n', tmp, ipiv);
 
         *dm_ = tmp;
+
+        update_index_++;
     }
 }
 
@@ -370,8 +364,7 @@ double DensityMatrix<MatrixType>::computeEntropy() const
 }
 
 template <class MatrixType>
-void DensityMatrix<MatrixType>::setto2InvS(
-    const MatrixType& invS, const int orbitals_index)
+void DensityMatrix<MatrixType>::setto2InvS(const MatrixType& invS)
 {
     *dm_ = invS;
     dm_->scal(orbital_occupation_);
@@ -382,8 +375,8 @@ void DensityMatrix<MatrixType>::setto2InvS(
             occupation_[st] = 1.;
         occ_uptodate_ = true;
     }
-    uniform_occ_    = false;
-    orbitals_index_ = orbitals_index;
+    uniform_occ_ = false;
+    update_index_++;
 }
 
 template <class MatrixType>
@@ -400,8 +393,7 @@ void DensityMatrix<MatrixType>::stripS(const MatrixType& ls)
 }
 
 template <class MatrixType>
-void DensityMatrix<MatrixType>::dressUpS(
-    const MatrixType& ls, const int new_orbitals_index)
+void DensityMatrix<MatrixType>::dressUpS(const MatrixType& ls)
 {
     assert(stripped_);
 
@@ -410,10 +402,11 @@ void DensityMatrix<MatrixType>::dressUpS(
     *dm_ = *work_;
     ls.trtrs('l', 't', 'n', *dm_);
 
-    orbitals_index_ = new_orbitals_index;
-    occ_uptodate_   = false;
-    uniform_occ_    = false;
-    stripped_       = false;
+    update_index_++;
+
+    occ_uptodate_ = false;
+    uniform_occ_  = false;
+    stripped_     = false;
 }
 
 // dm_ -> u*dm_*u^T
@@ -424,6 +417,8 @@ void DensityMatrix<MatrixType>::transform(const MatrixType& u)
 {
     work_->gemm('n', 't', 1., *dm_, u, 0.);
     dm_->gemm('n', 'n', 1., u, *work_, 0.);
+
+    update_index_++;
 }
 
 template <class MatrixType>
@@ -434,13 +429,21 @@ double DensityMatrix<MatrixType>::getExpectation(const MatrixType& A)
 }
 
 template <class MatrixType>
-void DensityMatrix<MatrixType>::mix(
-    const double mix, const MatrixType& matA, const int new_orbitals_index)
+void DensityMatrix<MatrixType>::mix(const double mix, const MatrixType& matA)
 {
     dm_->scal(1. - mix);
-
     dm_->axpy(mix, matA);
-    orbitals_index_ = new_orbitals_index;
+
+    update_index_++;
+}
+
+template <class MatrixType>
+void DensityMatrix<MatrixType>::linearExtrapolate(const MatrixType& previous_dm)
+{
+    dm_->scal(2.);
+    dm_->axpy(-1., previous_dm);
+
+    update_index_++;
 }
 
 template <class MatrixType>