diff --git a/scripts/build_condo-mod.sh b/scripts/build_condo-mod.sh index d9abe034..4dfe5fc6 100755 --- a/scripts/build_condo-mod.sh +++ b/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 \ .. diff --git a/scripts/build_condo.sh b/scripts/build_condo.sh index 3a349cbe..68aa16dc 100755 --- a/scripts/build_condo.sh +++ b/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 diff --git a/scripts/modules.condo b/scripts/modules.condo index 6f27e62b..de36562f 100644 --- a/scripts/modules.condo +++ b/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 diff --git a/scripts/modules.condo-mod b/scripts/modules.condo-mod index 18377905..de36562f 100644 --- a/scripts/modules.condo-mod +++ b/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 diff --git a/src/Control.cc b/src/Control.cc index cbf2eadd..980c1abe 100644 --- a/src/Control.cc +++ b/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 = " diff --git a/src/DavidsonSolver.cc b/src/DavidsonSolver.cc index a049fde7..c34b6396 100644 --- a/src/DavidsonSolver.cc +++ b/src/DavidsonSolver.cc @@ -206,8 +206,8 @@ double DavidsonSolver::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::buildTarget2N_MVP( // if( onpe0 )os_<<"Compute X2N..."<setHB2H(); - proj_mat2N_->updateDM(orbitals_index); + proj_mat2N_->updateDM(); target = proj_mat2N_->dm(); @@ -584,7 +584,7 @@ int DavidsonSolver::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::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::solve( = dynamic_cast*>( orbitals.getProjMatrices()); assert(pmat); - pmat->buildDM(new_occ, orbitals.getIterativeIndex()); + pmat->buildDM(new_occ); if (retval == 0) break; diff --git a/src/DensityMatrix.cc b/src/DensityMatrix.cc index 78d7008e..83d715d5 100644 --- a/src/DensityMatrix.cc +++ b/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 DensityMatrix::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::~DensityMatrix() } template -void DensityMatrix::build(const MatrixType& zmat, - const std::vector& occ, const int new_orbitals_index) +void DensityMatrix::build( + const MatrixType& zmat, const std::vector& occ) { #ifdef PRINT_OPERATIONS MGmol_MPI& mmpi = *(MGmol_MPI::instance()); @@ -98,32 +102,30 @@ void DensityMatrix::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 -void DensityMatrix::build( - const MatrixType& zmat, const int new_orbitals_index) +void DensityMatrix::build(const MatrixType& zmat) { - build(zmat, occupation_, new_orbitals_index); + build(zmat, occupation_); } // build diagonal matrix template -void DensityMatrix::build( - const std::vector& occ, const int new_orbitals_index) +void DensityMatrix::build(const std::vector& occ) { assert(dm_ != nullptr); assert(!occ.empty()); setOccupations(occ); - build(new_orbitals_index); + build(); } template -void DensityMatrix::build(const int new_orbitals_index) +void DensityMatrix::build() { MGmol_MPI& mmpi = *(MGmol_MPI::instance()); #ifdef PRINT_OPERATIONS @@ -148,13 +150,12 @@ void DensityMatrix::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 -void DensityMatrix::setUniform( - const double nel, const int new_orbitals_index) +void DensityMatrix::setUniform(const double nel) { assert(!occupation_.empty()); @@ -171,22 +172,13 @@ void DensityMatrix::setUniform( uniform_occ_ = true; - build(occupation_, new_orbitals_index); -} - -template -void DensityMatrix::buildFromBlock(const MatrixType& block00) -{ - dm_->clear(); - dm_->assign(block00, 0, 0); - dm_->print(std::cout, 0, 0, 25, 25); + build(occupation_); } template void DensityMatrix::rotate( const MatrixType& rotation_matrix, const bool flag_eigen) { - if (!flag_eigen) { MatrixType invU(rotation_matrix); @@ -204,6 +196,8 @@ void DensityMatrix::rotate( invU.getrs('n', tmp, ipiv); *dm_ = tmp; + + update_index_++; } } @@ -370,8 +364,7 @@ double DensityMatrix::computeEntropy() const } template -void DensityMatrix::setto2InvS( - const MatrixType& invS, const int orbitals_index) +void DensityMatrix::setto2InvS(const MatrixType& invS) { *dm_ = invS; dm_->scal(orbital_occupation_); @@ -382,8 +375,8 @@ void DensityMatrix::setto2InvS( occupation_[st] = 1.; occ_uptodate_ = true; } - uniform_occ_ = false; - orbitals_index_ = orbitals_index; + uniform_occ_ = false; + update_index_++; } template @@ -400,8 +393,7 @@ void DensityMatrix::stripS(const MatrixType& ls) } template -void DensityMatrix::dressUpS( - const MatrixType& ls, const int new_orbitals_index) +void DensityMatrix::dressUpS(const MatrixType& ls) { assert(stripped_); @@ -410,10 +402,11 @@ void DensityMatrix::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::transform(const MatrixType& u) { work_->gemm('n', 't', 1., *dm_, u, 0.); dm_->gemm('n', 'n', 1., u, *work_, 0.); + + update_index_++; } template @@ -434,13 +429,21 @@ double DensityMatrix::getExpectation(const MatrixType& A) } template -void DensityMatrix::mix( - const double mix, const MatrixType& matA, const int new_orbitals_index) +void DensityMatrix::mix(const double mix, const MatrixType& matA) { dm_->scal(1. - mix); - dm_->axpy(mix, matA); - orbitals_index_ = new_orbitals_index; + + update_index_++; +} + +template +void DensityMatrix::linearExtrapolate(const MatrixType& previous_dm) +{ + dm_->scal(2.); + dm_->axpy(-1., previous_dm); + + update_index_++; } template diff --git a/src/DensityMatrix.h b/src/DensityMatrix.h index 1d09e79e..8ce74542 100644 --- a/src/DensityMatrix.h +++ b/src/DensityMatrix.h @@ -31,7 +31,10 @@ class DensityMatrix MatrixType* kernel4dot_; MatrixType* work_; - int orbitals_index_; + /*! + * Keep track of changes, incremented every time dm_ is updated + */ + int update_index_; bool occ_uptodate_; bool uniform_occ_; @@ -46,16 +49,16 @@ class DensityMatrix DensityMatrix& operator=(const DensityMatrix&); DensityMatrix(const DensityMatrix&); - void build(const int new_orbitals_index); + void build(); public: DensityMatrix(const int ndim); ~DensityMatrix(); - void setUniform(const double nel, const int new_orbitals_index); + void setUniform(const double nel); - int getOrbitalsIndex() const { return orbitals_index_; } + int getIndex() const { return update_index_; } bool occupationsUptodate() const { return occ_uptodate_; } bool fromUniformOccupations() const { return uniform_occ_; } @@ -84,10 +87,10 @@ class DensityMatrix const MatrixType& kernel4dot() const { return *kernel4dot_; } - void setMatrix(const MatrixType& mat, const int orbitals_index) + void setMatrix(const MatrixType& mat) { - *dm_ = mat; - orbitals_index_ = orbitals_index; + *dm_ = mat; + update_index_++; setDummyOcc(); @@ -112,7 +115,7 @@ class DensityMatrix uniform_occ_ = false; stripped_ = false; - orbitals_index_ = 0; + update_index_ = 0; } void getOccupations(std::vector& occ) const @@ -126,22 +129,19 @@ class DensityMatrix void setOccupations(const std::vector& occ); - void setto2InvS(const MatrixType& invS, const int orbitals_index); + void setto2InvS(const MatrixType& invS); void stripS(const MatrixType& ls); - void dressUpS(const MatrixType& ls, const int new_orbitals_index); + void dressUpS(const MatrixType& ls); // dm_ -> u*dm_*u^T void transform(const MatrixType& u); - void buildFromBlock(const MatrixType& block00); - double computeEntropy() const; void computeOccupations(const MatrixType& ls); - void build(const std::vector& occ, const int new_orbitals_index); - void build(const MatrixType& z, const int new_orbitals_index); - void build(const MatrixType& z, const std::vector& occ, - const int new_orbitals_index); + void build(const std::vector& occ); + void build(const MatrixType& z); + void build(const MatrixType& z, const std::vector& occ); void rotate(const MatrixType& rotation_matrix, const bool flag_eigen); void printOccupations(std::ostream& os) const; @@ -149,8 +149,12 @@ class DensityMatrix void diagonalize( const char eigv, std::vector& occ, MatrixType& vect); double getExpectation(const MatrixType& A); - void mix( - const double mix, const MatrixType& matA, const int new_orbitals_index); + void mix(const double mix, const MatrixType& matA); + + /*! + * dm <- dm + (dm-previous_dm) = 2.*dm - previous_dm + */ + void linearExtrapolate(const MatrixType& previous_dm); int write(HDFrestart& h5f_file, std::string& name); int read(HDFrestart& h5f_file, std::string& name); diff --git a/src/DensityMatrixSparse.cc b/src/DensityMatrixSparse.cc index 4217b8ac..e554501a 100644 --- a/src/DensityMatrixSparse.cc +++ b/src/DensityMatrixSparse.cc @@ -26,7 +26,7 @@ DensityMatrixSparse::DensityMatrixSparse( MGmol_MPI& mmpi = *(MGmol_MPI::instance()); orbital_occupation_ = mmpi.nspin() > 1 ? 1. : 2.; - orbitals_index_ = -1; + update_index_ = -1; if (dim_ > 0) { @@ -44,27 +44,24 @@ DensityMatrixSparse::~DensityMatrixSparse() } } -void DensityMatrixSparse::setUniform(const double nel, const int orbitals_index) +void DensityMatrixSparse::setUniform(const double nel) { const double occ = (double)((double)nel / (double)dim_); assert(occ < 1.01); - orbitals_index_ = orbitals_index; + update_index_++; (*dm_).reset(); const double uval = (double)occ * orbital_occupation_; for (std::vector::const_iterator st = locvars_.begin(); st != locvars_.end(); ++st) (*dm_).insertMatrixElement(*st, *st, uval, INSERT, true); - - return; } -void DensityMatrixSparse::setto2InvS( - const VariableSizeMatrix& invS, const int orbitals_index) +void DensityMatrixSparse::setto2InvS(const VariableSizeMatrix& invS) { *dm_ = invS; dm_->scale(orbital_occupation_); - orbitals_index_ = orbitals_index; + update_index_++; } // build density matrix, given computed locally centered data void DensityMatrixSparse::assembleMatrixFromCenteredData( @@ -92,7 +89,7 @@ void DensityMatrixSparse::assembleMatrixFromCenteredData( dtor_DM.updateLocalRows((*dm_)); gather_DM_tm_.stop(); - orbitals_index_ = orbitals_index; + update_index_ = orbitals_index; } // compute trace of dot product dm_ . vsmat double DensityMatrixSparse::getTraceDotProductWithMat( diff --git a/src/DensityMatrixSparse.h b/src/DensityMatrixSparse.h index b574bd0b..170057fa 100644 --- a/src/DensityMatrixSparse.h +++ b/src/DensityMatrixSparse.h @@ -32,7 +32,7 @@ class DensityMatrixSparse VariableSizeMatrix* dm_; - int orbitals_index_; + int update_index_; double orbital_occupation_; @@ -43,16 +43,14 @@ class DensityMatrixSparse ~DensityMatrixSparse(); - void setUniform(const double nel, const int new_orbitals_index); + void setUniform(const double nel); - void setto2InvS( - const VariableSizeMatrix& invS, const int orbitals_index); - int getOrbitalsIndex() const { return orbitals_index_; } - void setMatrix( - const VariableSizeMatrix& mat, const int orbitals_index) + void setto2InvS(const VariableSizeMatrix& invS); + int getIndex() const { return update_index_; } + void setMatrix(const VariableSizeMatrix& mat) { - *dm_ = mat; - orbitals_index_ = orbitals_index; + *dm_ = mat; + update_index_++; } void assembleMatrixFromCenteredData(const std::vector& data, const std::vector& locRowIds, const int* globalColIds, diff --git a/src/EigenDMStrategy.cc b/src/EigenDMStrategy.cc index 8fd8443a..358ed5f9 100644 --- a/src/EigenDMStrategy.cc +++ b/src/EigenDMStrategy.cc @@ -37,7 +37,7 @@ int EigenDMStrategy::update(OrbitalsType& orbitals) = dynamic_cast< ProjectedMatrices>*>( proj_matrices_); - pmat->updateDMwithEigenstatesAndRotate(orbitals.getIterativeIndex(), zz); + pmat->updateDMwithEigenstatesAndRotate(zz); // if( onpe0 && ct.verbose>2 ) // (*MPIdata::sout)<<"get_dm_diag: rotate orbitals "<::solve( ProjMatrixType* projmatrices = dynamic_cast(orbitals.getProjMatrices()); - int iterative_index = 0; - // save computed vh for a fair energy "comparison" with vh computed // in close neigborhood const pb::GridFunc vh_init(electrostat_->getVh()); @@ -134,13 +132,11 @@ int HamiltonianMVPSolver::solve( // // evaluate energy at origin // - iterative_index++; - projmatrices->assignH(*hmatrix_); projmatrices->setHB2H(); // update DM and compute entropy - projmatrices->updateDM(iterative_index); + projmatrices->updateDM(); double ts0 = evalEntropyMVP(projmatrices, true, os_); // Update density rho_->update(orbitals); @@ -168,10 +164,8 @@ int HamiltonianMVPSolver::solve( // MatrixType htarget(projmatrices->getH()); - iterative_index++; - // update DM and compute entropy - projmatrices->updateDM(iterative_index); + projmatrices->updateDM(); double ts1 = evalEntropyMVP(projmatrices, true, os_); // Update density rho_->update(orbitals); @@ -204,10 +198,8 @@ int HamiltonianMVPSolver::solve( projmatrices->assignH(h11); projmatrices->setHB2H(); - iterative_index++; - // update DM and entropy - projmatrices->updateDM(iterative_index); + projmatrices->updateDM(); double tsi = evalEntropyMVP(projmatrices, true, os_); // Update density @@ -262,10 +254,8 @@ int HamiltonianMVPSolver::solve( projmatrices->assignH(h11); projmatrices->setHB2H(); - iterative_index++; - // update DM and entropy - projmatrices->updateDM(iterative_index); + projmatrices->updateDM(); tsi = evalEntropyMVP(projmatrices, true, os_); // Update density @@ -325,9 +315,7 @@ int HamiltonianMVPSolver::solve( projmatrices->assignH(*hmatrix_); projmatrices->setHB2H(); - iterative_index++; - - projmatrices->updateDM(iterative_index); + projmatrices->updateDM(); // Generate new density rho_->update(orbitals); diff --git a/src/MGmol.cc b/src/MGmol.cc index ec780d0e..ac022d91 100644 --- a/src/MGmol.cc +++ b/src/MGmol.cc @@ -365,9 +365,7 @@ int MGmol::initial() // initialize Rho if (ct.verbose > 0) printWithTimeStamp("Initialize Rho...", os_); - if (ct.restart_info <= 1) - proj_matrices_->setDMuniform( - ct.getNelSpin(), current_orbitals_->getIterativeIndex()); + if (ct.restart_info <= 1) proj_matrices_->setDMuniform(ct.getNelSpin()); rho_->setup(ct.getOrthoType(), current_orbitals_->getOverlappingGids()); diff --git a/src/MVPSolver.cc b/src/MVPSolver.cc index 1c093e53..4f9165a5 100644 --- a/src/MVPSolver.cc +++ b/src/MVPSolver.cc @@ -104,8 +104,8 @@ double MVPSolver::evaluateDerivative( const double dbeta = 0.001; *work_ = dmInit; work_->axpy(dbeta, delta_dm); - // proj_mat_work_->setDM(*work_,orbitals.getIterativeIndex()); - proj_mat_work_->setDM(*work_, -1); + + proj_mat_work_->setDM(*work_); proj_mat_work_->computeOccupationsFromDM(); const double tsd0e = evalEntropyMVP(proj_mat_work_, false, os_); @@ -145,7 +145,7 @@ void MVPSolver::buildTarget_MVP( // proj_mat_work_->setHB2H(); - proj_mat_work_->updateDM(orbitals_index); + proj_mat_work_->updateDM(); target = proj_mat_work_->dm(); if (ct.verbose > 2) @@ -301,7 +301,7 @@ int MVPSolver::solve(OrbitalsType& orbitals) // if (onpe0 && ct.verbose > 2) std::cout << "MVP --- Target energy..." << std::endl; - proj_mat_work_->setDM(target, orbitals.getIterativeIndex()); + proj_mat_work_->setDM(target); proj_mat_work_->computeOccupationsFromDM(); if (ct.verbose > 2) proj_mat_work_->printOccupations(os_); const double nel = proj_mat_work_->getNel(); @@ -353,7 +353,7 @@ int MVPSolver::solve(OrbitalsType& orbitals) { *work_ = target; } - proj_mat_work_->setDM(*work_, orbitals.getIterativeIndex()); + proj_mat_work_->setDM(*work_); if (inner_it < n_inner_steps_ - 1) { @@ -373,7 +373,7 @@ int MVPSolver::solve(OrbitalsType& orbitals) ProjectedMatrices* projmatrices = dynamic_cast*>( orbitals.getProjMatrices()); - projmatrices->setDM(*work_, orbitals.getIterativeIndex()); + projmatrices->setDM(*work_); projmatrices->setEigenvalues(proj_mat_work_->getEigenvalues()); projmatrices->assignH(proj_mat_work_->getH()); projmatrices->setHB2H(); diff --git a/src/NonOrthoDMStrategy.cc b/src/NonOrthoDMStrategy.cc index 3835b039..28a8f1d5 100644 --- a/src/NonOrthoDMStrategy.cc +++ b/src/NonOrthoDMStrategy.cc @@ -31,7 +31,7 @@ void NonOrthoDMStrategy::initialize(OrbitalsType& orbitals) (*MPIdata::sout) << "NonOrthoDMStrategy::initialize()..." << std::endl; } - proj_matrices_->updateDM(orbitals.getIterativeIndex()); + proj_matrices_->updateDM(); } template @@ -59,11 +59,11 @@ int NonOrthoDMStrategy::update(OrbitalsType& orbitals) if (mix_ < 1.) proj_matrices_->saveDM(); // compute new density matrix - proj_matrices_->updateDM(orbitals.getIterativeIndex()); + proj_matrices_->updateDM(); if (mix_ < 1.) { - proj_matrices_->updateDMwithRelax(mix_, orbitals.getIterativeIndex()); + proj_matrices_->updateDMwithRelax(mix_); } if (ct.verbose > 2) diff --git a/src/ProjectedMatrices.cc b/src/ProjectedMatrices.cc index 74ff01ed..87c9fd48 100644 --- a/src/ProjectedMatrices.cc +++ b/src/ProjectedMatrices.cc @@ -271,7 +271,7 @@ void ProjectedMatrices::setDMto2InvS() if (mmpi.instancePE0() && ct.verbose > 1) std::cout << "ProjectedMatrices::setDMto2InvS()..." << std::endl; - dm_->setto2InvS(gm_->getInverse(), gm_->getAssociatedOrbitalsIndex()); + dm_->setto2InvS(gm_->getInverse()); } template @@ -296,30 +296,28 @@ void ProjectedMatrices::solveGenEigenProblem( } template -void ProjectedMatrices::buildDM( - const MatrixType& z, const int orbitals_index) +void ProjectedMatrices::buildDM(const MatrixType& z) { - dm_->build(z, orbitals_index); + dm_->build(z); } template -void ProjectedMatrices::buildDM(const MatrixType& z, - const std::vector& occ, const int orbitals_index) +void ProjectedMatrices::buildDM( + const MatrixType& z, const std::vector& occ) { - dm_->build(z, occ, orbitals_index); + dm_->build(z, occ); } template void ProjectedMatrices::buildDM( - const std::vector& occ, const int orbitals_index) + const std::vector& occ) { - dm_->build(occ, orbitals_index); + dm_->build(occ); } // Use Chebyshev approximation to compute chemical potential and density matrix template -void ProjectedMatrices::updateDMwithChebApproximation( - const int iterative_index) +void ProjectedMatrices::updateDMwithChebApproximation() { MGmol_MPI& mmpi = *(MGmol_MPI::instance()); Control& ct = *(Control::instance()); @@ -354,15 +352,14 @@ void ProjectedMatrices::updateDMwithChebApproximation( } // compute chemical potential and density matrix with Chebyshev // approximation. - double final_mu = computeChemicalPotentialAndDMwithChebyshev( - order, emin, emax, iterative_index); + double final_mu + = computeChemicalPotentialAndDMwithChebyshev(order, emin, emax); if (mmpi.instancePE0() && ct.verbose > 1) std::cout << "Final mu_ = " << final_mu << " [Ha]" << std::endl; } template -void ProjectedMatrices::updateDMwithEigenstates( - const int iterative_index) +void ProjectedMatrices::updateDMwithEigenstates() { MGmol_MPI& mmpi = *(MGmol_MPI::instance()); Control& ct = *(Control::instance()); @@ -381,14 +378,14 @@ void ProjectedMatrices::updateDMwithEigenstates( // Build the density matrix X // X = Z * gamma * Z^T - buildDM(zz, iterative_index); + buildDM(zz); } //"replicated" implementation of SP2. // Theta is replicated on each MPI task, and SP2 solve run independently // by each MPI task template -void ProjectedMatrices::updateDMwithSP2(const int iterative_index) +void ProjectedMatrices::updateDMwithSP2() { MGmol_MPI& mmpi = *(MGmol_MPI::instance()); Control& ct = *(Control::instance()); @@ -433,21 +430,21 @@ void ProjectedMatrices::updateDMwithSP2(const int iterative_index) MatrixType dm("dm", dim_, dim_); sp2.getDM(dm, gm_->getInverse()); - dm_->setMatrix(dm, iterative_index); + dm_->setMatrix(dm); } template -void ProjectedMatrices::updateDM(const int iterative_index) +void ProjectedMatrices::updateDM() { Control& ct = *(Control::instance()); MGmol_MPI& mmpi = *(MGmol_MPI::instance()); if (ct.DMEigensolver() == DMEigensolverType::Eigensolver) - updateDMwithEigenstates(iterative_index); + updateDMwithEigenstates(); else if (ct.DMEigensolver() == DMEigensolverType::Chebyshev) - updateDMwithChebApproximation(iterative_index); + updateDMwithChebApproximation(); else if (ct.DMEigensolver() == DMEigensolverType::SP2) - updateDMwithSP2(iterative_index); + updateDMwithSP2(); else { std::cerr << "Eigensolver not available in " @@ -470,7 +467,7 @@ void ProjectedMatrices::updateDM(const int iterative_index) template void ProjectedMatrices::updateDMwithEigenstatesAndRotate( - const int iterative_index, MatrixType& zz) + MatrixType& zz) { // solves generalized eigenvalue problem // and return solution in zz @@ -479,7 +476,7 @@ void ProjectedMatrices::updateDMwithEigenstatesAndRotate( rotateAll(zz, true); - dm_->build(zz, iterative_index); + dm_->build(zz); } template @@ -614,7 +611,7 @@ void ProjectedMatrices::dressupDM() if (mmpi.instancePE0()) std::cout << "ProjectedMatrices::dressupDM()" << std::endl; #endif - dm_->dressUpS(gm_->getCholeskyL(), gm_->getAssociatedOrbitalsIndex()); + dm_->dressUpS(gm_->getCholeskyL()); } template @@ -999,8 +996,7 @@ double ProjectedMatrices::computeTraceInvSmultMatMultTheta( template double ProjectedMatrices::computeChemicalPotentialAndDMwithChebyshev( - const int order, const double emin, const double emax, - const int iterative_index) + const int order, const double emin, const double emax) { assert(emax > emin); assert(nel_ >= 0.); @@ -1175,7 +1171,7 @@ ProjectedMatrices::computeChemicalPotentialAndDMwithChebyshev( dm.gemm('N', 'N', 1., tmp, gm_->getInverse(), 0.); double orbital_occupation = mmpi.nspin() > 1 ? 1. : 2.; dm.scal(orbital_occupation); - dm_->setMatrix(dm, iterative_index); + dm_->setMatrix(dm); return mu_; } diff --git a/src/ProjectedMatrices.h b/src/ProjectedMatrices.h index f2881763..142dbf70 100644 --- a/src/ProjectedMatrices.h +++ b/src/ProjectedMatrices.h @@ -80,8 +80,8 @@ class ProjectedMatrices : public ProjectedMatricesInterface { computeChemicalPotentialAndOccupations(width_, dim_); } - double computeChemicalPotentialAndDMwithChebyshev(const int order, - const double emin, const double emax, const int iterative_index); + double computeChemicalPotentialAndDMwithChebyshev( + const int order, const double emin, const double emax); protected: // indexes corresponding to valid function in each subdomain @@ -197,12 +197,9 @@ class ProjectedMatrices : public ProjectedMatricesInterface int getDMMatrixIndex() const override { assert(dm_); - return dm_->getOrbitalsIndex(); - } - void setDMuniform(const double nel, const int orbitals_index) override - { - dm_->setUniform(nel, orbitals_index); + return dm_->getIndex(); } + void setDMuniform(const double nel) override { dm_->setUniform(nel); } int dim() const { return dim_; } void computeInvS() override; @@ -280,10 +277,9 @@ class ProjectedMatrices : public ProjectedMatricesInterface } void setDMto2InvS() override; - void buildDM(const MatrixType& z, const int orbitals_index); - void buildDM(const MatrixType& z, const std::vector&, - const int orbitals_index); - void buildDM(const std::vector&, const int orbitals_index); + void buildDM(const MatrixType& z); + void buildDM(const MatrixType& z, const std::vector&); + void buildDM(const std::vector&); double getEigSum() override; double getExpectation(const MatrixType& A); @@ -323,12 +319,11 @@ class ProjectedMatrices : public ProjectedMatricesInterface int readDM(HDFrestart& h5f_file) override; int readWFDM(HDFrestart& h5f_file); void printEigenvalues(std::ostream& os) const; - void updateDM(const int iterative_index) override; - void updateDMwithEigenstates(const int iterative_index); - void updateDMwithSP2(const int iterative_index); - void updateDMwithEigenstatesAndRotate( - const int iterative_index, MatrixType& zz); - void updateDMwithChebApproximation(const int iterative_index) override; + void updateDM() override; + void updateDMwithEigenstates(); + void updateDMwithSP2(); + void updateDMwithEigenstatesAndRotate(MatrixType& zz); + void updateDMwithChebApproximation() override; void computeChemicalPotentialAndOccupations( const double width, const int max_numst) { @@ -354,16 +349,16 @@ class ProjectedMatrices : public ProjectedMatricesInterface void resetDM() override { - dm_->setMatrix(*mat_X_old_, 0); + dm_->setMatrix(*mat_X_old_); dm_->stripS(*mat_L_old_); } - void updateDMwithRelax(const double mix, const int itindex) override + void updateDMwithRelax(const double mix) override { // cout<<"ProjectedMatrices::updateDMwithRelax()..."<mix(mix, *mat_X_old_, itindex); + dm_->mix(mix, *mat_X_old_); } SquareLocalMatrices getReplicatedDM(); @@ -414,10 +409,7 @@ class ProjectedMatrices : public ProjectedMatricesInterface pmat.gemm('n', 'n', 1.0, mat, *theta_, 0.); } MatrixType& getMatHB() { return *matHB_; } - void setDM(const MatrixType& mat, const int orbitals_index) - { - dm_->setMatrix(mat, orbitals_index); - } + void setDM(const MatrixType& mat) { dm_->setMatrix(mat); } void setEigenvalues(const std::vector& eigenvalues) { memcpy(eigenvalues_.data(), eigenvalues.data(), diff --git a/src/ProjectedMatrices2N.cc b/src/ProjectedMatrices2N.cc index b0696f94..479ee4f0 100644 --- a/src/ProjectedMatrices2N.cc +++ b/src/ProjectedMatrices2N.cc @@ -39,7 +39,7 @@ void ProjectedMatrices2N::assignBlocksH( template void ProjectedMatrices2N::iterativeUpdateDMwithEigenstates( - const double occ_width, const int iterative_index, const bool flag_reduce_T) + const double occ_width, const bool flag_reduce_T) { MGmol_MPI& mmpi = *(MGmol_MPI::instance()); @@ -64,7 +64,7 @@ void ProjectedMatrices2N::iterativeUpdateDMwithEigenstates( (*MPIdata::sout) << "MVP target with mu = " << ProjectedMatricesInterface::mu_ << " [Ry]" << std::endl; - ProjectedMatrices::buildDM(*work2N_, iterative_index); + ProjectedMatrices::buildDM(*work2N_); } template class ProjectedMatrices2N>; diff --git a/src/ProjectedMatrices2N.h b/src/ProjectedMatrices2N.h index 69ca2838..a87bf83f 100644 --- a/src/ProjectedMatrices2N.h +++ b/src/ProjectedMatrices2N.h @@ -26,8 +26,8 @@ class ProjectedMatrices2N : public ProjectedMatrices void assignBlocksH(MatrixType&, MatrixType&, MatrixType&, MatrixType&); - void iterativeUpdateDMwithEigenstates(const double occ_width, - const int iterative_index, const bool flag_reduce_T = true); + void iterativeUpdateDMwithEigenstates( + const double occ_width, const bool flag_reduce_T = true); void diagonalizeDM(std::vector& occ, MatrixType& vect) { // we are assuming Gram matrix=identity diff --git a/src/ProjectedMatricesInterface.h b/src/ProjectedMatricesInterface.h index c52595e1..53508a2e 100644 --- a/src/ProjectedMatricesInterface.h +++ b/src/ProjectedMatricesInterface.h @@ -199,10 +199,10 @@ class ProjectedMatricesInterface : public ChebyshevApproximationFunction virtual double computeEntropy() = 0; virtual double getEigSum() = 0; - virtual void updateTheta() = 0; - virtual void computeInvB() = 0; - virtual void printGramMM(std::ofstream& tfile) = 0; - virtual void setDMuniform(const double nel, const int orbitals_index) = 0; + virtual void updateTheta() = 0; + virtual void computeInvB() = 0; + virtual void printGramMM(std::ofstream& tfile) = 0; + virtual void setDMuniform(const double nel) = 0; virtual double dotProductWithInvS( const SquareLocalMatrices& ss) @@ -258,10 +258,9 @@ class ProjectedMatricesInterface : public ChebyshevApproximationFunction } virtual void saveDM() { exitWithErrorMessage("saveDM"); } virtual void resetDM() { exitWithErrorMessage("resetDM"); } - virtual void updateDMwithRelax(const double mix, const int itindex) + virtual void updateDMwithRelax(const double mix) { (void)mix; - (void)itindex; exitWithErrorMessage("updateDMwithRelax"); } @@ -297,18 +296,11 @@ class ProjectedMatricesInterface : public ChebyshevApproximationFunction return 0; } - virtual void updateDMwithChebApproximation(const int iterative_index) + virtual void updateDMwithChebApproximation() { - (void)iterative_index; - exitWithErrorMessage("updateDMwithChebApproximation"); } - virtual void updateDM(const int iterative_index) - { - (void)iterative_index; - - exitWithErrorMessage("updateDM"); - } + virtual void updateDM() { exitWithErrorMessage("updateDM"); } virtual void setDMto2InvS() { exitWithErrorMessage("setDMto2InvS"); } virtual void initializeMatB( diff --git a/src/ProjectedMatricesSparse.h b/src/ProjectedMatricesSparse.h index fd3f6918..fdbb4f88 100644 --- a/src/ProjectedMatricesSparse.h +++ b/src/ProjectedMatricesSparse.h @@ -170,15 +170,12 @@ class ProjectedMatricesSparse : public ProjectedMatricesInterface /* scale H */ // (*matHB_).scale(vel_); } - void setDMuniform(const double nel, const int orbitals_index) override - { - dm_->setUniform(nel, orbitals_index); - } + void setDMuniform(const double nel) override { dm_->setUniform(nel); } int getDMMatrixIndex() const override { assert(dm_ != nullptr); - return dm_->getOrbitalsIndex(); + return dm_->getIndex(); } int getGramMatrixIndex() const { @@ -253,7 +250,7 @@ class ProjectedMatricesSparse : public ProjectedMatricesInterface void setDMto2InvS() override { assert(invS_ != nullptr); - dm_->setto2InvS(invS_->getInvS(), invS_->getGramMatrixOrbitalsIndex()); + dm_->setto2InvS(invS_->getInvS()); } void computeInvS() override diff --git a/src/md.cc b/src/md.cc index e8d06133..ac5aa8c1 100644 --- a/src/md.cc +++ b/src/md.cc @@ -420,7 +420,7 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) std::shared_ptr projmatrices = getProjectedMatrices(); - projmatrices->setDMuniform(ct.getNelSpin(), 0); + projmatrices->setDMuniform(ct.getNelSpin()); projmatrices->printDM(os_); } #endif diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index f55ce699..38cdc4c1 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -558,7 +558,7 @@ add_test(NAME ReplicatedSP2 ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) add_test(NAME testMD_D72 COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/MD_D72/test.py - ${MPIEXEC} --oversubscribe ${MPIEXEC_NUMPROC_FLAG} 5 ${MPIEXEC_PREFLAGS} + ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 5 ${MPIEXEC_PREFLAGS} ${CMAKE_CURRENT_BINARY_DIR}/../src/mgmol-opt ${CMAKE_CURRENT_SOURCE_DIR}/MD_D72/mgmol_quench.cfg ${CMAKE_CURRENT_SOURCE_DIR}/MD_D72/mgmol_md.cfg @@ -612,7 +612,7 @@ add_test(NAME ChebyshevMVP if(NOT ${MGMOL_WITH_MAGMA}) add_test(NAME testShortSighted COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/ShortSighted/test.py - ${MPIEXEC} --oversubscribe ${MPIEXEC_NUMPROC_FLAG} 5 ${MPIEXEC_PREFLAGS} + ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 5 ${MPIEXEC_PREFLAGS} ${CMAKE_CURRENT_BINARY_DIR}/../src/mgmol-opt ${CMAKE_CURRENT_SOURCE_DIR}/ShortSighted/quench.cfg ${CMAKE_CURRENT_SOURCE_DIR}/ShortSighted/md.cfg diff --git a/tests/DMandEnergyAndForces/testDMandEnergyAndForces.cc b/tests/DMandEnergyAndForces/testDMandEnergyAndForces.cc index 1cf1cf88..a5a086fe 100644 --- a/tests/DMandEnergyAndForces/testDMandEnergyAndForces.cc +++ b/tests/DMandEnergyAndForces/testDMandEnergyAndForces.cc @@ -172,7 +172,7 @@ int main(int argc, char** argv) // // reset initial DM to test iterative solve for it - projmatrices->setDMuniform(ct.getNelSpin(), 0); + projmatrices->setDMuniform(ct.getNelSpin()); ct.dm_inner_steps = 50; eks = mgmol->evaluateDMandEnergyAndForces( &orbitals, positions, anumbers, forces); diff --git a/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc b/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc index 83bfe8a6..9cb58f33 100644 --- a/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc +++ b/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc @@ -166,7 +166,7 @@ int main(int argc, char** argv) orbitals.setIterativeIndex(10); // set initial DM with uniform occupations - projmatrices->setDMuniform(ct.getNelSpin(), 0); + projmatrices->setDMuniform(ct.getNelSpin()); projmatrices->printDM(std::cout); // diff --git a/tests/RestartEnergyAndForces/testRestartEnergyAndForces.cc b/tests/RestartEnergyAndForces/testRestartEnergyAndForces.cc index b1fb0094..058b0a8e 100644 --- a/tests/RestartEnergyAndForces/testRestartEnergyAndForces.cc +++ b/tests/RestartEnergyAndForces/testRestartEnergyAndForces.cc @@ -162,7 +162,7 @@ int main(int argc, char** argv) orbitals.setIterativeIndex(1); // set initial DM with uniform occupations - projmatrices->setDMuniform(ct.getNelSpin(), 0); + projmatrices->setDMuniform(ct.getNelSpin()); projmatrices->printDM(std::cout); // swap H and O to make sure order of atoms in list does not matter diff --git a/tests/testDensityMatrix.cc b/tests/testDensityMatrix.cc index 8a3649c1..6a65c995 100644 --- a/tests/testDensityMatrix.cc +++ b/tests/testDensityMatrix.cc @@ -95,10 +95,10 @@ TEST_CASE( // setup density matrix DensityMatrix dm(n); - dm.setMatrix(matK, 0); + dm.setMatrix(matK); dm.stripS(ls); - dm.dressUpS(ls, 1); + dm.dressUpS(ls); const MatrixType& newM = dm.getMatrix(); if (myrank == 0) std::cout << "new M" << std::endl;