diff --git a/examples/PinnedH2O/job.basis_1_50 b/examples/PinnedH2O/job.basis_1_50 index bf0a4c49..5291a7b7 100644 --- a/examples/PinnedH2O/job.basis_1_50 +++ b/examples/PinnedH2O/job.basis_1_50 @@ -10,7 +10,7 @@ setenv OMP_NUM_THREADS 1 set ncpus = 8 -set maindir = /p/lustre2/cheung26/mgmol-20250219 +set maindir = /p/lustre2/cheung26/mgmol setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH diff --git a/examples/PinnedH2O/job.offline b/examples/PinnedH2O/job.offline index 64997e09..9d7150f8 100644 --- a/examples/PinnedH2O/job.offline +++ b/examples/PinnedH2O/job.offline @@ -10,7 +10,7 @@ setenv OMP_NUM_THREADS 1 set ncpus = 8 -set maindir = /p/lustre2/cheung26/mgmol-20250219 +set maindir = /p/lustre2/cheung26/mgmol setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH diff --git a/examples/PinnedH2O/job.ref b/examples/PinnedH2O/job.ref index 2c4660cf..32024e6f 100644 --- a/examples/PinnedH2O/job.ref +++ b/examples/PinnedH2O/job.ref @@ -11,7 +11,7 @@ setenv OMP_NUM_THREADS 1 set ncpus = 8 set case = 2 -set maindir = /p/lustre2/cheung26/mgmol-20250219 +set maindir = /p/lustre2/cheung26/mgmol setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH diff --git a/examples/PinnedH2O/job.rom b/examples/PinnedH2O/job.rom index fc4ae1b0..45285136 100644 --- a/examples/PinnedH2O/job.rom +++ b/examples/PinnedH2O/job.rom @@ -10,7 +10,7 @@ setenv OMP_NUM_THREADS 1 set ncpus = 8 -set maindir = /p/lustre2/cheung26/mgmol-20250219 +set maindir = /p/lustre2/cheung26/mgmol setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH diff --git a/examples/PinnedH2O/job.rom_3DOF b/examples/PinnedH2O/job.rom_3DOF index 8b328199..0e38688b 100644 --- a/examples/PinnedH2O/job.rom_3DOF +++ b/examples/PinnedH2O/job.rom_3DOF @@ -11,7 +11,7 @@ setenv OMP_NUM_THREADS 1 set ncpus = 8 set case = 2 -set maindir = /p/lustre2/cheung26/mgmol-20250219 +set maindir = /p/lustre2/cheung26/mgmol setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH diff --git a/examples/PinnedH2O/mgmol_rom_3DOF_test1.cfg b/examples/PinnedH2O/mgmol_rom_3DOF_test1.cfg index adfdfc43..b7a78a0d 100644 --- a/examples/PinnedH2O/mgmol_rom_3DOF_test1.cfg +++ b/examples/PinnedH2O/mgmol_rom_3DOF_test1.cfg @@ -1,4 +1,4 @@ -verbosity=1 +verbosity=2 xcFunctional=PBE FDtype=4th [Mesh] @@ -34,7 +34,9 @@ initial_width=2. nempty=30 [DensityMatrix] solver=MVP -nb_inner_it=100 +nb_inner_it=20 +mixing=0.8 +tol=1.e-8 [Restart] output_level=4 [ROM] diff --git a/examples/PinnedH2O/mgmol_rom_3DOF_test2.cfg b/examples/PinnedH2O/mgmol_rom_3DOF_test2.cfg index f0368790..24273c40 100644 --- a/examples/PinnedH2O/mgmol_rom_3DOF_test2.cfg +++ b/examples/PinnedH2O/mgmol_rom_3DOF_test2.cfg @@ -1,4 +1,4 @@ -verbosity=1 +verbosity=2 xcFunctional=PBE FDtype=4th [Mesh] @@ -30,7 +30,9 @@ initial_width=2. nempty=30 [DensityMatrix] solver=MVP -nb_inner_it=100 +nb_inner_it=20 +mixing=0.8 +tol=1.e-8 [Restart] output_level=4 [ROM] diff --git a/src/Control.cc b/src/Control.cc index 151fe920..cbf2eadd 100644 --- a/src/Control.cc +++ b/src/Control.cc @@ -224,6 +224,7 @@ 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; if (DMEigensolver() == DMEigensolverType::Eigensolver) { os << " Density matrix computation algorithm = " @@ -439,7 +440,7 @@ void Control::sync(void) memset(&int_buffer[0], 0, size_int_buffer * sizeof(int)); } - const short size_float_buffer = 43; + const short size_float_buffer = 44; float* float_buffer = new float[size_float_buffer]; if (mype_ == 0) { @@ -485,6 +486,7 @@ void Control::sync(void) float_buffer[40] = threshold_eigenvalue_gram_quench_; float_buffer[41] = pair_mlwf_distance_threshold_; float_buffer[42] = e0_; + float_buffer[43] = dm_tol; } else { @@ -680,6 +682,7 @@ void Control::sync(void) threshold_eigenvalue_gram_quench_ = float_buffer[40]; pair_mlwf_distance_threshold_ = float_buffer[41]; e0_ = float_buffer[42]; + dm_tol = float_buffer[43]; max_electronic_steps_loose_ = max_electronic_steps; delete[] short_buffer; @@ -699,8 +702,8 @@ void Control::setDefaultValues() void Control::adjust() { - // change dm_mix default to 1. if not using Davidson - if (it_algo_type_ != 2 && dm_mix < 0.) dm_mix = 1.; + // change dm_mix default to 1. if not using Davidson or MVP + if ((it_algo_type_ != 2 && DM_solver_ != 1) && dm_mix < 0.) dm_mix = 1.; if (nel_ - 2 * numst == 0) { @@ -1720,6 +1723,7 @@ void Control::setOptions(const boost::program_options::variables_map& vm) lrs_extrapolation = 10; dm_mix = vm["DensityMatrix.mixing"].as(); + dm_tol = vm["DensityMatrix.tol"].as(); dm_inner_steps = vm["DensityMatrix.nb_inner_it"].as(); dm_use_old_ = vm["DensityMatrix.use_old"].as() ? 1 : 0; str = vm["DensityMatrix.algo"].as(); @@ -1739,8 +1743,6 @@ void Control::setOptions(const boost::program_options::variables_map& vm) else dm_algo_ = 2; - dm_tol = vm["DensityMatrix.tol"].as(); - str = vm["DensityMatrix.solver"].as(); if (str.compare("Mixing") == 0) DM_solver_ = 0; if (str.compare("MVP") == 0) DM_solver_ = 1; diff --git a/src/Ions.cc b/src/Ions.cc index 072f0950..1879e1ba 100644 --- a/src/Ions.cc +++ b/src/Ions.cc @@ -1315,22 +1315,6 @@ void Ions::setLocalForces( { if (ion->compareName(*s)) { - //std::cout << "Ion found: " << ion->name() << std::endl; - //std::cout << "Ion force: (" << *it << ", " << *(it + 1) << ", " << *(it + 2) << ")" << std::endl; - //std::cout << "names: "; - //for (int i = 0; i < names.size(); i++) - //{ - // std::cout << names[i]; - // if (i == forces.size() - 1) std::cout << std::endl; - // else std::cout << ", "; - //} - //std::cout << "forces: "; - //for (int i = 0; i < forces.size(); i++) - //{ - // std::cout << forces[i]; - // if (i == forces.size() - 1) std::cout << ")" << std::endl; - // else std::cout << ", "; - //} ion->set_force(0, *it); ion->set_force(1, *(it + 1)); ion->set_force(2, *(it + 2)); @@ -1364,9 +1348,6 @@ void Ions::setLocalForces( double d = std::sqrt(d2); if (d < tol) { - std::cout << "Ion found: " << ion->name() << std::endl; - std::cout << "Ion position:( " << *cit << ", " << *(cit + 1) << ", " << *(cit + 2) << ")" << std::endl; - std::cout << "Ion force: (" << *fit << ", " << *(fit + 1) << ", " << *(fit + 2) << ")" << std::endl; ion->set_force(0, *fit); ion->set_force(1, *(fit + 1)); ion->set_force(2, *(fit + 2)); diff --git a/src/MGmol.cc b/src/MGmol.cc index facd362e..ec780d0e 100644 --- a/src/MGmol.cc +++ b/src/MGmol.cc @@ -102,6 +102,7 @@ extern Timer md_moveVnuc_tm; extern Timer md_updateMasks_tm; extern Timer md_extrapolateOrbitals_tm; extern Timer md_updateRhoAndPot_tm; +extern Timer md_updateDMandEnergy_tm; extern Timer quench_tm; extern Timer ions_setupInteractingIons_tm; extern Timer ions_setup_tm; @@ -928,6 +929,7 @@ void MGmol::printTimers() init_nuc_tm_.print(os_); md_updateMasks_tm.print(os_); md_extrapolateOrbitals_tm.print(os_); + md_updateDMandEnergy_tm.print(os_); quench_tm.print(os_); evnl_tm_.print(os_); ions_setupInteractingIons_tm.print(os_); diff --git a/src/MVPSolver.cc b/src/MVPSolver.cc index 6990e29d..1c093e53 100644 --- a/src/MVPSolver.cc +++ b/src/MVPSolver.cc @@ -47,13 +47,16 @@ MVPSolver::MVPSolver(MPI_Comm comm, std::ostream& os, Electrostatic* electrostat, MGmol* mgmol_strategy, const int numst, const double kbT, const std::vector>& global_indexes, - const short n_inner_steps, const bool use_old_dm) + const short n_inner_steps, const double mixing, const double tol_de0, + const bool use_old_dm) : comm_(comm), os_(os), n_inner_steps_(n_inner_steps), use_old_dm_(use_old_dm), ions_(ions), - numst_(numst) + numst_(numst), + mixing_(mixing), + tol_de0_(tol_de0) { Control& ct = *(Control::instance()); if (onpe0 && ct.verbose > 0) @@ -207,7 +210,6 @@ int MVPSolver::solve(OrbitalsType& orbitals) kbpsi.computeHvnlMatrix(&kbpsi, ions_, h11_nl); - const double tol_de0 = 1.e-12; for (int inner_it = 0; inner_it < n_inner_steps_; inner_it++) { if (onpe0 && ct.verbose > 1) @@ -270,7 +272,8 @@ int MVPSolver::solve(OrbitalsType& orbitals) double de0 = evaluateDerivative(dmInit, delta_dm, ts0); - if (std::abs(de0) < tol_de0 && inner_it > 0) + // check for convergence + if (std::abs(de0) < tol_de0_ && inner_it > 0) { if (onpe0 && ct.verbose > 0) std::cout << "MVP: de0 = " << de0 @@ -278,54 +281,70 @@ int MVPSolver::solve(OrbitalsType& orbitals) break; } - // - // evaluate free energy at beta=1 - // - if (onpe0 && ct.verbose > 2) - std::cout << "MVP --- Target energy..." << std::endl; - proj_mat_work_->setDM(target, orbitals.getIterativeIndex()); - proj_mat_work_->computeOccupationsFromDM(); - if (ct.verbose > 2) proj_mat_work_->printOccupations(os_); - const double nel = proj_mat_work_->getNel(); - if (onpe0 && ct.verbose > 1) - os_ << "MVP --- Number of electrons at beta=1 : " << nel - << std::endl; - - rho_->computeRho(orbitals, target); - - mgmol_strategy_->update_pot(vh_init, ions_); - - energy_->saveVofRho(); - - // update h11 + double beta = 0.; + if (mixing_ > 0.) { - h11 = h11_nl; - mgmol_strategy_->addHlocal2matrix(orbitals, orbitals, h11); + beta = mixing_; + if (onpe0 && ct.verbose > 0) + { + if (ct.verbose > 1) + os_ << "MVP with beta = " << beta << std::endl; + os_ << std::setprecision(12); + os_ << std::fixed << "MVP inner iteration " << inner_it + << ", E0=" << e0 << std::endl; + } } - - proj_mat_work_->assignH(h11); - proj_mat_work_->setHB2H(); - - const double ts1 - = evalEntropyMVP(proj_mat_work_, (ct.verbose > 2), os_); - const double e1 = energy_->evaluateTotal( - ts1, proj_mat_work_, ions_, orbitals, ct.verbose - 1, os_); - - // line minimization - const double beta - = minQuadPolynomial(e0, e1, de0, (ct.verbose > 2), os_); - assert(!std::isnan(beta)); - - if (onpe0 && ct.verbose > 0) + else { - os_ << std::setprecision(12); - os_ << std::fixed << "MVP inner iteration " << inner_it - << ", E0=" << e0 << ", E1=" << e1; - os_ << std::scientific << ", E0'=" << de0 - << " -> beta=" << beta; - os_ << std::endl; - } + // + // evaluate free energy at beta=1 + // + if (onpe0 && ct.verbose > 2) + std::cout << "MVP --- Target energy..." << std::endl; + proj_mat_work_->setDM(target, orbitals.getIterativeIndex()); + proj_mat_work_->computeOccupationsFromDM(); + if (ct.verbose > 2) proj_mat_work_->printOccupations(os_); + const double nel = proj_mat_work_->getNel(); + if (onpe0 && ct.verbose > 1) + os_ << "MVP --- Number of electrons at beta=1 : " << nel + << std::endl; + + rho_->computeRho(orbitals, target); + + mgmol_strategy_->update_pot(vh_init, ions_); + + energy_->saveVofRho(); + + // update h11 + { + h11 = h11_nl; + mgmol_strategy_->addHlocal2matrix( + orbitals, orbitals, h11); + } + + proj_mat_work_->assignH(h11); + proj_mat_work_->setHB2H(); + + const double ts1 + = evalEntropyMVP(proj_mat_work_, (ct.verbose > 2), os_); + const double e1 = energy_->evaluateTotal(ts1, + proj_mat_work_, ions_, orbitals, ct.verbose - 1, os_); + + // line minimization + beta + = minQuadPolynomial(e0, e1, de0, (ct.verbose > 2), os_); + assert(!std::isnan(beta)); + if (onpe0 && ct.verbose > 0) + { + os_ << std::setprecision(12); + os_ << std::fixed << "MVP inner iteration " << inner_it + << ", E0=" << e0 << ", E1=" << e1; + os_ << std::scientific << ", E0'=" << de0 + << " -> beta=" << beta; + os_ << std::endl; + } + } // update DM *work_ = dmInit; work_->axpy(beta, delta_dm); diff --git a/src/MVPSolver.h b/src/MVPSolver.h index 8e7a7c70..77a6a345 100644 --- a/src/MVPSolver.h +++ b/src/MVPSolver.h @@ -33,6 +33,12 @@ class MVPSolver Ions& ions_; int numst_; + double mixing_; + + /*! + * tolerance on energy slope in inner iterations + */ + double tol_de0_; Rho* rho_; Energy* energy_; @@ -56,7 +62,8 @@ class MVPSolver Electrostatic* electrostat, MGmol* mgmol_strategy, const int numst, const double kbT, const std::vector>& global_indexes, - const short n_inner_steps, const bool use_old_dm); + const short n_inner_steps, const double mixing, const double tol_de0, + const bool use_old_dm); ~MVPSolver(); int solve(OrbitalsType& orbitals); diff --git a/src/MVP_DMStrategy.cc b/src/MVP_DMStrategy.cc index 8712ae92..3fb53dee 100644 --- a/src/MVP_DMStrategy.cc +++ b/src/MVP_DMStrategy.cc @@ -54,7 +54,7 @@ int MVP_DMStrategy::update(OrbitalsType& orbitals) MVPSolver solver(comm_, os_, ions_, rho_, energy_, electrostat_, mgmol_strategy_, ct.numst, ct.occ_width, global_indexes_, - ct.dm_inner_steps, use_old_dm_); + ct.dm_inner_steps, ct.dm_mix, ct.dm_tol, use_old_dm_); return solver.solve(orbitals); } diff --git a/src/md.cc b/src/md.cc index bc986955..e8d06133 100644 --- a/src/md.cc +++ b/src/md.cc @@ -50,6 +50,7 @@ Timer md_tau_tm("md_tau"); Timer md_moveVnuc_tm("md_moveVnuc"); Timer md_updateMasks_tm("md_updateMasks"); Timer md_extrapolateOrbitals_tm("md_extrapolateOrbitals"); +Timer md_updateDMandEnergy_tm("md_updateDMandEnergy"); #define DUMP_MAX_NUM_TRY 5 @@ -485,7 +486,9 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) if (ROM_MVP) { + md_updateDMandEnergy_tm.start(); updateDMandEnergy(**orbitals, *ROM_ions, eks); + md_updateDMandEnergy_tm.stop(); } else { diff --git a/src/read_config.cc b/src/read_config.cc index 884aa8d3..9dbd56d3 100644 --- a/src/read_config.cc +++ b/src/read_config.cc @@ -302,7 +302,9 @@ int read_config(int argc, char** argv, po::variables_map& vm, po::value()->default_value(0), "Flag for computing new centers from extrapolated orbitals.")( "DensityMatrix.mixing", po::value()->default_value(-1.), - "Mixing coefficient for Density Matrix")("DensityMatrix.solver", + "Mixing coefficient for Density Matrix")("DensityMatrix.tol", + po::value()->default_value(1.e-12), + "Tolerance for Density Matrix convergence")("DensityMatrix.solver", po::value()->default_value("Mixing"), "Algorithm for updating Density Matrix: Mixing, MVP, HMVP")( "DensityMatrix.nb_inner_it", po::value()->default_value(3), @@ -324,10 +326,7 @@ int read_config(int argc, char** argv, po::variables_map& vm, po::value()->default_value(100), "Maximum number of iterations for power method " "to compute interval for Chebyshev " - "approximation of density matrix. ")("DensityMatrix.tol", - po::value()->default_value(1.e-7), - "tolerance, used in iterative DM computation convergence " - "criteria"); + "approximation of density matrix. "); po::options_description cmdline_options; cmdline_options.add(generic); diff --git a/src/rom_workflows.cc b/src/rom_workflows.cc index 8fbd4cb3..f8d70eec 100644 --- a/src/rom_workflows.cc +++ b/src/rom_workflows.cc @@ -841,6 +841,8 @@ void testROMIonDensity(MGmolInterface *mgmol_) mmpi.bcastGlobal(cfgs[idx].data(), 3 * num_ions, 0); } + /* Artificial ions object to avoid repeated setPositions */ + Ions* new_ions; /* Collect fictitious ion density based on each configuration */ std::vector> fom_rhoc(num_snap); /* Sanity check for overlappingVL_ions */ @@ -848,24 +850,25 @@ void testROMIonDensity(MGmolInterface *mgmol_) for (int idx = 0; idx < num_snap; idx++) { /* set ion positions */ - ions->setPositions(cfgs[idx], atnumbers); + new_ions = new Ions(cfgs[idx], atnumbers, lattice, ions->getSpecies()); /* save overlapping ions for sanity check */ - fom_overlap_ions[idx].resize(ions->overlappingVL_ions().size()); - for (int k = 0; k < ions->overlappingVL_ions().size(); k++) + fom_overlap_ions[idx].resize(new_ions->overlappingVL_ions().size()); + for (int k = 0; k < new_ions->overlappingVL_ions().size(); k++) { fom_overlap_ions[idx][k].resize(3); for (int d = 0; d < 3; d++) - fom_overlap_ions[idx][k][d] = ions->overlappingVL_ions()[k]->position(d); + fom_overlap_ions[idx][k][d] = new_ions->overlappingVL_ions()[k]->position(d); } /* compute resulting ion density */ /* NOTE: we exclude rescaling for the sake of verification */ - pot.initialize(*ions); + mgmol->setupPotentials(*new_ions); - mgmol->electrostat_->setupRhoc(pot.rho_comp()); fom_rhoc[idx].resize(dim); mgmol->electrostat_->getRhoc()->init_vect(fom_rhoc[idx].data(), 'd'); + + delete new_ions; } /* Initialize libROM classes */ @@ -911,23 +914,32 @@ void testROMIonDensity(MGmolInterface *mgmol_) if (rank == 0) printf("test index: %d\n", test_idx); /* set ion positions */ - ions->setPositions(cfgs[test_idx], atnumbers); + new_ions = new Ions(cfgs[test_idx], atnumbers, lattice, ions->getSpecies()); /* Sanity check for overlapping ions */ - CAROM_VERIFY(fom_overlap_ions[test_idx].size() == ions->overlappingVL_ions().size()); - for (int k = 0; k < ions->overlappingVL_ions().size(); k++) + CAROM_VERIFY(fom_overlap_ions[test_idx].size() == new_ions->overlappingVL_ions().size()); + for (int k = 0; k < new_ions->overlappingVL_ions().size(); k++) for (int d = 0; d < 3; d++) - CAROM_VERIFY(abs(fom_overlap_ions[test_idx][k][d] - ions->overlappingVL_ions()[k]->position(d)) < 1.0e-12); + CAROM_VERIFY(abs(fom_overlap_ions[test_idx][k][d] - new_ions->overlappingVL_ions()[k]->position(d)) < 1.0e-12); + + /* set up potentials */ + mgmol->setupPotentials(*new_ions); /* eval ion density on sample grid points */ std::vector sampled_rhoc(sampled_row.size()); - pot.evalIonDensityOnSamplePts(*ions, sampled_row, sampled_rhoc); + pot.evalIonDensityOnSamplePts(*new_ions, sampled_row, sampled_rhoc); + // For now, we relax the threshold to allow the slight difference of the values at the sampled indices. + // This is because the rescaling procedure after getting the radial data on mesh requires global information, + // thus cannot be done in the ROM level only with sampled values. + // After we have the DEIM reconstruction, we will do rescaling and revisit this. for (int d = 0; d < sampled_row.size(); d++) { printf("rank %d, fom rhoc[%d]: %.3e, rom rhoc: %.3e\n", rank, sampled_row[d], fom_rhoc[test_idx][sampled_row[d]], sampled_rhoc[d]); - CAROM_VERIFY(abs(fom_rhoc[test_idx][sampled_row[d]] - sampled_rhoc[d]) < 1.0e-12); + CAROM_VERIFY(abs(fom_rhoc[test_idx][sampled_row[d]] - sampled_rhoc[d]) < 1.0e-2); } + + delete new_ions; } template void readRestartFiles(MGmolInterface *mgmol_); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index b4844557..f55ce699 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -511,6 +511,13 @@ add_test(NAME testMVP ${CMAKE_CURRENT_SOURCE_DIR}/MVP/mvp.cfg ${CMAKE_CURRENT_SOURCE_DIR}/MVP/coords.in ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) +add_test(NAME testMVPmix + COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/MVPmix/test.py + ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} + ${CMAKE_CURRENT_BINARY_DIR}/../src/mgmol-opt + ${CMAKE_CURRENT_SOURCE_DIR}/MVPmix/mgmol.cfg + ${CMAKE_CURRENT_SOURCE_DIR}/MVPmix/h2o.xyz + ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) add_test(NAME BandGapN2 COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/BandGapN2/test.py ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 3 ${MPIEXEC_PREFLAGS} diff --git a/tests/MVP/test.py b/tests/MVP/test.py index bd5d708e..57d9522e 100644 --- a/tests/MVP/test.py +++ b/tests/MVP/test.py @@ -51,6 +51,8 @@ energies=[] print("Check forces are smaller than tol = {}".format(tol)) for line in lines: + if line.count(b'MVP'): + print(line) if line.count(b'%%'): print(line) words=line.split() diff --git a/tests/MVPmix/h2o.xyz b/tests/MVPmix/h2o.xyz new file mode 100644 index 00000000..d5171c8b --- /dev/null +++ b/tests/MVPmix/h2o.xyz @@ -0,0 +1,6 @@ +3 + +O 0.00 0.00 0.00 +H -0.76 0.59 0.00 +H 0.76 0.59 0.00 + diff --git a/tests/MVPmix/mgmol.cfg b/tests/MVPmix/mgmol.cfg new file mode 100644 index 00000000..10edc464 --- /dev/null +++ b/tests/MVPmix/mgmol.cfg @@ -0,0 +1,36 @@ +verbosity=2 +xcFunctional=PBE +FDtype=4th +[Mesh] +nx=48 +ny=48 +nz=48 +[Domain] +ox=-4.5 +oy=-4.5 +oz=-4.5 +lx=9. +ly=9. +lz=9. +[Potentials] +pseudopotential=pseudo.O_ONCV_PBE_SG15 +pseudopotential=pseudo.H_ONCV_PBE_SG15 +[Run] +type=QUENCH +[Quench] +solver=PSD +max_steps=100 +atol=1.e-8 +step_length=2. +ortho_freq=10 +[Orbitals] +initial_type=random +initial_width=1.5 +temperature=10. +nempty=5 +[Restart] +output_level=0 +[DensityMatrix] +solver=MVP +nb_inner_it=2 +mixing=0.5 diff --git a/tests/MVPmix/test.py b/tests/MVPmix/test.py new file mode 100644 index 00000000..8ffafdbb --- /dev/null +++ b/tests/MVPmix/test.py @@ -0,0 +1,103 @@ +#!/usr/bin/env python +import sys +import os +import subprocess +import string + +print("Test MVP solver with mixing coefficient...") + +nargs=len(sys.argv) + +mpicmd = sys.argv[1]+" "+sys.argv[2]+" "+sys.argv[3] +for i in range(4,nargs-4): + mpicmd = mpicmd + " "+sys.argv[i] +print("MPI run command: {}".format(mpicmd)) + +exe = sys.argv[nargs-4] +inp = sys.argv[nargs-3] +coords = sys.argv[nargs-2] +print("coordinates file: %s"%coords) + +#create links to potentials files +cwd = os.getcwd() + +dst = 'pseudo.O_ONCV_PBE_SG15' +src = sys.argv[nargs-1] + '/' + dst +if not os.path.exists(cwd+'/'+dst): + print("Create link to %s"%dst) + os.symlink(src, dst) + +dst = 'pseudo.H_ONCV_PBE_SG15' +src = sys.argv[nargs-1] + '/' + dst +if not os.path.exists(cwd+'/'+dst): + print("Create link to %s"%dst) + os.symlink(src, dst) + +#run mgmol +command = "{} {} -c {} -i {}".format(mpicmd,exe,inp,coords) +print("Run command: {}".format(command)) + +output = subprocess.check_output(command,stderr=subprocess.STDOUT,shell=True) + +#analyse mgmol standard output +#make sure force is below tolerance +lines=output.split(b'\n') + +convergence=0 +for line in lines: + if line.count(b'DFTsolver:') and line.count(b'convergence'): + convergence=1 + break + +if convergence==0: + print("MVP Solver did not converge") + sys.exit(1) + +flag = 0 +eigenvalues=[] +energies=[] +ecount=0 +for line in lines: + if line.count(b'FERMI'): + flag = 0 + if flag==1: + words=line.split() + for w in words: + eigenvalues.append(eval(w)) + if line.count(b'Eigenvalues'): + flag = 1 + eigenvalues=[] + if line.count(b'%%'): + words=line.split() + e=words[5][0:-1] + print(e) + ecount=ecount+1 + energies.append(eval(e)) +print(energies) + +print(eigenvalues) +tol = 1.e-4 +eigenvalue0 = -0.916 +if abs(eigenvalues[0]-eigenvalue0)>tol: + print("Expected eigenvalue 0 to be {}".format(eigenvalue0)) + sys.exit(1) +eigenvalue8 = 0.219 +if abs(eigenvalues[8]-eigenvalue8)>tol: + print("Expected eigenvalue 8 to be {}".format(eigenvalue8)) + sys.exit(1) + +niterations = ecount +print("MVP solver ran for {} iterations".format(niterations)) +if niterations>180: + print("MVP test FAILED for taking too many iterations") + sys.exit(1) + +print("Check energy...") +last_energy = energies[-1] +print("Energy = {}".format(last_energy)) +if last_energy>-17.16269: + print("Last energy = {}".format(last_energy)) + sys.exit(1) + +print("Test PASSED") +sys.exit(0) diff --git a/tests/PinnedH2O_3DOF/job.basis b/tests/PinnedH2O_3DOF/job.basis index cab929d1..35f564a3 100644 --- a/tests/PinnedH2O_3DOF/job.basis +++ b/tests/PinnedH2O_3DOF/job.basis @@ -10,7 +10,7 @@ setenv OMP_NUM_THREADS 1 set ncpus = 8 -set maindir = /p/lustre2/cheung26/mgmol-20250219 +set maindir = /p/lustre2/cheung26/mgmol setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH diff --git a/tests/PinnedH2O_3DOF/job.offline b/tests/PinnedH2O_3DOF/job.offline index 243f563a..5b4224f0 100644 --- a/tests/PinnedH2O_3DOF/job.offline +++ b/tests/PinnedH2O_3DOF/job.offline @@ -10,7 +10,7 @@ setenv OMP_NUM_THREADS 1 set ncpus = 8 -set maindir = /p/lustre2/cheung26/mgmol-20250219 +set maindir = /p/lustre2/cheung26/mgmol setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH diff --git a/tests/PinnedH2O_3DOF/job.online b/tests/PinnedH2O_3DOF/job.online index 1373e011..ea75e59f 100644 --- a/tests/PinnedH2O_3DOF/job.online +++ b/tests/PinnedH2O_3DOF/job.online @@ -10,7 +10,7 @@ setenv OMP_NUM_THREADS 1 set ncpus = 8 -set maindir = /p/lustre2/cheung26/mgmol-20250219 +set maindir = /p/lustre2/cheung26/mgmol setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH diff --git a/tests/PinnedH2O_3DOF/job.online_rotate b/tests/PinnedH2O_3DOF/job.online_rotate index 3f8a65ad..a2065b72 100644 --- a/tests/PinnedH2O_3DOF/job.online_rotate +++ b/tests/PinnedH2O_3DOF/job.online_rotate @@ -10,7 +10,7 @@ setenv OMP_NUM_THREADS 1 set ncpus = 8 -set maindir = /p/lustre2/cheung26/mgmol-20250219 +set maindir = /p/lustre2/cheung26/mgmol setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH diff --git a/tests/WFEnergyAndForces/testWFEnergyAndForces.cc b/tests/WFEnergyAndForces/testWFEnergyAndForces.cc index 64765c7a..f9039abe 100644 --- a/tests/WFEnergyAndForces/testWFEnergyAndForces.cc +++ b/tests/WFEnergyAndForces/testWFEnergyAndForces.cc @@ -190,6 +190,9 @@ int main(int argc, char** argv) } } + // rename output restart file + ct.out_restart_file = ct.out_restart_file + "1"; + mgmol->cleanup(); delete mgmol;