From 864695710bd42b11504962da6633819387c508bb Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Thu, 2 May 2024 10:07:58 -0400 Subject: [PATCH 01/11] Fix build with HDF5P on (#231) --- CMakeLists.txt | 5 ++++- scripts/build_condo-mod.sh | 1 + src/Ions.cc | 26 +++++++++++++------------- 3 files changed, 18 insertions(+), 14 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 177d29bf..aee0ddb5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -62,7 +62,6 @@ else() #search for HDF5 if(MGMOL_USE_HDF5P) message(STATUS "Use HDF5 parallel capability") set(HDF5_PREFER_PARALLEL True) - add_definitions(-DMGMOL_USE_HDF5P) endif() message(STATUS "HDF5_ROOT: ${HDF5_ROOT}") find_package(HDF5 REQUIRED COMPONENTS C HL) @@ -73,6 +72,10 @@ else() #search for HDF5 message(FATAL_ERROR "Required HDF5 package not found.") endif (${HDF5_FOUND}) endif() +if(MGMOL_USE_HDF5P) + add_definitions(-DMGMOL_USE_HDF5P) +endif() + set(MGMOL_WITH_LIBXC FALSE CACHE BOOL "Compile with LIBXC") if(${MGMOL_WITH_LIBXC}) diff --git a/scripts/build_condo-mod.sh b/scripts/build_condo-mod.sh index f7be0c3b..d9abe034 100755 --- a/scripts/build_condo-mod.sh +++ b/scripts/build_condo-mod.sh @@ -21,6 +21,7 @@ cmake -DCMAKE_INSTALL_PREFIX=${INSTALL_DIR} \ -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" \ diff --git a/src/Ions.cc b/src/Ions.cc index ca5120ab..e2cc9916 100644 --- a/src/Ions.cc +++ b/src/Ions.cc @@ -47,9 +47,10 @@ double Ions::max_Vl_radius_ = -1.; double Ions::max_Vnl_radius_ = -1.; template -void writeData2d(hid_t file_id, std::string datasetname, std::vector& data, - const int n, T element) +void writeData2d(HDFrestart& h5f_file, std::string datasetname, + std::vector& data, const size_t n, T element) { + hid_t file_id = h5f_file.file_id(); #ifdef MGMOL_USE_HDF5P if (h5f_file.useHdf5p()) { @@ -67,7 +68,7 @@ void writeData2d(hid_t file_id, std::string datasetname, std::vector& data, else #endif { - size_t dims[2] = { data.size()/n, n }; + size_t dims[2] = { data.size() / n, n }; mgmol_tools::write2d(file_id, datasetname, data, dims); } } @@ -753,7 +754,7 @@ void Ions::writeAtomicNumbers(HDFrestart& h5f_file) if (file_id >= 0) { std::string datasetname("/Atomic_numbers"); - writeData2d(file_id, datasetname, data, 1, -1); + writeData2d(h5f_file, datasetname, data, 1, -1); } } @@ -788,12 +789,11 @@ void Ions::writeAtomNames(HDFrestart& h5f_file) // write data hid_t file_id = h5f_file.file_id(); - if (file_id >= 0) { std::string datasetname("/Atomic_names"); std::string empty; - writeData2d(file_id, datasetname, data, 1, empty); + writeData2d(h5f_file, datasetname, data, 1, empty); } } @@ -864,7 +864,7 @@ void Ions::writeLockedAtomNames(HDFrestart& h5f_file) { std::string datasetname("/LockedAtomsNames"); std::string empty; - writeData2d(file_id, datasetname, data, 1, empty); + writeData2d(h5f_file, datasetname, data, 1, empty); } } @@ -900,7 +900,7 @@ void Ions::writeAtomicIDs(HDFrestart& h5f_file) if (file_id >= 0) { std::string datasetname("/Atomic_IDs"); - writeData2d(file_id, datasetname, data, 1, -1); + writeData2d(h5f_file, datasetname, data, 1, -1); } } @@ -937,7 +937,7 @@ void Ions::writeAtomicNLprojIDs(HDFrestart& h5f_file) if (file_id >= 0) { std::string datasetname("/AtomicNLproj_IDs"); - writeData2d(file_id, datasetname, data, 1, -1); + writeData2d(h5f_file, datasetname, data, 1, -1); } } @@ -975,7 +975,7 @@ void Ions::writePositions(HDFrestart& h5f_file) if (file_id >= 0) { std::string datasetname("/Ionic_positions"); - writeData2d(file_id, datasetname, data, 3, 1.e32); + writeData2d(h5f_file, datasetname, data, 3, 1.e32); } } @@ -1138,7 +1138,7 @@ void Ions::writeVelocities(HDFrestart& h5f_file) if (file_id >= 0) { std::string datasetname("/Ionic_velocities"); - writeData2d(file_id, datasetname, data, 3, 1.e32); + writeData2d(h5f_file, datasetname, data, 3, 1.e32); } } @@ -1184,7 +1184,7 @@ void Ions::writeRandomStates(HDFrestart& h5f_file) if (file_id >= 0) { std::string datasetname("/Ionic_RandomStates"); - writeData2d(file_id, datasetname, data, 3, (unsigned short)0); + writeData2d(h5f_file, datasetname, data, 3, (unsigned short)0); } } @@ -1343,7 +1343,7 @@ void Ions::writeForces(HDFrestart& h5f_file) if (file_id >= 0) { std::string datasetname("/Ionic_forces"); - writeData2d(file_id, datasetname, data, 3, 1.e32); + writeData2d(h5f_file, datasetname, data, 3, 1.e32); } } From d46c35b57cf2658455f4296f8cf91d99444f071b Mon Sep 17 00:00:00 2001 From: "Kevin\" Seung Whan Chung" Date: Thu, 2 May 2024 11:43:33 -0700 Subject: [PATCH 02/11] LoadOrbitalsFromRestartFile and miscellaneous changes (#226) * Signal.h -> mgmol_Signal.h * Carbyne example files. * fixed FindSCALAPACK.cmake to search for default paths * MGmol::loadOrbitalFromRestartFile * nullptr initialization --- cmake_modules/FindSCALAPACK.cmake | 21 ++++++- examples/Carbyne/carbyne.cfg | 26 ++++++++ examples/Carbyne/carbyne.in | 14 +++++ src/MGmol.cc | 2 +- src/MGmol.h | 4 +- src/OrbitalsExtrapolation.h | 2 +- src/md.cc | 86 +++++++++++++++++++++++++- src/tools/Timeout.h | 2 +- src/tools/{Signal.h => mgmol_Signal.h} | 2 +- 9 files changed, 150 insertions(+), 9 deletions(-) create mode 100644 examples/Carbyne/carbyne.cfg create mode 100644 examples/Carbyne/carbyne.in rename src/tools/{Signal.h => mgmol_Signal.h} (98%) diff --git a/cmake_modules/FindSCALAPACK.cmake b/cmake_modules/FindSCALAPACK.cmake index dde12c26..74e58e7d 100644 --- a/cmake_modules/FindSCALAPACK.cmake +++ b/cmake_modules/FindSCALAPACK.cmake @@ -7,7 +7,7 @@ if(DEFINED ENV{SCALAPACK_ROOT}) endif(DEFINED ENV{SCALAPACK_ROOT}) if(SCALAPACK_ROOT) - set(_SCALAPACK_SEARCH_DIR ${SCALAPACK_ROOT}) + set(_SCALAPACK_SEARCH_DIR ${SCALAPACK_ROOT} ${SCALAPACK_ROOT}/lib/intel64) list(APPEND _SCALAPACK_SEARCHES ${_SCALAPACK_SEARCH_DIR}) endif(SCALAPACK_ROOT) @@ -29,13 +29,28 @@ if(NOT SCALAPACK_LIBRARY) endforeach() endif() -unset(SCALAPACK_NAMES) - mark_as_advanced(SCALAPACK_LIBRARY SCALAPACK_INCLUDE_DIR) include(FindPackageHandleStandardArgs) FIND_PACKAGE_HANDLE_STANDARD_ARGS(SCALAPACK REQUIRED_VARS SCALAPACK_LIBRARY) +# Search for some default library paths +if (NOT SCALAPACK_FOUND) + find_library(SCALAPACK_LIBRARY + NAMES ${SCALAPACK_NAMES} + PATHS /usr/lib64 /usr/lib /usr/local/lib64 /usr/local/lib + /opt/local/lib /opt/sw/lib /sw/lib + ENV LD_LIBRARY_PATH + ENV DYLD_FALLBACK_LIBRARY_PATH + ENV DYLD_LIBRARY_PATH + ENV SCALAPACKDIR + ENV BLACSDIR) + + FIND_PACKAGE_HANDLE_STANDARD_ARGS(SCALAPACK REQUIRED_VARS SCALAPACK_LIBRARY) +endif() + +unset(SCALAPACK_NAMES) + if(SCALAPACK_FOUND) # Only Intel's scalapack requires an include directory if(SCALAPACK_INCLUDE_DIR) diff --git a/examples/Carbyne/carbyne.cfg b/examples/Carbyne/carbyne.cfg new file mode 100644 index 00000000..d153e6f6 --- /dev/null +++ b/examples/Carbyne/carbyne.cfg @@ -0,0 +1,26 @@ +verbosity=2 +xcFunctional=PBE +FDtype=4th +[Mesh] +nx= 96 +ny= 96 +nz= 192 +[Domain] +ox= -10. +oy= -10. +oz= -20. +lx= 20. +ly= 20. +lz= 40. +[Potentials] +pseudopotential=pseudo.H_ONCV_PBE_SG15 +pseudopotential=pseudo.C_ONCV_PBE_SG15 +[Run] +type=QUENCH +[Quench] +max_steps=200 +atol=1.e-8 +[Orbitals] +initial_type=Fourier +[Restart] +output_level=3 diff --git a/examples/Carbyne/carbyne.in b/examples/Carbyne/carbyne.in new file mode 100644 index 00000000..9f4e975a --- /dev/null +++ b/examples/Carbyne/carbyne.in @@ -0,0 +1,14 @@ +H00 1 -0.0000 -0.0000 15.2674 +C01 2 -0.0000 0.0000 13.2519 +C02 2 -0.0000 0.0000 10.9495 +C03 2 -0.0000 -0.0000 8.4221 +C04 2 0.0000 0.0000 6.0897 +C05 2 -0.0000 0.0000 3.5892 +C06 2 -0.0000 -0.0000 1.2470 +C07 2 0.0000 -0.0000 -1.2469 +C08 2 0.0000 -0.0000 -3.5891 +C09 2 -0.0000 -0.0000 -6.0897 +C10 2 -0.0000 0.0000 -8.4221 +C11 2 0.0000 -0.0000 -10.9495 +C12 2 0.0000 0.0000 -13.2520 +H13 1 0.0000 0.0000 -15.2675 diff --git a/src/MGmol.cc b/src/MGmol.cc index b43c70e3..0af62ebc 100644 --- a/src/MGmol.cc +++ b/src/MGmol.cc @@ -106,7 +106,7 @@ extern Timer ions_setupInteractingIons_tm; extern Timer ions_setup_tm; extern Timer updateCenters_tm; -#include "Signal.h" +#include "mgmol_Signal.h" std::set Signal::recv_; template diff --git a/src/MGmol.h b/src/MGmol.h index 23b507c4..e63e7069 100644 --- a/src/MGmol.h +++ b/src/MGmol.h @@ -102,7 +102,7 @@ class MGmol : public MGmolInterface double total_energy_; ConstraintSet* constraints_; - OrbitalsExtrapolation* orbitals_extrapol_; + OrbitalsExtrapolation* orbitals_extrapol_ = nullptr; float md_time_; int md_iteration_; @@ -301,6 +301,8 @@ class MGmol : public MGmolInterface { forces_->force(orbitals, ions); } + + OrbitalsType* loadOrbitalFromRestartFile(const std::string filename); }; // Instantiate static variables here to avoid clang warnings template diff --git a/src/OrbitalsExtrapolation.h b/src/OrbitalsExtrapolation.h index 9f2348f7..47c1a21b 100644 --- a/src/OrbitalsExtrapolation.h +++ b/src/OrbitalsExtrapolation.h @@ -44,7 +44,7 @@ class OrbitalsExtrapolation virtual short getNumOrbitalExtrapolations() { return 0; } protected: - OrbitalsType* orbitals_minus1_; + OrbitalsType* orbitals_minus1_ = nullptr; }; #endif diff --git a/src/md.cc b/src/md.cc index 1adab46b..ec288839 100644 --- a/src/md.cc +++ b/src/md.cc @@ -31,7 +31,7 @@ #include "ProjectedMatricesMehrstellen.h" #include "ProjectedMatricesSparse.h" #include "Rho.h" -#include "Signal.h" +#include "mgmol_Signal.h" #include "SpreadsAndCenters.h" #include "tools.h" @@ -674,5 +674,89 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) delete orbitals_extrapol_; } +template +OrbitalsType* MGmol::loadOrbitalFromRestartFile(const std::string filename) +{ + MGmol_MPI& mmpi(*(MGmol_MPI::instance())); + Control& ct = *(Control::instance()); + Mesh* mymesh = Mesh::instance(); + const pb::PEenv& myPEenv = mymesh->peenv(); + + /* For now, we only consider double-precision hdf5 I/O. */ + assert(ct.restart_info > 3); + assert((ct.AtomsDynamic() == AtomsDynamicType::MD) || (ct.AtomsDynamic() == AtomsDynamicType::Quench)); + + HDFrestart h5file(filename, myPEenv, ct.out_restart_file_type); + int ierr; + + OrbitalsType *restart_orbitals = nullptr; + + /* This corresponds to MGmol::initial */ + { + // Copy from current orbital, instead of constructing brand-new one + restart_orbitals = new OrbitalsType("ForLoading", *current_orbitals_, false); + + /* This corresponds to MGmol::read_restart_data */ + { + ierr = restart_orbitals->read_func_hdf5(h5file); + } // read_restart_data + } // initial() + + /* This corresponds to MGmol::md */ + { + int flag_extrapolated_data = 0; + if (onpe0) + { + flag_extrapolated_data + = h5file.dset_exists("ExtrapolatedFunction0000"); + if (flag_extrapolated_data == 0) + flag_extrapolated_data + = h5file.dset_exists("ExtrapolatedFunction0"); + } + MPI_Bcast(&flag_extrapolated_data, 1, MPI_INT, 0, comm_); + + /* + If extrapolated function exists, + then function is set as previous orbitals, + while the extrapolated function is set as the current orbitals. + This is how the restart file is saved via dumprestartFile. + */ + if (flag_extrapolated_data) + { + orbitals_extrapol_ = OrbitalsExtrapolationFactory::create( + ct.WFExtrapolation()); + + if (onpe0) os_ << "Create new orbitals_minus1..." << std::endl; + + orbitals_extrapol_->setupPreviousOrbitals(&restart_orbitals, + proj_matrices_, lrs_, local_cluster_, currentMasks_, + corrMasks_, h5file); + } + + /* main workflow delete h5f_file_ here, meaning the loading is over. */ + } // md() + + ierr = h5file.close(); + mmpi.allreduce(&ierr, 1, MPI_MIN); + if (ierr < 0) + { + if (onpe0) + (*MPIdata::serr) + << "loadRestartFile: cannot close file..." << std::endl; + return nullptr; + } + + /* + In returning restart_orbitals, + we hope that the wavefunctions in restart_orbitals are all set. + At least the following functions should return proper data loaded from the file: + + restart_orbitals->getLocNumpt() + restart_orbitals->chromatic_number() + restart_orbitals->getPsi(idx) (for int idx) + */ + return restart_orbitals; +} + template class MGmol; template class MGmol; diff --git a/src/tools/Timeout.h b/src/tools/Timeout.h index 18940a78..25751b2a 100644 --- a/src/tools/Timeout.h +++ b/src/tools/Timeout.h @@ -17,7 +17,7 @@ #include #include "MPIdata.h" -#include "Signal.h" +#include "mgmol_Signal.h" #if PCS #include diff --git a/src/tools/Signal.h b/src/tools/mgmol_Signal.h similarity index 98% rename from src/tools/Signal.h rename to src/tools/mgmol_Signal.h index 8e8d0a21..33c47b3e 100644 --- a/src/tools/Signal.h +++ b/src/tools/mgmol_Signal.h @@ -15,7 +15,7 @@ // access the flag set, reset flags, or interrogate flags. // The Signal class can be used in an application by declaring // -// #include "Signal.h" +// #include "mgmol_Signal.h" // set Signal::recv_; // // A signal can be registered using, e.g. From e4c102bb7b4777aac3b18bc728daa4e36a07ac21 Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Wed, 22 May 2024 09:00:32 -0400 Subject: [PATCH 03/11] Clean up class Rho (#234) --- src/DavidsonSolver.cc | 2 - src/MVPSolver.cc | 2 - src/Rho.cc | 371 +++--------------------------------------- src/Rho.h | 29 +--- 4 files changed, 30 insertions(+), 374 deletions(-) diff --git a/src/DavidsonSolver.cc b/src/DavidsonSolver.cc index f21cf7fa..c53e11c0 100644 --- a/src/DavidsonSolver.cc +++ b/src/DavidsonSolver.cc @@ -598,7 +598,6 @@ int DavidsonSolver::solve( rho_->computeRho( orbitals, work_orbitals, dm11, dm12, dm21, dm22); - rho_->rescaleTotalCharge(); mgmol_strategy_->update_pot(vh_init, ions_); @@ -663,7 +662,6 @@ int DavidsonSolver::solve( rho_->computeRho( orbitals, work_orbitals, dm11, dm12, dm21, dm22); - rho_->rescaleTotalCharge(); } } // inner iterations diff --git a/src/MVPSolver.cc b/src/MVPSolver.cc index 2f8060b7..fc30d27c 100644 --- a/src/MVPSolver.cc +++ b/src/MVPSolver.cc @@ -302,7 +302,6 @@ int MVPSolver::solve(OrbitalsType& orbitals) // if( onpe0 )os_<<"Rho..."<computeRho(orbitals, target); - rho_->rescaleTotalCharge(); mgmol_strategy_->update_pot(vh_init, ions_); @@ -358,7 +357,6 @@ int MVPSolver::solve(OrbitalsType& orbitals) // if( onpe0 )os_<<"Rho..."<computeRho(orbitals, *work_); - rho_->rescaleTotalCharge(); } } // inner iterations diff --git a/src/Rho.cc b/src/Rho.cc index 823e4478..b456f6c1 100644 --- a/src/Rho.cc +++ b/src/Rho.cc @@ -28,8 +28,6 @@ template Timer Rho::compute_tm_("Rho::compute"); template Timer Rho::compute_blas_tm_("Rho::compute_usingBlas"); -// template -// Timer Rho::compute_offdiag_tm_("Rho::compute_offdiag"); #ifdef HAVE_MAGMA template @@ -46,10 +44,6 @@ Rho::Rho() myspin_ = mmpi.myspin(); nspin_ = mmpi.nspin(); - // default values for block sizes - // block_functions_ = 8; - // block_space_ = 256; - Mesh* mymesh = Mesh::instance(); const pb::Grid& mygrid = mymesh->grid(); np_ = mygrid.size(); @@ -77,31 +71,6 @@ void Rho::setup(const OrthoType orbitals_type, orbitals_indexes_ = orbitals_indexes; } -template -void Rho::extrapolate() -{ - double minus = -1; - double two = 2.; - if (rho_minus1_.empty()) - { - rho_minus1_.resize(nspin_); - rho_minus1_[myspin_].resize(np_); - memcpy(&rho_minus1_[myspin_][0], &rho_[myspin_][0], - np_ * sizeof(RHODTYPE)); - return; - } - RHODTYPE* tmp = new RHODTYPE[np_]; - memcpy(tmp, &rho_[myspin_][0], np_ * sizeof(RHODTYPE)); - - LinearAlgebraUtils::MPscal(np_, two, &rho_[myspin_][0]); - LinearAlgebraUtils::MPaxpy( - np_, minus, &rho_minus1_[myspin_][0], &rho_[myspin_][0]); - - memcpy(&rho_minus1_[myspin_][0], tmp, np_ * sizeof(RHODTYPE)); - - delete[] tmp; -} - template void Rho::axpyRhoc(const double alpha, RHODTYPE* rhoc) { @@ -144,8 +113,6 @@ void Rho::update(OrbitalsType& current_orbitals) computeRho(current_orbitals); - rescaleTotalCharge(); - assert(iterative_index_ >= 0); update_tm_.stop(); @@ -234,68 +201,6 @@ void Rho::rescaleTotalCharge() #endif } -// template -// int Rho::setupSubdomainData(const int iloc, -// const vector& vorbitals, -// const ProjectedMatricesInterface* const projmatrices, -// vector& melements, vector>& vmpsi) -//{ -// // printWithTimeStamp("Rho::setupSubdomainData()...",cout); -// -// const short norb = (short)vorbitals.size(); -// vmpsi.resize(norb); -// -// const vector& loc_indexes(orbitals_indexes_[iloc]); -// const int n_colors = (int)loc_indexes.size(); -// -// if (n_colors == 0) return 0; -// -// vector mycolors; -// mycolors.reserve(n_colors); -// for (int icolor = 0; icolor < n_colors; icolor++) -// if (loc_indexes[icolor] != -1) mycolors.push_back(icolor); -// const int nmycolors = (int)mycolors.size(); -// -// for (short j = 0; j < norb; j++) -// vmpsi[j].resize(nmycolors); -// -// short j = 0; -// for (typename vector::const_iterator it = vorbitals.begin(); -// it != vorbitals.end(); ++it) -// { -// for (int color = 0; color < nmycolors; color++) -// { -// vmpsi[j][color] = (*it)->getPsi(mycolors[color], iloc); -// assert(vmpsi[j][color] != NULL); -// } -// j++; -// } -// -// SquareLocalMatrices& -// localX(projmatrices->getLocalX()); const MATDTYPE* const localX_iloc = -// localX.getSubMatrix(iloc); melements.clear(); melements.resize(nmycolors * -// nmycolors); -// -// for (int i = 0; i < nmycolors; i++) -// { -// const int icolor = mycolors[i]; -// if (norb == 1) -// { -// melements[i * nmycolors + i] -// = localX_iloc[icolor + n_colors * icolor]; -// } -// const int jmax = (norb == 1) ? i : nmycolors; -// for (int j = 0; j < jmax; j++) -// { -// int jcolor = mycolors[j]; -// melements[j * nmycolors + i] -// = localX_iloc[icolor + n_colors * jcolor]; -// } -// } -// -// return nmycolors; -//} - template void Rho::accumulateCharge(const double alpha, const short ix_max, const ORBDTYPE* const psii, const ORBDTYPE* const psij, @@ -305,235 +210,6 @@ void Rho::accumulateCharge(const double alpha, const short ix_max, plrho[ix] += (RHODTYPE)(alpha * (double)psii[ix] * (double)psij[ix]); } -// template -// void Rho::computeRhoSubdomain( -// const int iloc_init, const int iloc_end, const OrbitalsType& orbitals) -//{ -// assert(orbitals_type_ == OrthoType::Eigenfunctions -// || orbitals_type_ == OrthoType::Nonorthogonal); -// -// compute_tm_.start(); -// -// Mesh* mymesh = Mesh::instance(); -// -// const int loc_numpt = mymesh->locNumpt(); -// -// RHODTYPE* const prho = &rho_[myspin_][0]; -// -// vector vorbitals; -// vorbitals.push_back(&orbitals); -// -// const double nondiagfactor = 2.; -// -// for (int iloc = iloc_init; iloc < iloc_end; iloc++) -// { -// const int istart = iloc * loc_numpt; -// -// RHODTYPE* const lrho = &prho[istart]; -// -// vector melements; -// vector> vmpsi; -// -// const int nmycolors = setupSubdomainData( -// iloc, vorbitals, orbitals.projMatrices(), melements, vmpsi); -// assert(vmpsi.size() == 1); -// vector mpsi = vmpsi[0]; -// -// const int nblocks_color = nmycolors / block_functions_; -// const int max_icolor = (nmycolors % block_functions_ == 0) -// ? nmycolors -// : nblocks_color * block_functions_; -// const int missed_rows = nmycolors - max_icolor; -// //(*MPIdata::sout)<<"max_icolor="<::computeRhoSubdomain: -// //missed_rows="<::computeRhoSubdomainOffDiagBlock()...",cout); -// -// Mesh* mymesh = Mesh::instance(); -// -// const int loc_numpt = mymesh->locNumpt(); -// -// RHODTYPE* const prho = &rho_[myspin_][0]; -// -// for (int iloc = iloc_init; iloc < iloc_end; iloc++) -// { -// const int istart = iloc * loc_numpt; -// -// RHODTYPE* const lrho = &prho[istart]; -// -// vector melements; -// vector> vmpsi; -// -// const int nmycolors = setupSubdomainData( -// iloc, vorbitals, projmatrices, melements, vmpsi); -// const int nblocks_color = nmycolors / block_functions_; -// const int max_icolor = (nmycolors % block_functions_ == 0) -// ? nmycolors -// : nblocks_color * block_functions_; -// const int missed_rows = nmycolors - max_icolor; -// //(*MPIdata::sout)<<"max_icolor="<::computeRhoSubdomainOffDiagBlock: -// //missed_rows="<& occ); - // void computeRhoSubdomainOffDiagBlock(const int iloc_init, - // const int iloc_end, - // const std::vector& vorbitals, - // const ProjectedMatricesInterface* const); void computeRhoSubdomainUsingBlas3(const int iloc_init, const int iloc_end, const OrbitalsType& orbitals1, const OrbitalsType& orbitals2); void computeRhoSubdomainUsingBlas3( @@ -64,11 +53,6 @@ class Rho void accumulateCharge(const double alpha, const short ix_max, const ORBDTYPE* const psii, const ORBDTYPE* const psij, RHODTYPE* const plrho); - // int setupSubdomainData(const int iloc, - // const std::vector& vorbitals, - // const ProjectedMatricesInterface* const projmatrices, - // std::vector& melements, - // std::vector>& mpsi); void computeRho(OrbitalsType& orbitals); void computeRho( @@ -76,15 +60,17 @@ class Rho void gatherSpin(); + void resetValues(); + + void rescaleTotalCharge(); + public: // electronic density on grid std::vector> rho_; - std::vector> rho_minus1_; Rho(); ~Rho(){}; - void rescaleTotalCharge(); void setup( const OrthoType orbitals_type, const std::vector>&); void setVerbosityLevel(const int vlevel) { verbosity_level_ = vlevel; } @@ -105,18 +91,11 @@ class Rho void initUniform(); int getIterativeIndex() const { return iterative_index_; } int readRestart(HDFrestart& file); - void extrapolate(); void axpyRhoc(const double alpha, RHODTYPE* rhoc); template double dotWithRho(const T2* const func) const; - // void setupBlockSizes(const int block_functions, const int block_space) - //{ - // block_functions_ = block_functions; - // block_space_ = block_space; - //} - static void printTimers(std::ostream& os); }; From 0f093c5a188d7adc7656206e026fd6d9137d59c0 Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Wed, 22 May 2024 09:00:58 -0400 Subject: [PATCH 04/11] Clean up class MVPSolver (#236) --- src/MVPSolver.cc | 85 ++++++++++++++++++------------------------------ src/MVPSolver.h | 15 +++------ 2 files changed, 36 insertions(+), 64 deletions(-) diff --git a/src/MVPSolver.cc b/src/MVPSolver.cc index fc30d27c..8ac6a0be 100644 --- a/src/MVPSolver.cc +++ b/src/MVPSolver.cc @@ -52,18 +52,22 @@ MVPSolver::MVPSolver(MPI_Comm comm, std::ostream& os, os_(os), n_inner_steps_(n_inner_steps), use_old_dm_(use_old_dm), - ions_(ions) + ions_(ions), + numst_(numst) { - history_length_ = 2; - eks_history_.resize(history_length_, 100000.); + Control& ct = *(Control::instance()); + if (onpe0 && ct.verbose > 0) + { + os_ << "MVPSolver..." << std::endl; + if (use_old_dm_) os_ << "MVPSolver uses old DM..." << std::endl; + } rho_ = rho; energy_ = energy; electrostat_ = electrostat; mgmol_strategy_ = mgmol_strategy; - numst_ = numst; - work_ = new MatrixType("workMVP", numst_, numst_); + work_ = new MatrixType("workMVP", numst_, numst_); proj_mat_work_ = new ProjectedMatrices(numst_, false, kbT); proj_mat_work_->setup(global_indexes); @@ -123,9 +127,9 @@ void MVPSolver::buildTarget_MVP( { target_tm_.start(); - if (onpe0) std::cout << "buildTarget_MVP()..." << std::endl; - Control& ct = *(Control::instance()); + if (onpe0 && ct.verbose > 1) + std::cout << "buildTarget_MVP()..." << std::endl; proj_mat_work_->assignH(h11); @@ -135,10 +139,8 @@ void MVPSolver::buildTarget_MVP( // // compute target density matrix, with occupations>0 in numst only // - // if( onpe0 )os_<<"Compute XN..."<setHB2H(); - // if( onpe0 )os_<<"Compute DM..."<updateDM(orbitals_index); target = proj_mat_work_->dm(); @@ -168,18 +170,12 @@ int MVPSolver::solve(OrbitalsType& orbitals) os_ << "---------------------------------------------------------------" "-" << std::endl; - os_ << "Update DM functions using MVP Solver..." << std::endl; + os_ << "Update DM using MVP Solver..." << std::endl; os_ << "---------------------------------------------------------------" "-" << std::endl; } - // step for numerical derivative - - // double nel=orbitals.projMatrices()->getNel(); - // if( onpe0 ) - // os_<<"NEW algorithm: Number of electrons at start = "<::solve(OrbitalsType& orbitals) energy_->saveVofRho(); } - ProjectedMatricesInterface* current_proj_mat - = (inner_it == 0) ? orbitals.getProjMatrices() : proj_mat_work_; + ProjectedMatrices* current_proj_mat + = (inner_it == 0) + ? dynamic_cast*>( + orbitals.getProjMatrices()) + : proj_mat_work_; - double ts0; - double e0 = 0.; const int printE = (ct.verbose > 1) ? 1 : 0; // compute h11 for the current potential by adding local part to @@ -239,49 +236,32 @@ int MVPSolver::solve(OrbitalsType& orbitals) MatrixType h11(h11_nl); mgmol_strategy_->addHlocal2matrix(orbitals, orbitals, h11); + current_proj_mat->assignH(h11); + current_proj_mat->setHB2H(); + if (inner_it == 0) { - // orbitals are new, so a few things need to recomputed - ProjectedMatrices* projmatrices - = dynamic_cast*>( - orbitals.getProjMatrices()); - - projmatrices->assignH(h11); - projmatrices->setHB2H(); - if (use_old_dm_) { - dmInit = projmatrices->dm(); + dmInit = current_proj_mat->dm(); } - - ts0 = evalEntropyMVP(projmatrices, true, os_); - e0 = energy_->evaluateTotal( - ts0, projmatrices, orbitals, printE, os_); } else { - dmInit = proj_mat_work_->dm(); - - proj_mat_work_->assignH(h11); - proj_mat_work_->setHB2H(); - - ts0 = evalEntropyMVP(proj_mat_work_, true, os_); - e0 = energy_->evaluateTotal( - ts0, proj_mat_work_, orbitals, printE, os_); } - // N x N target... - // proj_mat_work_->setHiterativeIndex(orbitals.getIterativeIndex(), - // pot.getIterativeIndex()); + const double ts0 = evalEntropyMVP(current_proj_mat, true, os_); + const double e0 = energy_->evaluateTotal( + ts0, current_proj_mat, orbitals, printE, os_); MatrixType target("target", numst_, numst_); - MatrixType delta_dm("delta_dm", numst_, numst_); buildTarget_MVP(h11, s11, target); if (use_old_dm_ || inner_it > 0) { + MatrixType delta_dm("delta_dm", numst_, numst_); delta_dm = target; delta_dm -= dmInit; @@ -291,16 +271,15 @@ int MVPSolver::solve(OrbitalsType& orbitals) // evaluate free energy at beta=1 // if (onpe0 && ct.verbose > 2) - std::cout << "Target energy..." << std::endl; + std::cout << "MVP --- Target energy..." << std::endl; proj_mat_work_->setDM(target, orbitals.getIterativeIndex()); proj_mat_work_->computeOccupationsFromDM(); if (ct.verbose > 2) current_proj_mat->printOccupations(os_); - double nel = proj_mat_work_->getNel(); + const double nel = proj_mat_work_->getNel(); if (onpe0 && ct.verbose > 1) - os_ << "Number of electrons at beta=1 : " << nel + os_ << "MVP --- Number of electrons at beta=1 : " << nel << std::endl; - // if( onpe0 )os_<<"Rho..."<computeRho(orbitals, target); mgmol_strategy_->update_pot(vh_init, ions_); @@ -322,13 +301,13 @@ int MVPSolver::solve(OrbitalsType& orbitals) ts1, proj_mat_work_, orbitals, ct.verbose - 1, os_); // line minimization - double beta + const double beta = minQuadPolynomial(e0, e1, de0, (ct.verbose > 2), os_); if (onpe0 && ct.verbose > 0) { os_ << std::setprecision(12); - os_ << std::fixed << "Inner iteration " << inner_it + os_ << std::fixed << "MVP inner iteration " << inner_it << ", E0=" << e0 << ", E1=" << e1; os_ << std::scientific << ", E0'=" << de0 << " -> beta=" << beta; @@ -354,13 +333,13 @@ int MVPSolver::solve(OrbitalsType& orbitals) os_ << "Number of electrons for interpolated DM = " << pnel << std::endl; } - // if( onpe0 )os_<<"Rho..."<computeRho(orbitals, *work_); } } // inner iterations + // set DM to latest iteration value ProjectedMatrices* projmatrices = dynamic_cast*>( orbitals.getProjMatrices()); @@ -373,11 +352,9 @@ int MVPSolver::solve(OrbitalsType& orbitals) if (onpe0 && ct.verbose > 1) { os_ << "---------------------------------------------------------------" - "-" << std::endl; os_ << "End MVP Solver..." << std::endl; os_ << "---------------------------------------------------------------" - "-" << std::endl; } solve_tm_.stop(); diff --git a/src/MVPSolver.h b/src/MVPSolver.h index a3c29d1d..8e7a7c70 100644 --- a/src/MVPSolver.h +++ b/src/MVPSolver.h @@ -24,27 +24,22 @@ template class MVPSolver { private: - MPI_Comm comm_; + const MPI_Comm comm_; std::ostream& os_; - short n_inner_steps_; + const short n_inner_steps_; - bool use_old_dm_; + const bool use_old_dm_; Ions& ions_; + int numst_; + Rho* rho_; Energy* energy_; Electrostatic* electrostat_; - int history_length_; - std::vector eks_history_; - MGmol* mgmol_strategy_; - double de_old_; - double de_; - - int numst_; MatrixType* work_; ProjectedMatrices* proj_mat_work_; From 12bbed07dba16b3a91cfb4ff58c0c406fb64e515 Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Mon, 27 May 2024 15:25:33 -0400 Subject: [PATCH 05/11] Enable orthogonalization with unoccupied states (#238) Previously, runs with unoccupied states would not be possible with orthonormalization every n steps. This fixes it. Now tested in MVP test. --- src/DFTsolver.cc | 11 +++++++++-- src/ProjectedMatrices.cc | 15 +++++++++++---- tests/MVP/mvp.cfg | 1 + 3 files changed, 21 insertions(+), 6 deletions(-) diff --git a/src/DFTsolver.cc b/src/DFTsolver.cc index 034366fd..5e512090 100644 --- a/src/DFTsolver.cc +++ b/src/DFTsolver.cc @@ -316,7 +316,7 @@ int DFTsolver::solve(OrbitalsType& orbitals, const bool ortho = (ct.getOrthoType() == OrthoType::Eigenfunctions || orthof); - if (!ortho) + if (!ortho || !ct.fullyOccupied()) { // strip dm from the overlap contribution // dm <- Ls**T * dm * Ls @@ -337,7 +337,14 @@ int DFTsolver::solve(OrbitalsType& orbitals, } else { - orbitals.orthonormalizeLoewdin(); + bool updateDM = false; + if (!ct.fullyOccupied()) + { + orbitals.computeGramAndInvS(); + dm_strategy_->dressDM(); + updateDM = true; + } + orbitals.orthonormalizeLoewdin(true, nullptr, updateDM); orbitals_stepper_->restartMixing(); } diff --git a/src/ProjectedMatrices.cc b/src/ProjectedMatrices.cc index d3707657..6d6b6339 100644 --- a/src/ProjectedMatrices.cc +++ b/src/ProjectedMatrices.cc @@ -448,12 +448,14 @@ void ProjectedMatrices::updateDM(const int iterative_index) #ifndef NDEBUG double nel = getNel(); - std::cout << "ProjectedMatrices::updateDM(), nel = " << nel - << std::endl; + if (mmpi.instancePE0()) + std::cout << "ProjectedMatrices::updateDM(), nel = " << nel + << std::endl; assert(std::isfinite(nel)); double energy = getExpectationH(); - std::cout << "ProjectedMatrices::updateDM(), energy = " - << energy << std::endl; + if (mmpi.instancePE0()) + std::cout << "ProjectedMatrices::updateDM(), energy = " + << energy << std::endl; #endif } @@ -935,6 +937,11 @@ void ProjectedMatrices::computeLoewdinTransform( if (transform_matrices) { + Control& ct = *(Control::instance()); + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + if (mmpi.instancePE0() && ct.verbose > 1) + std::cout << "Transform DM to reflect Loewdin orthonormalization" + << std::endl; assert(sqrtMat); // transform DM to reflect Loewdin orthonormalization diff --git a/tests/MVP/mvp.cfg b/tests/MVP/mvp.cfg index e34d3513..25d51ed2 100644 --- a/tests/MVP/mvp.cfg +++ b/tests/MVP/mvp.cfg @@ -20,6 +20,7 @@ type=QUENCH solver=PSD max_steps=300 atol=1.e-7 +ortho_freq=10 [Orbitals] nempty=10 initial_type=random From 86fd52c5e32711bcdcfd32ae58cca7061a365c15 Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Thu, 30 May 2024 11:15:46 -0400 Subject: [PATCH 06/11] MD with unoccupied states (#239) * Fix MD with unoccupied states * Fix interface to DM solvers to use updated orbitals --- src/DFTsolver.cc | 5 +-- src/DFTsolver.h | 8 ++--- src/DMStrategy.h | 5 +-- src/DMStrategyFactory.cc | 33 ++++++++++-------- src/DMStrategyFactory.h | 27 ++++++++------- src/EigenDMStrategy.cc | 27 +++++++-------- src/EigenDMStrategy.h | 12 +++---- src/ExtendedGridOrbitals.cc | 2 ++ src/ExtendedGridOrbitals.h | 6 +++- src/FullyOccupiedNonOrthoDMStrategy.cc | 20 +++++++++-- src/FullyOccupiedNonOrthoDMStrategy.h | 7 ++-- src/HamiltonianMVP_DMStrategy.cc | 11 +++--- src/HamiltonianMVP_DMStrategy.h | 8 ++--- src/MGmol.cc | 2 +- src/MGmol.h | 4 +-- src/MVPSolver.cc | 1 + src/MVP_DMStrategy.cc | 12 +++---- src/MVP_DMStrategy.h | 11 +++--- src/NonOrthoDMStrategy.cc | 22 ++++++------ src/NonOrthoDMStrategy.h | 13 ++++--- src/Orbitals.h | 3 ++ src/OrbitalsExtrapolationOrder2.cc | 2 +- src/PolakRibiereSolver.cc | 5 +-- src/PolakRibiereSolver.h | 8 ++--- src/ProjectedMatrices.cc | 12 +++++-- src/ProjectedMatrices.h | 8 +++-- src/md.cc | 48 ++++++++++++++------------ src/quench.cc | 2 +- 28 files changed, 182 insertions(+), 142 deletions(-) diff --git a/src/DFTsolver.cc b/src/DFTsolver.cc index 5e512090..c1484160 100644 --- a/src/DFTsolver.cc +++ b/src/DFTsolver.cc @@ -24,7 +24,8 @@ template DFTsolver::DFTsolver(Hamiltonian* hamiltonian, ProjectedMatricesInterface* proj_matrices, Energy* energy, Electrostatic* electrostat, MGmol* mgmol_strategy, Ions& ions, - Rho* rho, DMStrategy* dm_strategy, std::ostream& os) + Rho* rho, DMStrategy* dm_strategy, + std::ostream& os) : mgmol_strategy_(mgmol_strategy), hamiltonian_(hamiltonian), proj_matrices_(proj_matrices), @@ -391,7 +392,7 @@ int DFTsolver::solve(OrbitalsType& orbitals, mgmol_strategy_->updateHmatrix(orbitals, ions); // compute new density matrix - dm_strategy_->update(); + dm_strategy_->update(orbitals); incInnerIt(); diff --git a/src/DFTsolver.h b/src/DFTsolver.h index e348fee7..9ca72d17 100644 --- a/src/DFTsolver.h +++ b/src/DFTsolver.h @@ -10,6 +10,7 @@ #ifndef MGMOL_DFTSOLVER_H #define MGMOL_DFTSOLVER_H +#include "DMStrategy.h" #include "DielectricControl.h" #include "Energy.h" #include "Hamiltonian.h" @@ -23,7 +24,6 @@ class Ions; class Electrostatic; class ProjectedMatricesInterface; -class DMStrategy; template class DFTsolver @@ -40,7 +40,7 @@ class DFTsolver Electrostatic* electrostat_; Ions& ions_; Rho* rho_; - DMStrategy* dm_strategy_; + DMStrategy* dm_strategy_; OrbitalsStepper* orbitals_stepper_; @@ -72,8 +72,8 @@ class DFTsolver DFTsolver(Hamiltonian* hamiltonian, ProjectedMatricesInterface* proj_matrices, Energy* energy, Electrostatic* electrostat, MGmol* mgmol_strategy, - Ions& ions, Rho* rho, DMStrategy* dm_strategy, - std::ostream& os); + Ions& ions, Rho* rho, + DMStrategy* dm_strategy, std::ostream& os); ~DFTsolver(); diff --git a/src/DMStrategy.h b/src/DMStrategy.h index a80aae1d..69322347 100644 --- a/src/DMStrategy.h +++ b/src/DMStrategy.h @@ -10,11 +10,12 @@ #ifndef DMSTRATEGY_H #define DMSTRATEGY_H +template class DMStrategy { public: - virtual void initialize() = 0; - virtual int update() = 0; + virtual void initialize(OrbitalsType& orbitals) = 0; + virtual int update(OrbitalsType& orbitals) = 0; virtual ~DMStrategy(){}; diff --git a/src/DMStrategyFactory.cc b/src/DMStrategyFactory.cc index 3ebd3e89..6dae6870 100644 --- a/src/DMStrategyFactory.cc +++ b/src/DMStrategyFactory.cc @@ -2,7 +2,7 @@ #include "ReplicatedMatrix.h" template <> -DMStrategy* DMStrategyFactory* DMStrategyFactory>::createHamiltonianMVP_DMStrategy(MPI_Comm comm, std::ostream& os, Ions& ions, Rho* rho, @@ -13,7 +13,7 @@ DMStrategy* DMStrategyFactory* dm_strategy = new HamiltonianMVP_DMStrategy, ProjectedMatricesSparse, LocGridOrbitals>(comm, os, ions, rho, energy, electrostat, mgmol_strategy, orbitals); @@ -22,18 +22,19 @@ DMStrategy* DMStrategyFactory, - ProjectedMatrices>, - LocGridOrbitals>( - comm, os, ions, rho, energy, electrostat, mgmol_strategy, orbitals); + DMStrategy* dm_strategy + = new HamiltonianMVP_DMStrategy< + dist_matrix::DistMatrix, + ProjectedMatrices>, + LocGridOrbitals>(comm, os, ions, rho, energy, electrostat, + mgmol_strategy, orbitals); return dm_strategy; } } template <> -DMStrategy* DMStrategyFactory* DMStrategyFactory>::createHamiltonianMVP_DMStrategy(MPI_Comm comm, std::ostream& os, Ions& ions, Rho* rho, @@ -44,7 +45,7 @@ DMStrategy* DMStrategyFactory* dm_strategy = new HamiltonianMVP_DMStrategy, ProjectedMatrices>, ExtendedGridOrbitals>( @@ -55,19 +56,21 @@ DMStrategy* DMStrategyFactory -DMStrategy* DMStrategyFactory* DMStrategyFactory::createHamiltonianMVP_DMStrategy(MPI_Comm comm, std::ostream& os, Ions& ions, Rho* rho, Energy* energy, Electrostatic* electrostat, MGmol* mgmol_strategy, - ProjectedMatricesInterface* /*proj_matrices*/, - ExtendedGridOrbitals* orbitals, const bool short_sighted) + ProjectedMatricesInterface* /*proj_matrices*/, LocGridOrbitals* orbitals, + const bool short_sighted) { (void)short_sighted; - DMStrategy* dm_strategy = new HamiltonianMVP_DMStrategy, ExtendedGridOrbitals>( - comm, os, ions, rho, energy, electrostat, mgmol_strategy, orbitals); + DMStrategy* dm_strategy + = new HamiltonianMVP_DMStrategy, ExtendedGridOrbitals>(comm, os, + ions, rho, energy, electrostat, mgmol_strategy, + orbitals->getOverlappingGids()); return dm_strategy; } diff --git a/src/DMStrategyFactory.h b/src/DMStrategyFactory.h index aa8d6e27..24f4f272 100644 --- a/src/DMStrategyFactory.h +++ b/src/DMStrategyFactory.h @@ -24,20 +24,20 @@ template class DMStrategyFactory { public: - static DMStrategy* create(MPI_Comm comm, std::ostream& os, Ions& ions, - Rho* rho, Energy* energy, + static DMStrategy* create(MPI_Comm comm, std::ostream& os, + Ions& ions, Rho* rho, Energy* energy, Electrostatic* electrostat, MGmol* mgmol_strategy, ProjectedMatricesInterface* proj_matrices, OrbitalsType* orbitals) { Control& ct = *(Control::instance()); MGmol_MPI& mmpi = *(MGmol_MPI::instance()); - DMStrategy* dm_strategy = nullptr; + DMStrategy* dm_strategy = nullptr; if (ct.DM_solver() == DMNonLinearSolverType::MVP) { dm_strategy = new MVP_DMStrategy(comm, os, - ions, rho, energy, electrostat, mgmol_strategy, orbitals, - proj_matrices, ct.use_old_dm()); + ions, rho, energy, electrostat, mgmol_strategy, + orbitals->getOverlappingGids(), proj_matrices, ct.use_old_dm()); } else if (ct.DM_solver() == DMNonLinearSolverType::HMVP) { @@ -51,8 +51,8 @@ class DMStrategyFactory { if (mmpi.instancePE0()) std::cout << "Fully occupied strategy" << std::endl; - dm_strategy - = new FullyOccupiedNonOrthoDMStrategy(proj_matrices); + dm_strategy = new FullyOccupiedNonOrthoDMStrategy( + proj_matrices); } else { @@ -60,8 +60,8 @@ class DMStrategyFactory { if (mmpi.instancePE0()) std::cout << "EigenDMStrategy..." << std::endl; - dm_strategy = new EigenDMStrategy( - orbitals, proj_matrices); + dm_strategy + = new EigenDMStrategy(proj_matrices); } else { @@ -70,7 +70,7 @@ class DMStrategyFactory if (mmpi.instancePE0()) std::cout << "NonOrthoDMStrategy..." << std::endl; dm_strategy = new NonOrthoDMStrategy( - orbitals, proj_matrices, ct.dm_mix); + proj_matrices, ct.dm_mix); } } } @@ -81,11 +81,12 @@ class DMStrategyFactory } private: - static DMStrategy* createHamiltonianMVP_DMStrategy(MPI_Comm comm, - std::ostream& os, Ions& ions, Rho* rho, + static DMStrategy* createHamiltonianMVP_DMStrategy( + MPI_Comm comm, std::ostream& os, Ions& ions, Rho* rho, Energy* energy, Electrostatic* electrostat, MGmol* mgmol_strategy, - ProjectedMatricesInterface* proj_matrices, OrbitalsType*, const bool); + ProjectedMatricesInterface* proj_matrices, OrbitalsType* orbitals, + const bool); }; #endif diff --git a/src/EigenDMStrategy.cc b/src/EigenDMStrategy.cc index b6563679..8fd8443a 100644 --- a/src/EigenDMStrategy.cc +++ b/src/EigenDMStrategy.cc @@ -13,21 +13,21 @@ #include "LocGridOrbitals.h" #include "ProjectedMatrices.h" -template -EigenDMStrategy::EigenDMStrategy( - T* current_orbitals, ProjectedMatricesInterface* proj_matrices) - : current_orbitals_(current_orbitals), proj_matrices_(proj_matrices) +template +EigenDMStrategy::EigenDMStrategy( + ProjectedMatricesInterface* proj_matrices) + : proj_matrices_(proj_matrices) { } -template -void EigenDMStrategy::initialize() +template +void EigenDMStrategy::initialize(OrbitalsType& orbitals) { - update(); + update(orbitals); } -template -int EigenDMStrategy::update() +template +int EigenDMStrategy::update(OrbitalsType& orbitals) { Control& ct = *(Control::instance()); @@ -37,14 +37,13 @@ int EigenDMStrategy::update() = dynamic_cast< ProjectedMatrices>*>( proj_matrices_); - pmat->updateDMwithEigenstatesAndRotate( - current_orbitals_->getIterativeIndex(), zz); + pmat->updateDMwithEigenstatesAndRotate(orbitals.getIterativeIndex(), zz); // if( onpe0 && ct.verbose>2 ) // (*MPIdata::sout)<<"get_dm_diag: rotate orbitals "<multiply_by_matrix(zz); - current_orbitals_->setDataWithGhosts(); - current_orbitals_->trade_boundaries(); + orbitals.multiply_by_matrix(zz); + orbitals.setDataWithGhosts(); + orbitals.trade_boundaries(); return 0; } diff --git a/src/EigenDMStrategy.h b/src/EigenDMStrategy.h index f2a61cba..619456e3 100644 --- a/src/EigenDMStrategy.h +++ b/src/EigenDMStrategy.h @@ -13,19 +13,17 @@ #include "DMStrategy.h" class ProjectedMatricesInterface; -template -class EigenDMStrategy : public DMStrategy +template +class EigenDMStrategy : public DMStrategy { private: - T* current_orbitals_; ProjectedMatricesInterface* proj_matrices_; public: - EigenDMStrategy( - T* current_orbitals, ProjectedMatricesInterface* proj_matrices); + EigenDMStrategy(ProjectedMatricesInterface* proj_matrices); - void initialize() override; - int update() override; + void initialize(OrbitalsType& orbitals) override; + int update(OrbitalsType& orbitals) override; bool needH() const override { return true; } diff --git a/src/ExtendedGridOrbitals.cc b/src/ExtendedGridOrbitals.cc index a11b4230..12e5b04c 100644 --- a/src/ExtendedGridOrbitals.cc +++ b/src/ExtendedGridOrbitals.cc @@ -188,6 +188,8 @@ void ExtendedGridOrbitals::reset(MasksSet* masks, MasksSet* corrmasks, void ExtendedGridOrbitals::assign(const ExtendedGridOrbitals& orbitals) { + assert(proj_matrices_ != nullptr); + assign_tm_.start(); setIterativeIndex(orbitals); diff --git a/src/ExtendedGridOrbitals.h b/src/ExtendedGridOrbitals.h index 14753830..fb7a4b42 100644 --- a/src/ExtendedGridOrbitals.h +++ b/src/ExtendedGridOrbitals.h @@ -184,7 +184,11 @@ class ExtendedGridOrbitals : public Orbitals virtual void assign(const ExtendedGridOrbitals& orbitals); void copyDataFrom(const ExtendedGridOrbitals& src); - ProjectedMatricesInterface* getProjMatrices() { return proj_matrices_; } + ProjectedMatricesInterface* getProjMatrices() + { + assert(proj_matrices_ != nullptr); + return proj_matrices_; + } const ProjectedMatricesInterface* projMatrices() const { diff --git a/src/FullyOccupiedNonOrthoDMStrategy.cc b/src/FullyOccupiedNonOrthoDMStrategy.cc index 707514f5..8887bc69 100644 --- a/src/FullyOccupiedNonOrthoDMStrategy.cc +++ b/src/FullyOccupiedNonOrthoDMStrategy.cc @@ -8,19 +8,30 @@ // Please also read this link https://github.com/llnl/mgmol/LICENSE #include "FullyOccupiedNonOrthoDMStrategy.h" +#include "ExtendedGridOrbitals.h" +#include "LocGridOrbitals.h" #include "ProjectedMatricesInterface.h" -FullyOccupiedNonOrthoDMStrategy::FullyOccupiedNonOrthoDMStrategy( +template +FullyOccupiedNonOrthoDMStrategy::FullyOccupiedNonOrthoDMStrategy( ProjectedMatricesInterface* proj_matrices) : proj_matrices_(proj_matrices) { } -void FullyOccupiedNonOrthoDMStrategy::initialize() { update(); } +template +void FullyOccupiedNonOrthoDMStrategy::initialize( + OrbitalsType& orbitals) +{ + update(orbitals); +} -int FullyOccupiedNonOrthoDMStrategy::update() +template +int FullyOccupiedNonOrthoDMStrategy::update( + OrbitalsType& orbitals) { assert(proj_matrices_ != nullptr); + (void)orbitals; proj_matrices_->setDMto2InvS(); @@ -34,3 +45,6 @@ int FullyOccupiedNonOrthoDMStrategy::update() return 0; // success } + +template class FullyOccupiedNonOrthoDMStrategy; +template class FullyOccupiedNonOrthoDMStrategy; diff --git a/src/FullyOccupiedNonOrthoDMStrategy.h b/src/FullyOccupiedNonOrthoDMStrategy.h index 5a3e6d80..01a66830 100644 --- a/src/FullyOccupiedNonOrthoDMStrategy.h +++ b/src/FullyOccupiedNonOrthoDMStrategy.h @@ -13,15 +13,16 @@ #include "DMStrategy.h" class ProjectedMatricesInterface; -class FullyOccupiedNonOrthoDMStrategy : public DMStrategy +template +class FullyOccupiedNonOrthoDMStrategy : public DMStrategy { ProjectedMatricesInterface* proj_matrices_; public: FullyOccupiedNonOrthoDMStrategy(ProjectedMatricesInterface* proj_matrices); - void initialize() override; - int update() override; + void initialize(OrbitalsType& orbitals) override; + int update(OrbitalsType& orbitals) override; bool needH() const override { return false; } diff --git a/src/HamiltonianMVP_DMStrategy.cc b/src/HamiltonianMVP_DMStrategy.cc index 3a953e91..2b5f8913 100644 --- a/src/HamiltonianMVP_DMStrategy.cc +++ b/src/HamiltonianMVP_DMStrategy.cc @@ -24,8 +24,7 @@ HamiltonianMVP_DMStrategy* rho, Energy* energy, Electrostatic* electrostat, MGmol* mgmol_strategy, OrbitalsType* orbitals) - : orbitals_(orbitals), - comm_(comm), + : comm_(comm), os_(os), ions_(ions), rho_(rho), @@ -58,13 +57,13 @@ HamiltonianMVP_DMStrategy void HamiltonianMVP_DMStrategy::initialize() + OrbitalsType>::initialize(OrbitalsType& orbitals) { } template -int HamiltonianMVP_DMStrategy::update() +int HamiltonianMVP_DMStrategy::update( + OrbitalsType& orbitals) { assert(solver_ != nullptr); @@ -75,7 +74,7 @@ int HamiltonianMVP_DMStrategysolve(*orbitals_); + return solver_->solve(orbitals); } template diff --git a/src/HamiltonianMVP_DMStrategy.h b/src/HamiltonianMVP_DMStrategy.h index 2ddc91af..da0e58b1 100644 --- a/src/HamiltonianMVP_DMStrategy.h +++ b/src/HamiltonianMVP_DMStrategy.h @@ -22,11 +22,9 @@ template class MGmol; template -class HamiltonianMVP_DMStrategy : public DMStrategy +class HamiltonianMVP_DMStrategy : public DMStrategy { private: - OrbitalsType* orbitals_; - MPI_Comm comm_; std::ostream& os_; @@ -47,8 +45,8 @@ class HamiltonianMVP_DMStrategy : public DMStrategy ~HamiltonianMVP_DMStrategy() override; - void initialize() override; - int update() override; + void initialize(OrbitalsType& orbitals) override; + int update(OrbitalsType& orbitals) override; // H is updated with HamiltonianMVP loop, so no need to compute it outside bool needH() const override { return false; } diff --git a/src/MGmol.cc b/src/MGmol.cc index 0af62ebc..fb6c64f6 100644 --- a/src/MGmol.cc +++ b/src/MGmol.cc @@ -511,7 +511,7 @@ int MGmol::initial() // theta = invB * Hij proj_matrices_->updateThetaAndHB(); - dm_strategy_->initialize(); + dm_strategy_->initialize(*current_orbitals_); if (ct.verbose > 1) { diff --git a/src/MGmol.h b/src/MGmol.h index e63e7069..d43be062 100644 --- a/src/MGmol.h +++ b/src/MGmol.h @@ -35,13 +35,13 @@ class ProjectedMatricesInterface; class KBPsiMatrix; class KBPsiMatrixSparse; class MasksSet; -class DMStrategy; template class IonicAlgorithm; #include "AOMMprojector.h" #include "ClusterOrbitals.h" +#include "DMStrategy.h" #include "ExtendedGridOrbitals.h" #include "Forces.h" #include "Ions.h" @@ -95,7 +95,7 @@ class MGmol : public MGmolInterface SpreadPenaltyInterface* spread_penalty_; - DMStrategy* dm_strategy_; + DMStrategy* dm_strategy_; HDFrestart* h5f_file_; diff --git a/src/MVPSolver.cc b/src/MVPSolver.cc index 8ac6a0be..63ddf690 100644 --- a/src/MVPSolver.cc +++ b/src/MVPSolver.cc @@ -184,6 +184,7 @@ int MVPSolver::solve(OrbitalsType& orbitals) ProjectedMatrices* projmatrices = dynamic_cast*>( orbitals.getProjMatrices()); + assert(projmatrices); s11 = projmatrices->getGramMatrix(); } diff --git a/src/MVP_DMStrategy.cc b/src/MVP_DMStrategy.cc index 2a6bebe1..8712ae92 100644 --- a/src/MVP_DMStrategy.cc +++ b/src/MVP_DMStrategy.cc @@ -24,17 +24,17 @@ template MVP_DMStrategy::MVP_DMStrategy(MPI_Comm comm, ostream& os, Ions& ions, Rho* rho, Energy* energy, Electrostatic* electrostat, - MGmol* mgmol_strategy, OrbitalsType* orbitals, + MGmol* mgmol_strategy, + const std::vector>& overlappingGids, ProjectedMatricesInterface* proj_matrices, const bool use_old_dm) - : orbitals_(orbitals), - proj_matrices_(proj_matrices), + : proj_matrices_(proj_matrices), comm_(comm), os_(os), ions_(ions), rho_(rho), energy_(energy), electrostat_(electrostat), - global_indexes_(orbitals->getOverlappingGids()), + global_indexes_(overlappingGids), mgmol_strategy_(mgmol_strategy), use_old_dm_(use_old_dm) { @@ -43,7 +43,7 @@ MVP_DMStrategy::MVP_DMStrategy(MPI_Comm comm, } template -int MVP_DMStrategy::update() +int MVP_DMStrategy::update(OrbitalsType& orbitals) { Control& ct = *(Control::instance()); if (onpe0 && ct.verbose > 2) @@ -56,7 +56,7 @@ int MVP_DMStrategy::update() electrostat_, mgmol_strategy_, ct.numst, ct.occ_width, global_indexes_, ct.dm_inner_steps, use_old_dm_); - return solver.solve(*orbitals_); + return solver.solve(orbitals); } template diff --git a/src/MVP_DMStrategy.h b/src/MVP_DMStrategy.h index 61bb6157..4b7e3be0 100644 --- a/src/MVP_DMStrategy.h +++ b/src/MVP_DMStrategy.h @@ -24,10 +24,9 @@ class Ions; class Electrostatic; template -class MVP_DMStrategy : public DMStrategy +class MVP_DMStrategy : public DMStrategy { private: - OrbitalsType* orbitals_; ProjectedMatricesInterface* proj_matrices_; MPI_Comm comm_; @@ -46,11 +45,11 @@ class MVP_DMStrategy : public DMStrategy MVP_DMStrategy(MPI_Comm comm, std::ostream& os, Ions& ions, Rho* rho, Energy* energy, Electrostatic* electrostat, MGmol* mgmol_strategy, - OrbitalsType* orbitals, ProjectedMatricesInterface* proj_matrices, - const bool use_old_dm); + const std::vector>& overlappingGids, + ProjectedMatricesInterface* proj_matrices, const bool use_old_dm); - void initialize() override{}; - int update() override; + void initialize(OrbitalsType& orbitals) override{}; + int update(OrbitalsType& orbitals) override; // H is updated with MVP loop, so no need to compute it outside bool needH() const override { return false; } diff --git a/src/NonOrthoDMStrategy.cc b/src/NonOrthoDMStrategy.cc index 1d8a5bbc..5452592f 100644 --- a/src/NonOrthoDMStrategy.cc +++ b/src/NonOrthoDMStrategy.cc @@ -13,15 +13,15 @@ #include "LocGridOrbitals.h" #include "ProjectedMatricesInterface.h" -template -NonOrthoDMStrategy::NonOrthoDMStrategy( - T* orbitals, ProjectedMatricesInterface* proj_matrices, const double mix) - : orbitals_(orbitals), proj_matrices_(proj_matrices), mix_(mix) +template +NonOrthoDMStrategy::NonOrthoDMStrategy( + ProjectedMatricesInterface* proj_matrices, const double mix) + : proj_matrices_(proj_matrices), mix_(mix) { } -template -void NonOrthoDMStrategy::initialize() +template +void NonOrthoDMStrategy::initialize(OrbitalsType& orbitals) { Control& ct = *(Control::instance()); MGmol_MPI& mmpi = *(MGmol_MPI::instance()); @@ -31,11 +31,11 @@ void NonOrthoDMStrategy::initialize() (*MPIdata::sout) << "NonOrthoDMStrategy::initialize()..." << std::endl; } - proj_matrices_->updateDM(orbitals_->getIterativeIndex()); + proj_matrices_->updateDM(orbitals.getIterativeIndex()); } -template -int NonOrthoDMStrategy::update() +template +int NonOrthoDMStrategy::update(OrbitalsType& orbitals) { assert(proj_matrices_ != nullptr); @@ -59,11 +59,11 @@ int NonOrthoDMStrategy::update() if (mix_ < 1.) proj_matrices_->saveDM(); // compute new density matrix - proj_matrices_->updateDM(orbitals_->getIterativeIndex()); + proj_matrices_->updateDM(orbitals.getIterativeIndex()); if (mix_ < 1.) { - proj_matrices_->updateDMwithRelax(mix_, orbitals_->getIterativeIndex()); + proj_matrices_->updateDMwithRelax(mix_, orbitals.getIterativeIndex()); } if (ct.verbose > 2) diff --git a/src/NonOrthoDMStrategy.h b/src/NonOrthoDMStrategy.h index 44c2e788..a661ebee 100644 --- a/src/NonOrthoDMStrategy.h +++ b/src/NonOrthoDMStrategy.h @@ -13,20 +13,19 @@ #include "DMStrategy.h" #include "ProjectedMatricesInterface.h" -template -class NonOrthoDMStrategy : public DMStrategy +template +class NonOrthoDMStrategy : public DMStrategy { private: - T* orbitals_; ProjectedMatricesInterface* proj_matrices_; const double mix_; public: - NonOrthoDMStrategy(T* orbitals, ProjectedMatricesInterface* proj_matrices, - const double mix); + NonOrthoDMStrategy( + ProjectedMatricesInterface* proj_matrices, const double mix); - void initialize() override; - int update() override; + void initialize(OrbitalsType& orbitals) override; + int update(OrbitalsType& orbitals) override; bool needH() const override { return true; } diff --git a/src/Orbitals.h b/src/Orbitals.h index 3d769af0..d58f46bb 100644 --- a/src/Orbitals.h +++ b/src/Orbitals.h @@ -73,5 +73,8 @@ class Orbitals = out_restart_info > 3 ? H5T_NATIVE_DOUBLE : H5T_NATIVE_FLOAT; return dtype_id; } + +private: + Orbitals& operator=(const Orbitals& orbitals); }; #endif diff --git a/src/OrbitalsExtrapolationOrder2.cc b/src/OrbitalsExtrapolationOrder2.cc index 4550637f..869dcd47 100644 --- a/src/OrbitalsExtrapolationOrder2.cc +++ b/src/OrbitalsExtrapolationOrder2.cc @@ -85,7 +85,7 @@ void OrbitalsExtrapolationOrder2::extrapolate_orbitals( { // DM (if not recomputed from scratch) // is consistant with orthonormal set of orbitals... - (*orbitals)->orthonormalizeLoewdin(); + if (ct.fullyOccupied()) (*orbitals)->orthonormalizeLoewdin(); } } diff --git a/src/PolakRibiereSolver.cc b/src/PolakRibiereSolver.cc index a2a26ae9..6b5f88de 100644 --- a/src/PolakRibiereSolver.cc +++ b/src/PolakRibiereSolver.cc @@ -30,7 +30,8 @@ PolakRibiereSolver::PolakRibiereSolver( Hamiltonian* hamiltonian, ProjectedMatricesInterface* proj_matrices, Energy* energy, Electrostatic* electrostat, MGmol* mgmol_strategy, Ions& ions, - Rho* rho, DMStrategy* dm_strategy, std::ostream& os) + Rho* rho, DMStrategy* dm_strategy, + std::ostream& os) : mgmol_strategy_(mgmol_strategy), hamiltonian_(hamiltonian), proj_matrices_(proj_matrices), @@ -543,7 +544,7 @@ int PolakRibiereSolver::solve(OrbitalsType& orbitals, mgmol_strategy_->updateHmatrix(orbitals, ions); // compute new density matrix - dm_success = dm_strategy_->update(); + dm_success = dm_strategy_->update(orbitals); skip_wolfe_cond = false; } diff --git a/src/PolakRibiereSolver.h b/src/PolakRibiereSolver.h index 5f1a2421..fd96ded9 100644 --- a/src/PolakRibiereSolver.h +++ b/src/PolakRibiereSolver.h @@ -10,6 +10,7 @@ #ifndef MGMOL_PolakRibiereSolver_H #define MGMOL_PolakRibiereSolver_H +#include "DMStrategy.h" #include "Energy.h" #include "Hamiltonian.h" #include "MGmol.h" @@ -20,7 +21,6 @@ class Ions; class Electrostatic; class ProjectedMatricesInterface; -class DMStrategy; template class PolakRibiereSolver @@ -37,7 +37,7 @@ class PolakRibiereSolver Electrostatic* electrostat_; Ions& ions_; Rho* rho_; - DMStrategy* dm_strategy_; + DMStrategy* dm_strategy_; std::ostream& os_; @@ -87,8 +87,8 @@ class PolakRibiereSolver PolakRibiereSolver(Hamiltonian* hamiltonian, ProjectedMatricesInterface* proj_matrices, Energy* energy, Electrostatic* electrostat, MGmol* mgmol_strategy, - Ions& ions, Rho* rho, DMStrategy* dm_strategy, - std::ostream& os); + Ions& ions, Rho* rho, + DMStrategy* dm_strategy, std::ostream& os); ~PolakRibiereSolver(); diff --git a/src/ProjectedMatrices.cc b/src/ProjectedMatrices.cc index 6d6b6339..bdfe7786 100644 --- a/src/ProjectedMatrices.cc +++ b/src/ProjectedMatrices.cc @@ -104,7 +104,9 @@ ProjectedMatrices::ProjectedMatrices( gm_(new GramMatrix(ndim)) { MGmol_MPI& mmpi = *(MGmol_MPI::instance()); - if (mmpi.instancePE0()) + Control& ct = *(Control::instance()); + + if (mmpi.instancePE0() && ct.verbose > 1) { std::cout << "New ProjectedMatrices with MatrixType: " << getMatrixType() << std::endl; @@ -164,7 +166,7 @@ void ProjectedMatrices>::setupMPI( MPI_Comm comm = mmpi.commSpin(); DistMatrix2SquareLocalMatrices::setup( - comm, global_indexes, dm_->getMatrix()); + comm, global_indexes, gm_->getMatrix()); LocalMatrices2DistMatrix::setup(comm, global_indexes); } @@ -262,6 +264,12 @@ void ProjectedMatrices::applyInvS( template void ProjectedMatrices::setDMto2InvS() { + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + Control& ct = *(Control::instance()); + + if (mmpi.instancePE0() && ct.verbose > 1) + std::cout << "ProjectedMatrices::setDMto2InvS()..." << std::endl; + dm_->setto2InvS(gm_->getInverse(), gm_->getAssociatedOrbitalsIndex()); } diff --git a/src/ProjectedMatrices.h b/src/ProjectedMatrices.h index 3257ba7c..702637f4 100644 --- a/src/ProjectedMatrices.h +++ b/src/ProjectedMatrices.h @@ -252,7 +252,11 @@ class ProjectedMatrices : public ProjectedMatricesInterface gm_->computeInverse(); } - const MatrixType& getGramMatrix() { return gm_->getMatrix(); } + const MatrixType& getGramMatrix() + { + assert(gm_ != nullptr); + return gm_->getMatrix(); + } void getOccupations(std::vector& occ) const override; void setOccupations(const std::vector& occ); @@ -264,7 +268,7 @@ class ProjectedMatrices : public ProjectedMatricesInterface SquareLocalMatrices& mat) override; void applyInvS(MatrixType& mat) { - assert(gm_ != 0); + assert(gm_ != nullptr); gm_->applyInv(mat); } diff --git a/src/md.cc b/src/md.cc index ec288839..776886c4 100644 --- a/src/md.cc +++ b/src/md.cc @@ -31,8 +31,8 @@ #include "ProjectedMatricesMehrstellen.h" #include "ProjectedMatricesSparse.h" #include "Rho.h" -#include "mgmol_Signal.h" #include "SpreadsAndCenters.h" +#include "mgmol_Signal.h" #include "tools.h" #include @@ -72,7 +72,7 @@ void MGmol::preWFextrapolation() if (ct.OuterSolver() == OuterSolverType::ABPG || ct.OuterSolver() == OuterSolverType::NLCG) { - if (ct.dm_mix < 1.) proj_matrices_->stripDM(); + if (ct.dm_mix < 1. || !ct.fullyOccupied()) proj_matrices_->stripDM(); } } @@ -80,22 +80,21 @@ template void MGmol::postWFextrapolation(OrbitalsType* orbitals) { Control& ct = *(Control::instance()); - if (ct.isLocMode()) - { - orbitals->computeGramAndInvS(); - } + + orbitals->computeGramAndInvS(); if (ct.OuterSolver() == OuterSolverType::ABPG || ct.OuterSolver() == OuterSolverType::NLCG) { - if (ct.dm_mix < 1.) + if (ct.dm_mix < 1. || !ct.fullyOccupied()) proj_matrices_->dressupDM(); else proj_matrices_->setDMto2InvS(); } else { - if (orbitals_extrapol_->extrapolatedH()) dm_strategy_->update(); + if (orbitals_extrapol_->extrapolatedH()) + dm_strategy_->update(*orbitals); } } @@ -382,7 +381,7 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) // need to reset a few things as we just read new orbitals (*orbitals)->computeGramAndInvS(); - dm_strategy_->update(); + dm_strategy_->update(*current_orbitals_); } DFTsolver::setItCountLarge(); @@ -462,7 +461,7 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) move_orbitals(orbitals); (*orbitals)->computeGramAndInvS(); - dm_strategy_->update(); + dm_strategy_->update(**orbitals); // reduce number of steps to keep total run time about the // same @@ -675,7 +674,8 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) } template -OrbitalsType* MGmol::loadOrbitalFromRestartFile(const std::string filename) +OrbitalsType* MGmol::loadOrbitalFromRestartFile( + const std::string filename) { MGmol_MPI& mmpi(*(MGmol_MPI::instance())); Control& ct = *(Control::instance()); @@ -684,23 +684,25 @@ OrbitalsType* MGmol::loadOrbitalFromRestartFile(const std::string /* For now, we only consider double-precision hdf5 I/O. */ assert(ct.restart_info > 3); - assert((ct.AtomsDynamic() == AtomsDynamicType::MD) || (ct.AtomsDynamic() == AtomsDynamicType::Quench)); + assert((ct.AtomsDynamic() == AtomsDynamicType::MD) + || (ct.AtomsDynamic() == AtomsDynamicType::Quench)); HDFrestart h5file(filename, myPEenv, ct.out_restart_file_type); int ierr; - OrbitalsType *restart_orbitals = nullptr; + OrbitalsType* restart_orbitals = nullptr; /* This corresponds to MGmol::initial */ { // Copy from current orbital, instead of constructing brand-new one - restart_orbitals = new OrbitalsType("ForLoading", *current_orbitals_, false); + restart_orbitals + = new OrbitalsType("ForLoading", *current_orbitals_, false); /* This corresponds to MGmol::read_restart_data */ { ierr = restart_orbitals->read_func_hdf5(h5file); - } // read_restart_data - } // initial() + } // read_restart_data + } // initial() /* This corresponds to MGmol::md */ { @@ -723,18 +725,19 @@ OrbitalsType* MGmol::loadOrbitalFromRestartFile(const std::string */ if (flag_extrapolated_data) { - orbitals_extrapol_ = OrbitalsExtrapolationFactory::create( - ct.WFExtrapolation()); + orbitals_extrapol_ + = OrbitalsExtrapolationFactory::create( + ct.WFExtrapolation()); if (onpe0) os_ << "Create new orbitals_minus1..." << std::endl; orbitals_extrapol_->setupPreviousOrbitals(&restart_orbitals, - proj_matrices_, lrs_, local_cluster_, currentMasks_, - corrMasks_, h5file); + proj_matrices_, lrs_, local_cluster_, currentMasks_, corrMasks_, + h5file); } /* main workflow delete h5f_file_ here, meaning the loading is over. */ - } // md() + } // md() ierr = h5file.close(); mmpi.allreduce(&ierr, 1, MPI_MIN); @@ -749,7 +752,8 @@ OrbitalsType* MGmol::loadOrbitalFromRestartFile(const std::string /* In returning restart_orbitals, we hope that the wavefunctions in restart_orbitals are all set. - At least the following functions should return proper data loaded from the file: + At least the following functions should return proper data loaded from + the file: restart_orbitals->getLocNumpt() restart_orbitals->chromatic_number() diff --git a/src/quench.cc b/src/quench.cc index c02f1f73..de01dad3 100644 --- a/src/quench.cc +++ b/src/quench.cc @@ -176,7 +176,7 @@ void MGmol::resetProjectedMatricesAndDM( proj_matrices_->updateThetaAndHB(); // reset DM using chosen strategy - dm_strategy_->initialize(); + dm_strategy_->initialize(orbitals); } // try to use some rotations to avoid degeneracies From 2289afefc47c31a8dfd187e1e5f271b06937c77c Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Thu, 30 May 2024 14:25:24 -0400 Subject: [PATCH 07/11] Fix issue with single hdf5 output file (#241) * clean up some related loops * have MVP test write single hdf5 file for testing --- src/Ions.cc | 315 ++++++++-------------------------------------- tests/MVP/mvp.cfg | 3 +- tests/MVP/test.py | 8 +- 3 files changed, 59 insertions(+), 267 deletions(-) diff --git a/src/Ions.cc b/src/Ions.cc index e2cc9916..53ee0cf4 100644 --- a/src/Ions.cc +++ b/src/Ions.cc @@ -521,114 +521,6 @@ double Ions::energyDiff(const short bc[3]) const return energy; // Hartree } -// Associate each ion with the PE where they are centered -// by setting a flag "here" -// void Ions::associate2PE() -//{ -// Control& ct(*(Control::instance())); -// Mesh* mymesh = Mesh::instance(); -// const pb::PEenv& myPEenv = mymesh->peenv(); -// const pb::Grid& mygrid = mymesh->grid(); -// -// double div_lattice[3]; -// for (short i = 0; i < 3; i++) -// div_lattice[i] = lattice_[i] / (double)(myPEenv.n_mpi_task(i)); -// -// double offset[3]; -// for (short i = 0; i < 3; i++) -// offset[i] = (double)myPEenv.my_mpi(i) * div_lattice[i]; -// -//#ifdef DEBUG -// int isum2 = 0; -// (*MPIdata::sout) << " offset=(" << offset[0] << "," << offset[1] << "," -// << offset[2] << ")" << std::endl; -//#endif -// -// for (short i = 0; i < 3; i++) -// if (myPEenv.n_mpi_task(i) == (myPEenv.my_mpi(i) + 1)) -// div_lattice[i] = lattice_[i] - offset[i]; -// -// const double origin[3] -// = { mygrid.origin(0), mygrid.origin(1), mygrid.origin(2) }; -// const double end[3] = { origin[0] + lattice_[0], origin[1] + lattice_[1], -// origin[2] + lattice_[2] }; -// -// int isum = 0; -// -// // Loop over ions -// local_ions_.clear(); -// std::vector::iterator ion = list_ions_.begin(); -// while (ion != list_ions_.end()) -// { -// double t[3]; -// for (short i = 0; i < 3; i++) -// { -// t[i] = (*ion)->position(i); -// while (t[i] >= end[i]) -// t[i] -= lattice_[i]; -// while (t[i] < origin[i]) -// t[i] += lattice_[i]; -// t[i] -= origin[i]; -// t[i] -= offset[i]; -// } -// -//#if DEBUG -// (*MPIdata::sout) << " t=(" << t[0] << "," << t[1] << "," << t[2] << -// ")" -// << std::endl; -//#endif -// -// if ((t[0] >= 0. && t[0] < (div_lattice[0])) -// && (t[1] >= 0. && t[1] < (div_lattice[1])) -// && (t[2] >= 0. && t[2] < (div_lattice[2]))) -// { -// -// (*ion)->set_here(true); -// -// local_ions_.push_back(*ion); -// -// isum++; -// } -// else -// { -// (*ion)->set_here(false); -// } -//#ifndef NDEBUG -// if ((*ion)->here()) -// { -// (*MPIdata::sout) << " Ion " << (*ion)->name() << " centered on PE -// " -// << myPEenv.mytask() << std::endl; -// } -//#endif -// ion++; -// } -// // Test all atoms associated to each processors sum up to -// // total number of atoms -// -// MGmol_MPI& mmpi = *(MGmol_MPI::instance()); -// mmpi.allreduce(&isum, 1, MPI_SUM); -// if (onpe0) -// { -// (*MPIdata::sout) << " num_ions=" << num_ions_ << std::endl; -// (*MPIdata::sout) << " isum=" << isum << std::endl; -// } -// if (isum != num_ions_) -// { -// (*MPIdata::sout) << " num_ions != isum !!!!" << std::endl; -// ct.global_exit(2); -// } -// -//#ifdef DEBUG -// mmpi.allreduce(&isum2, 1, MPI_SUM); -// (*MPIdata::sout) << " Isum2=" << isum2 << std::endl; -// assert(isum2 == num_ions_); -//#endif -// -// // if( onpe0 ) -// // (*MPIdata::sout)<<"Ions::associate2PE() done..."<::const_iterator ion = local_ions_.begin(); - while (ion != local_ions_.end()) + for (auto& ion : local_ions_) { - data.push_back((*ion)->name()); - ion++; + data.push_back(ion->name()); } } @@ -851,20 +741,21 @@ void Ions::writeLockedAtomNames(HDFrestart& h5f_file) } else { - std::vector::const_iterator ion = local_ions_.begin(); - while (ion != local_ions_.end()) + for (auto& ion : local_ions_) { - if ((*ion)->locked()) data.push_back((*ion)->name()); - ion++; + if (ion->locked()) data.push_back(ion->name()); } } - hid_t file_id = h5f_file.file_id(); - if (file_id >= 0) + if (!data.empty()) { - std::string datasetname("/LockedAtomsNames"); - std::string empty; - writeData2d(h5f_file, datasetname, data, 1, empty); + hid_t file_id = h5f_file.file_id(); + if (file_id >= 0) + { + std::string datasetname("/LockedAtomsNames"); + std::string empty; + writeData2d(h5f_file, datasetname, data, 1, empty); + } } } @@ -888,11 +779,9 @@ void Ions::writeAtomicIDs(HDFrestart& h5f_file) } else { - std::vector::const_iterator ion = local_ions_.begin(); - while (ion != local_ions_.end()) + for (auto& ion : local_ions_) { - data.push_back((*ion)->index()); - ion++; + data.push_back(ion->index()); } } @@ -925,11 +814,9 @@ void Ions::writeAtomicNLprojIDs(HDFrestart& h5f_file) } else { - std::vector::const_iterator ion = local_ions_.begin(); - while (ion != local_ions_.end()) + for (auto& ion : local_ions_) { - data.push_back((*ion)->nlprojid()); - ion++; + data.push_back(ion->nlprojid()); } } @@ -961,13 +848,11 @@ void Ions::writePositions(HDFrestart& h5f_file) } else { - std::vector::const_iterator ion = local_ions_.begin(); - while (ion != local_ions_.end()) + for (auto& ion : local_ions_) { - data.push_back((*ion)->position(0)); - data.push_back((*ion)->position(1)); - data.push_back((*ion)->position(2)); - ion++; + data.push_back(ion->position(0)); + data.push_back(ion->position(1)); + data.push_back(ion->position(2)); } } @@ -1124,13 +1009,11 @@ void Ions::writeVelocities(HDFrestart& h5f_file) } else { - std::vector::const_iterator ion = local_ions_.begin(); - while (ion != local_ions_.end()) + for (auto& ion : local_ions_) { - data.push_back((*ion)->velocity(0)); - data.push_back((*ion)->velocity(1)); - data.push_back((*ion)->velocity(2)); - ion++; + data.push_back(ion->velocity(0)); + data.push_back(ion->velocity(1)); + data.push_back(ion->velocity(2)); } } @@ -1162,13 +1045,11 @@ void Ions::writeRandomStates(HDFrestart& h5f_file) } else { - std::vector::const_iterator ion = local_ions_.begin(); - while (ion != local_ions_.end()) + for (auto& ion : local_ions_) { - data.push_back((*ion)->randomState(0)); - data.push_back((*ion)->randomState(1)); - data.push_back((*ion)->randomState(2)); - ion++; + data.push_back(ion->randomState(0)); + data.push_back(ion->randomState(1)); + data.push_back(ion->randomState(2)); } } @@ -2556,11 +2437,9 @@ void Ions::gatherNames( { std::vector local_names; - std::vector::const_iterator ion = local_ions_.begin(); - while (ion != local_ions().end()) + for (auto& ion : local_ions_) { - local_names.push_back((*ion)->name()); - ++ion; + local_names.push_back(ion->name()); } // gather data to PE root @@ -2578,11 +2457,9 @@ void Ions::gatherLockedNames( { std::vector local_names; - std::vector::const_iterator ion = local_ions_.begin(); - while (ion != local_ions().end()) + for (auto& ion : local_ions_) { - if ((*ion)->locked()) local_names.push_back((*ion)->name()); - ++ion; + if (ion->locked()) local_names.push_back(ion->name()); } // gather data to PE root @@ -2600,11 +2477,9 @@ void Ions::gatherIndexes( { std::vector local_indexes; - std::vector::const_iterator ion = local_ions_.begin(); - while (ion != local_ions_.end()) + for (auto& ion : local_ions_) { - local_indexes.push_back((*ion)->index()); - ++ion; + local_indexes.push_back(ion->index()); } // gather data to PE root @@ -2622,11 +2497,9 @@ void Ions::gatherNLprojIds( { std::vector local_nlprojids; - std::vector::const_iterator ion = local_ions_.begin(); - while (ion != local_ions_.end()) + for (auto& ion : local_ions_) { - local_nlprojids.push_back((*ion)->nlprojid()); - ++ion; + local_nlprojids.push_back(ion->nlprojid()); } // gather data to PE root @@ -2644,11 +2517,9 @@ void Ions::gatherAtomicNumbers( { std::vector local_atnumbers; - std::vector::const_iterator ion = local_ions_.begin(); - while (ion != local_ions_.end()) + for (auto& ion : local_ions_) { - local_atnumbers.push_back((*ion)->atomic_number()); - ++ion; + local_atnumbers.push_back(ion->atomic_number()); } // gather data to PE root @@ -2666,13 +2537,11 @@ void Ions::gatherRandStates(std::vector& rstates, { std::vector local_rstates; - std::vector::const_iterator ion = local_ions_.begin(); - while (ion != local_ions().end()) + for (auto& ion : local_ions_) { - local_rstates.push_back((*ion)->randomState(0)); - local_rstates.push_back((*ion)->randomState(1)); - local_rstates.push_back((*ion)->randomState(2)); - ++ion; + local_rstates.push_back(ion->randomState(0)); + local_rstates.push_back(ion->randomState(1)); + local_rstates.push_back(ion->randomState(2)); } // gather data to PE root @@ -2690,17 +2559,14 @@ void Ions::gatherPositions( { std::vector local_positions; - std::vector::const_iterator ion = local_ions_.begin(); - while (ion != local_ions().end()) + for (auto& ion : local_ions_) { // get position of local ion double position[3]; - (*ion)->getPosition(&position[0]); + ion->getPosition(&position[0]); local_positions.push_back(position[0]); local_positions.push_back(position[1]); local_positions.push_back(position[2]); - - ++ion; } // gather data to PE root @@ -2718,17 +2584,14 @@ void Ions::gatherForces( { std::vector local_forces; - std::vector::const_iterator ion = local_ions_.begin(); - while (ion != local_ions().end()) + for (auto& ion : local_ions_) { // get position of local ion double force[3]; - (*ion)->getForce(&force[0]); + ion->getForce(&force[0]); local_forces.push_back(force[0]); local_forces.push_back(force[1]); local_forces.push_back(force[2]); - - ++ion; } // gather data to PE root @@ -2746,14 +2609,11 @@ void Ions::gatherVelocities( { std::vector local_velocities; - std::vector::const_iterator ion = local_ions_.begin(); - while (ion != local_ions().end()) + for (auto& ion : local_ions_) { - local_velocities.push_back((*ion)->velocity(0)); - local_velocities.push_back((*ion)->velocity(1)); - local_velocities.push_back((*ion)->velocity(2)); - - ++ion; + local_velocities.push_back(ion->velocity(0)); + local_velocities.push_back(ion->velocity(1)); + local_velocities.push_back(ion->velocity(2)); } // gather data to PE root @@ -2846,81 +2706,6 @@ bool Ions::hasLockedAtoms() const return (flag == 1); } -#if 0 -void Ions::syncNames(const int nions, std::vector& local_names, vector& names) -{ - MGmol_MPI& mmpi ( *(MGmol_MPI::instance()) ); - MPI_Comm& comm=mmpi.comm(); - // prepare to gather data - int npes = mmpi.size(); - int num_ions = nions; - std::vectorrecvcnts(npes); - std::vector::iterator rcv; - std::vectordisp(npes); - std::vector::iterator disp_it; - int tot = 0; - int totchars; - // Gather data for atom names - // first get length of each name - std::vectorloc_name_len; - std::vectornameLen(num_ions,0); - tot=0; - for(std::vector::iterator str=local_names.begin(); str!=local_names.end(); ++str){ - //tot+= *str.size(); - std::string s = *str; - tot+= s.length(); - loc_name_len.push_back(s.length()); - } - mmpi.allGatherV(loc_name_len, nameLen); - - std::vectorchar_loc_names(tot); - int idx = 0; - for(std::vector::iterator str=local_names.begin(); str!=local_names.end(); ++str){ - //tot+= *str.size(); - std::string s = *str; - memcpy(&char_loc_names[idx], s.c_str(), s.size()); - idx+=s.size(); - } - assert(idx == tot); - - // gather data for atom names - //get recv counts - rcv = recvcnts.begin(); - MPI_Allgather(&tot, 1, MPI_INT, &(*rcv), 1, MPI_INT, comm); - - //get displacements - totchars = recvcnts[npes-1]; - disp[0] = 0; - for(int pos = 1; pos < npes; pos++){ - disp[pos] = disp[pos-1] + recvcnts[pos-1]; - totchars+=recvcnts[pos-1]; - } - // gather atomic names - std::vector names_data(totchars); - std::vector::iterator locnames_it=char_loc_names.begin(); - std::vector::iterator names_it=names_data.begin(); - disp_it=disp.begin(); - rcv = recvcnts.begin(); - MPI_Allgatherv(&(*locnames_it), tot, MPI_CHAR, &(*names_it), &(*rcv), &(*disp_it), MPI_CHAR, comm); - -// std::vector names; - int pos = 0; - for(int i=0; itol: sys.exit(1) +for line in lines: + if line.count(b'HDF5-DIAG') and line.count(b'Error'): + print(line) + print("Found HDF5 error") + sys.exit(1) + niterations = len(energies) print("MVP solver ran for {} iterations".format(niterations)) if niterations>180: From 3ed3b624228bd194a0e7d5398af05580bc9d3e4b Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Thu, 30 May 2024 16:40:36 -0400 Subject: [PATCH 08/11] Clean up setup phase (#242) --- src/MGmol.cc | 13 ++++++++++++- src/MGmol.h | 14 ++++++++------ src/main.cc | 27 +++++++++------------------ src/mgmol_run.cc | 27 +-------------------------- src/mgmol_run.h | 2 +- src/setup.cc | 36 ++++++++++++++++++++++++++++++++++++ 6 files changed, 67 insertions(+), 52 deletions(-) diff --git a/src/MGmol.cc b/src/MGmol.cc index fb6c64f6..b6a95777 100644 --- a/src/MGmol.cc +++ b/src/MGmol.cc @@ -110,7 +110,10 @@ extern Timer updateCenters_tm; std::set Signal::recv_; template -MGmol::MGmol(MPI_Comm comm, std::ostream& os) : os_(os) +MGmol::MGmol(MPI_Comm comm, std::ostream& os, + std::string input_filename, std::string lrs_filename, + std::string constraints_filename) + : os_(os) { comm_ = comm; @@ -136,6 +139,14 @@ MGmol::MGmol(MPI_Comm comm, std::ostream& os) : os_(os) forces_ = nullptr; energy_ = nullptr; + + setupFromInput(input_filename); + + setupLRs(lrs_filename); + + setupConstraintsFromInput(constraints_filename); + + setup(); } template diff --git a/src/MGmol.h b/src/MGmol.h index d43be062..a44a37f8 100644 --- a/src/MGmol.h +++ b/src/MGmol.h @@ -145,7 +145,12 @@ class MGmol : public MGmolInterface template int initial(); void initialMasks(); - int setupLRsFromInput(const std::string input_file); + int setupLRsFromInput(const std::string filename); + + void setup(); + int setupLRs(const std::string input_file) override; + int setupFromInput(const std::string input_file) override; + int setupConstraintsFromInput(const std::string input_file) override; // timers static Timer total_tm_; @@ -166,7 +171,8 @@ class MGmol : public MGmolInterface public: Electrostatic* electrostat_; - MGmol(MPI_Comm comm, std::ostream& os); + MGmol(MPI_Comm comm, std::ostream& os, std::string input_filename, + std::string lrs_filename, std::string constraints_filename); ~MGmol() override; @@ -274,10 +280,6 @@ class MGmol : public MGmolInterface void set_forces(std::vector>& f); int nions() { return ions_->getNumIons(); } double getTotalEnergy(); - void setup(); - int setupLRs(const std::string input_file) override; - int setupFromInput(const std::string input_file) override; - int setupConstraintsFromInput(const std::string input_file) override; void cleanup(); void geomOptimSetup(); void geomOptimQuench(); diff --git a/src/main.cc b/src/main.cc index ae50c299..f6f396d2 100644 --- a/src/main.cc +++ b/src/main.cc @@ -108,6 +108,11 @@ int main(int argc, char** argv) int ret = ct.checkOptions(); if (ret < 0) return ret; + unsigned ngpts[3] = { ct.ngpts_[0], ct.ngpts_[1], ct.ngpts_[2] }; + double origin[3] = { ct.ox_, ct.oy_, ct.oz_ }; + const double cell[3] = { ct.lx_, ct.ly_, ct.lz_ }; + Mesh::setup(mmpi.commSpin(), ngpts, origin, cell, ct.lap_type); + mmpi.bcastGlobal(input_filename); mmpi.bcastGlobal(lrs_filename); @@ -115,28 +120,14 @@ int main(int argc, char** argv) { MGmolInterface* mgmol; if (ct.isLocMode()) - mgmol = new MGmol(global_comm, *MPIdata::sout); + mgmol = new MGmol(global_comm, *MPIdata::sout, + input_filename, lrs_filename, constraints_filename); else - mgmol - = new MGmol(global_comm, *MPIdata::sout); - - unsigned ngpts[3] = { ct.ngpts_[0], ct.ngpts_[1], ct.ngpts_[2] }; - double origin[3] = { ct.ox_, ct.oy_, ct.oz_ }; - const double cell[3] = { ct.lx_, ct.ly_, ct.lz_ }; - Mesh::setup(mmpi.commSpin(), ngpts, origin, cell, ct.lap_type); - - mgmol->setupFromInput(input_filename); - - if (ct.isLocMode() || ct.init_loc == 1) mgmol->setupLRs(lrs_filename); - - mgmol->setupConstraintsFromInput(constraints_filename); - - mgmol_setup(); + mgmol = new MGmol(global_comm, *MPIdata::sout, + input_filename, lrs_filename, constraints_filename); if (!tcheck) { - mgmol->setup(); - mgmol->run(); } else diff --git a/src/mgmol_run.cc b/src/mgmol_run.cc index 0a791d35..e746611c 100644 --- a/src/mgmol_run.cc +++ b/src/mgmol_run.cc @@ -86,38 +86,13 @@ int mgmol_init(MPI_Comm comm) return 0; } -int mgmol_setup() +int mgmol_check() { Control& ct = *(Control::instance()); Mesh* mymesh = Mesh::instance(); const pb::PEenv& myPEenv = mymesh->peenv(); MGmol_MPI& mmpi = *(MGmol_MPI::instance()); - ct.checkNLrange(); - - LocGridOrbitals::setDotProduct(ct.dot_product_type); - - if (!ct.short_sighted) - { - MatricesBlacsContext::instance().setup(mmpi.commSpin(), ct.numst); - - dist_matrix::DistMatrix::setBlockSize(64); - - dist_matrix::DistMatrix::setDefaultBlacsContext( - MatricesBlacsContext::instance().bcxt()); - - ReplicatedWorkSpace::instance().setup(ct.numst); - - dist_matrix::SparseDistMatrix::setNumTasksPerPartitioning( - 128); - - int npes = mmpi.size(); - setSparseDistMatriConsolidationNumber(npes); - } -#ifdef HAVE_MAGMA - ReplicatedMatrix::setMPIcomm(mmpi.commSpin()); -#endif - if (myPEenv.color() > 0) { std::cerr << "Code should be called with " << myPEenv.n_mpi_tasks() diff --git a/src/mgmol_run.h b/src/mgmol_run.h index 7e07e148..b091a17b 100644 --- a/src/mgmol_run.h +++ b/src/mgmol_run.h @@ -7,5 +7,5 @@ // This file is part of MGmol. For details, see https://github.com/llnl/mgmol. // Please also read this link https://github.com/llnl/mgmol/LICENSE int mgmol_init(MPI_Comm comm); -int mgmol_setup(); +int mgmol_check(); void mgmol_finalize(); diff --git a/src/setup.cc b/src/setup.cc index 035d7c8c..62d9f835 100644 --- a/src/setup.cc +++ b/src/setup.cc @@ -13,6 +13,10 @@ #include "LocGridOrbitals.h" #include "MGmol.h" #include "Potentials.h" +#include "ReplicatedWorkSpace.h" +#include "SparseDistMatrix.h" + +#include "mgmol_run.h" template int MGmol::setupFromInput(const std::string filename) @@ -51,6 +55,36 @@ int MGmol::setupFromInput(const std::string filename) ct.setTolEnergy(); ct.setSpreadRadius(); + ct.checkNLrange(); + + // now that we know the number of states, we can set a few other static + // data + if (!ct.short_sighted) + { + MatricesBlacsContext::instance().setup(mmpi.commSpin(), ct.numst); + + dist_matrix::DistMatrix::setBlockSize(64); + + dist_matrix::DistMatrix::setDefaultBlacsContext( + MatricesBlacsContext::instance().bcxt()); + + ReplicatedWorkSpace::instance().setup(ct.numst); + + dist_matrix::SparseDistMatrix::setNumTasksPerPartitioning( + 128); + + int npes = mmpi.size(); + setSparseDistMatriConsolidationNumber(npes); + } + +#ifdef HAVE_MAGMA + ReplicatedMatrix::setMPIcomm(mmpi.commSpin()); +#endif + + LocGridOrbitals::setDotProduct(ct.dot_product_type); + + mgmol_check(); + return 0; } @@ -59,6 +93,8 @@ int MGmol::setupLRs(const std::string filename) { Control& ct = *(Control::instance()); + if (!(ct.isLocMode() || ct.init_loc == 1)) return 0; + // create localization regions Mesh* mymesh = Mesh::instance(); const pb::Grid& mygrid = mymesh->grid(); From 98d0952fadaaea1e67642c696a335c7bcd45aec0 Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Thu, 30 May 2024 16:48:11 -0400 Subject: [PATCH 09/11] Reorder regression tests (#243) * run shorter ones first --- tests/CMakeLists.txt | 121 ++++++++++++++++++++++--------------------- 1 file changed, 62 insertions(+), 59 deletions(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 5ed818fa..4c68adb7 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -271,6 +271,7 @@ endif() target_compile_definitions(testAndersonMix PUBLIC TESTING) +# unit tests if(MGMOL_USE_HDF5P) add_test(NAME testHDF5P COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} @@ -345,15 +346,6 @@ if(${MAGMA_FOUND}) COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} ${CMAKE_CURRENT_BINARY_DIR}/testReplicatedMatrix) else() - add_test(NAME testShortSighted - COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/ShortSighted/test.py - ${MPIEXEC} --oversubscribe ${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 - ${CMAKE_CURRENT_SOURCE_DIR}/ShortSighted/coords.in - ${CMAKE_CURRENT_SOURCE_DIR}/ShortSighted/lrs.in - ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) add_test(NAME testProjectedMatrices COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/ProjectedMatrices/test.py ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} @@ -366,6 +358,28 @@ else() endif() #Regression tests +add_test(NAME testSiH4 + COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/SiH4/testSiH4.py + ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} + ${CMAKE_CURRENT_BINARY_DIR}/../src/mgmol-opt + ${CMAKE_CURRENT_SOURCE_DIR}/SiH4/mgmol.cfg + ${CMAKE_CURRENT_SOURCE_DIR}/SiH4/sih4.xyz + ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) +add_test(NAME testCl2 + COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/Cl2_ONCVPSP_LDA/test.py + ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} + ${CMAKE_CURRENT_BINARY_DIR}/../src/mgmol-opt + ${CMAKE_CURRENT_SOURCE_DIR}/Cl2_ONCVPSP_LDA/mgmol.cfg + ${CMAKE_CURRENT_SOURCE_DIR}/Cl2_ONCVPSP_LDA/coords.in + ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) +add_test(NAME testContinuumSolvent + COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/ContinuumSolvent/test.py + ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} + ${CMAKE_CURRENT_BINARY_DIR}/../src/mgmol-opt + ${CMAKE_CURRENT_SOURCE_DIR}/ContinuumSolvent/mgmol_diel.cfg + ${CMAKE_CURRENT_SOURCE_DIR}/ContinuumSolvent/coords.in + ${CMAKE_CURRENT_SOURCE_DIR}/ContinuumSolvent/lrs.in + ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) add_test(NAME Davidson COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/Davidson/test.py ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} @@ -373,33 +387,27 @@ add_test(NAME Davidson ${CMAKE_CURRENT_SOURCE_DIR}/Davidson/davidson.cfg ${CMAKE_CURRENT_SOURCE_DIR}/Davidson/coords.in ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) -add_test(NAME testMVP - COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/MVP/test.py - ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} - ${CMAKE_CURRENT_BINARY_DIR}/../src/mgmol-opt - ${CMAKE_CURRENT_SOURCE_DIR}/MVP/mvp.cfg - ${CMAKE_CURRENT_SOURCE_DIR}/MVP/coords.in - ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) -add_test(NAME testFatom - COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/Fatom/test.py +add_test(NAME testSpinO2 + COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/SpinO2/test.py ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} ${CMAKE_CURRENT_BINARY_DIR}/../src/mgmol-opt - ${CMAKE_CURRENT_SOURCE_DIR}/Fatom/mgmol.cfg - ${CMAKE_CURRENT_SOURCE_DIR}/Fatom/F.xyz + ${CMAKE_CURRENT_SOURCE_DIR}/SpinO2/mgmol_spin1.cfg + ${CMAKE_CURRENT_SOURCE_DIR}/SpinO2/coords.in ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) -add_test(NAME testSiH4 - COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/SiH4/testSiH4.py +add_test(NAME testSpinO2LDA + COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/SpinO2LDA/test.py ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} ${CMAKE_CURRENT_BINARY_DIR}/../src/mgmol-opt - ${CMAKE_CURRENT_SOURCE_DIR}/SiH4/mgmol.cfg - ${CMAKE_CURRENT_SOURCE_DIR}/SiH4/sih4.xyz + ${CMAKE_CURRENT_SOURCE_DIR}/SpinO2LDA/mgmol.cfg + ${CMAKE_CURRENT_SOURCE_DIR}/SpinO2LDA/coords.in + ${CMAKE_CURRENT_SOURCE_DIR}/SpinO2LDA/lrs.in ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) -add_test(NAME testCl2 - COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/Cl2_ONCVPSP_LDA/test.py +add_test(NAME testMVP + COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/MVP/test.py ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} ${CMAKE_CURRENT_BINARY_DIR}/../src/mgmol-opt - ${CMAKE_CURRENT_SOURCE_DIR}/Cl2_ONCVPSP_LDA/mgmol.cfg - ${CMAKE_CURRENT_SOURCE_DIR}/Cl2_ONCVPSP_LDA/coords.in + ${CMAKE_CURRENT_SOURCE_DIR}/MVP/mvp.cfg + ${CMAKE_CURRENT_SOURCE_DIR}/MVP/coords.in ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) add_test(NAME BandGapN2 COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/BandGapN2/test.py @@ -409,14 +417,6 @@ add_test(NAME BandGapN2 ${CMAKE_CURRENT_SOURCE_DIR}/BandGapN2/coords.in ${CMAKE_CURRENT_SOURCE_DIR}/BandGapN2/lrs.in ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) -add_test(NAME ReplicatedSP2 - COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/ReplicatedSP2/test.py - ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 3 ${MPIEXEC_PREFLAGS} - ${CMAKE_CURRENT_BINARY_DIR}/../src/mgmol-opt - ${CMAKE_CURRENT_SOURCE_DIR}/ReplicatedSP2/mgmol.cfg - ${CMAKE_CURRENT_SOURCE_DIR}/ReplicatedSP2/coords.in - ${CMAKE_CURRENT_SOURCE_DIR}/ReplicatedSP2/lrs.in - ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) add_test(NAME testMLWF COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/MLWF/test.py ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} @@ -432,6 +432,21 @@ add_test(NAME testSpreadPenalty ${CMAKE_CURRENT_SOURCE_DIR}/SpreadPenalty/coords.in ${CMAKE_CURRENT_SOURCE_DIR}/SpreadPenalty/lrs.in ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) +add_test(NAME testFatom + COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/Fatom/test.py + ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} + ${CMAKE_CURRENT_BINARY_DIR}/../src/mgmol-opt + ${CMAKE_CURRENT_SOURCE_DIR}/Fatom/mgmol.cfg + ${CMAKE_CURRENT_SOURCE_DIR}/Fatom/F.xyz + ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) +add_test(NAME ReplicatedSP2 + COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/ReplicatedSP2/test.py + ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 3 ${MPIEXEC_PREFLAGS} + ${CMAKE_CURRENT_BINARY_DIR}/../src/mgmol-opt + ${CMAKE_CURRENT_SOURCE_DIR}/ReplicatedSP2/mgmol.cfg + ${CMAKE_CURRENT_SOURCE_DIR}/ReplicatedSP2/coords.in + ${CMAKE_CURRENT_SOURCE_DIR}/ReplicatedSP2/lrs.in + ${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} @@ -466,30 +481,18 @@ add_test(NAME ChebyshevMVP ${CMAKE_CURRENT_SOURCE_DIR}/Chebyshev/cheb.cfg ${CMAKE_CURRENT_SOURCE_DIR}/Chebyshev/coords.in ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) -add_test(NAME testContinuumSolvent - COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/ContinuumSolvent/test.py - ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} - ${CMAKE_CURRENT_BINARY_DIR}/../src/mgmol-opt - ${CMAKE_CURRENT_SOURCE_DIR}/ContinuumSolvent/mgmol_diel.cfg - ${CMAKE_CURRENT_SOURCE_DIR}/ContinuumSolvent/coords.in - ${CMAKE_CURRENT_SOURCE_DIR}/ContinuumSolvent/lrs.in - ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) -add_test(NAME testSpinO2 - COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/SpinO2/test.py - ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} - ${CMAKE_CURRENT_BINARY_DIR}/../src/mgmol-opt - ${CMAKE_CURRENT_SOURCE_DIR}/SpinO2/mgmol_spin1.cfg - ${CMAKE_CURRENT_SOURCE_DIR}/SpinO2/coords.in - ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) -add_test(NAME testSpinO2LDA - COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/SpinO2LDA/test.py - ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} - ${CMAKE_CURRENT_BINARY_DIR}/../src/mgmol-opt - ${CMAKE_CURRENT_SOURCE_DIR}/SpinO2LDA/mgmol.cfg - ${CMAKE_CURRENT_SOURCE_DIR}/SpinO2LDA/coords.in - ${CMAKE_CURRENT_SOURCE_DIR}/SpinO2LDA/lrs.in - ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) +if(NOT ${MAGMA_FOUND}) + add_test(NAME testShortSighted + COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/ShortSighted/test.py + ${MPIEXEC} --oversubscribe ${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 + ${CMAKE_CURRENT_SOURCE_DIR}/ShortSighted/coords.in + ${CMAKE_CURRENT_SOURCE_DIR}/ShortSighted/lrs.in + ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) +endif() if(MGMOL_USE_HDF5P) target_include_directories(testHDF5P PRIVATE ${HDF5_INCLUDE_DIRS}) From 3fdd101483b8dfd93fb84bb3827fc5c7c7138740 Mon Sep 17 00:00:00 2001 From: "Kevin\" Seung Whan Chung" Date: Fri, 31 May 2024 07:18:57 -0700 Subject: [PATCH 10/11] add second to the time tag (#244) --- src/HDFrestart.cc | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/HDFrestart.cc b/src/HDFrestart.cc index 312de07c..f95292ac 100644 --- a/src/HDFrestart.cc +++ b/src/HDFrestart.cc @@ -161,6 +161,15 @@ void HDFrestart::addDateToFilename() filename_.append("0"); } filename_.append(extension); + + // seconds + filename_.append("_"); + extension = std::to_string(tt1->tm_sec); + if (tt1->tm_sec < 10) + { + filename_.append("0"); + } + filename_.append(extension); MGmol_MPI& mmpi = *(MGmol_MPI::instance()); // make sure all the PEs use the same filename std::vector name_buffer(filename_.begin(), filename_.end()); From 61276db032ca26ef451e208b85c7f16149c1f3a5 Mon Sep 17 00:00:00 2001 From: "Kevin\" Seung Whan Chung" Date: Mon, 3 Jun 2024 16:20:35 -0700 Subject: [PATCH 11/11] Hot-fix for `loadOrbitalsFromRestartFile` (#245) * delete orbital extrapolation after use. * removed extrapolation part --- src/md.cc | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/src/md.cc b/src/md.cc index 776886c4..c01589e5 100644 --- a/src/md.cc +++ b/src/md.cc @@ -722,18 +722,15 @@ OrbitalsType* MGmol::loadOrbitalFromRestartFile( then function is set as previous orbitals, while the extrapolated function is set as the current orbitals. This is how the restart file is saved via dumprestartFile. + + For now, we just enforce to not use the restart files with extrapolation. */ if (flag_extrapolated_data) { - orbitals_extrapol_ - = OrbitalsExtrapolationFactory::create( - ct.WFExtrapolation()); - - if (onpe0) os_ << "Create new orbitals_minus1..." << std::endl; + if (onpe0) + (*MPIdata::serr) << "loadRestartFile: does not support restart files with extrapolation." << std::endl; - orbitals_extrapol_->setupPreviousOrbitals(&restart_orbitals, - proj_matrices_, lrs_, local_cluster_, currentMasks_, corrMasks_, - h5file); + global_exit(0); } /* main workflow delete h5f_file_ here, meaning the loading is over. */