diff --git a/source/source_main/driver_run.cpp b/source/source_main/driver_run.cpp index 67da741a04..ccfa6e5595 100644 --- a/source/source_main/driver_run.cpp +++ b/source/source_main/driver_run.cpp @@ -129,7 +129,10 @@ void Driver::driver_run() { ModuleBase::WARNING_QUIT("Driver::driver_run", "cannot recognize the 'calculation' command"); } - p_esolver->after_all_runners(ucell); + if (cal != "scf" && cal != "relax" && cal != "cell-relax" && cal != "nscf") + { + p_esolver->after_all_runners(ucell); + } } delete p_esolver; diff --git a/source/source_relax/ions_move_bfgs2.cpp b/source/source_relax/ions_move_bfgs2.cpp index 25d731bcea..15ec27b970 100644 --- a/source/source_relax/ions_move_bfgs2.cpp +++ b/source/source_relax/ions_move_bfgs2.cpp @@ -6,6 +6,7 @@ #include "source_cell/print_cell.h" // lanshuyue add 2025-06-19 #include "source_cell/update_cell.h" #include "source_io/module_parameter/parameter.h" +#include "source_relax/relax_criteria.h" //! initialize H0、H、pos0、force0、force void Ions_Move_BFGS2::allocate(const int _size) @@ -38,7 +39,10 @@ void Ions_Move_BFGS2::reset() is_initialized = false; } -bool Ions_Move_BFGS2::relax_step(const ModuleBase::matrix& _force, UnitCell& ucell, std::ofstream& ofs_running) +bool Ions_Move_BFGS2::relax_step(const ModuleBase::matrix& _force, + UnitCell& ucell, + std::ofstream& ofs_running, + const Relax_Criteria& criteria) { if (!is_initialized) { @@ -46,7 +50,6 @@ bool Ions_Move_BFGS2::relax_step(const ModuleBase::matrix& _force, UnitCell& uce } GetPos(ucell, pos); GetPostaud(ucell, pos_taud); - ucell.ionic_position_updated = true; assert(_force.nr == force.size() && _force.nc == 3); for (int i = 0; i < _force.nr; i++) { @@ -79,14 +82,19 @@ bool Ions_Move_BFGS2::relax_step(const ModuleBase::matrix& _force, UnitCell& uce k += ucell.atoms[i].na; } + this->CalculateLargestGrad(_force, ucell); + if (this->IsRestrain(criteria.force_thr_ev)) + { + return true; + } + this->PrepareStep(force, pos, H, pos0, force0, steplength, dpos, size, ucell); this->DetermineStep(steplength, dpos, maxstep); this->UpdatePos(ucell); - this->CalculateLargestGrad(_force, ucell); - bool converged = this->IsRestrain(); + ucell.ionic_position_updated = true; // print out geometry information during bfgs_trad relax unitcell::print_tau(ucell.atoms, ucell.Coordinate, ucell.ntype, ucell.lat0, ofs_running); - return converged; + return false; } void Ions_Move_BFGS2::GetPos(UnitCell& ucell, std::vector>& pos) @@ -370,7 +378,7 @@ void Ions_Move_BFGS2::CalculateLargestGrad(const ModuleBase::matrix& _force, Uni } } -bool Ions_Move_BFGS2::IsRestrain() +bool Ions_Move_BFGS2::IsRestrain(const double force_thr_ev) { - return Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A < PARAM.inp.force_thr_ev; + return Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A < force_thr_ev; } diff --git a/source/source_relax/ions_move_bfgs2.h b/source/source_relax/ions_move_bfgs2.h index 69d1c30883..37bab5a9aa 100644 --- a/source/source_relax/ions_move_bfgs2.h +++ b/source/source_relax/ions_move_bfgs2.h @@ -11,6 +11,8 @@ #include #include +struct Relax_Criteria; + class Ions_Move_BFGS2 { public: @@ -18,7 +20,8 @@ class Ions_Move_BFGS2 void reset(); bool relax_step(const ModuleBase::matrix& _force, UnitCell& ucell, - std::ofstream& ofs_running); // a full iteration step + std::ofstream& ofs_running, + const Relax_Criteria& criteria); // a full iteration step /// @brief whether allocate() has run and the optimiser holds live state. /// reset() clears it; ions_move_methods_test checks that. @@ -57,7 +60,7 @@ class Ions_Move_BFGS2 std::vector>& dpos, int& size, UnitCell& ucell); // calculate the atomic displacement in one iteration step - bool IsRestrain(); // check if converged + bool IsRestrain(const double force_thr_ev); // check if converged void CalculateLargestGrad(const ModuleBase::matrix& _force, UnitCell& ucell); void GetPos(UnitCell& ucell, std::vector>& pos); void GetPostaud(UnitCell& ucell, std::vector>& pos_taud); diff --git a/source/source_relax/ions_move_lbfgs.cpp b/source/source_relax/ions_move_lbfgs.cpp index 5fbb01c96c..07d99fd47b 100644 --- a/source/source_relax/ions_move_lbfgs.cpp +++ b/source/source_relax/ions_move_lbfgs.cpp @@ -36,8 +36,6 @@ bool Ions_Move_LBFGS::relax_step(const ModuleBase::matrix _force, { get_pos(ucell, pos); get_pos_taud(ucell, pos_taud); - // solver=p_esolver; - ucell.ionic_position_updated = true; for (int i = 0; i < _force.nr; i++) { for (int j = 0; j < _force.nc; j++) @@ -65,14 +63,18 @@ bool Ions_Move_LBFGS::relax_step(const ModuleBase::matrix _force, } k += ucell.atoms[i].na; } + this->calculate_largest_grad(_force, ucell); + if (this->is_restrain(criteria)) + { + return true; + } this->prepare_step(force, pos, H, pos0, force0, dpos, ucell, etot); this->determine_step(steplength, dpos, maxstep); this->update_pos(ucell); - this->calculate_largest_grad(_force, ucell); - bool converged = this->is_restrain(criteria); + ucell.ionic_position_updated = true; // mohan add 2025-06-22 unitcell::print_tau(ucell.atoms, ucell.Coordinate, ucell.ntype, ucell.lat0, ofs_running); - return converged; + return false; } void Ions_Move_LBFGS::get_pos(UnitCell& ucell, std::vector>& pos) diff --git a/source/source_relax/ions_move_methods.cpp b/source/source_relax/ions_move_methods.cpp index 6929921380..1eabddd924 100644 --- a/source/source_relax/ions_move_methods.cpp +++ b/source/source_relax/ions_move_methods.cpp @@ -83,7 +83,7 @@ void Ions_Move_Methods::cal_movement(const int& istep, } else if (relax_method == "bfgs_trad") { - converged_ = bfgs_trad.relax_step(f, ucell, ofs); + converged_ = bfgs_trad.relax_step(f, ucell, ofs, criteria); } else if (relax_method == "lbfgs") { diff --git a/source/source_relax/relax_driver.cpp b/source/source_relax/relax_driver.cpp index 7a0b124ee0..cf56ffd2da 100644 --- a/source/source_relax/relax_driver.cpp +++ b/source/source_relax/relax_driver.cpp @@ -29,6 +29,7 @@ void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver, { Socket_Driver socket_driver; socket_driver.socket_driver(p_esolver, ucell, inp, ofs_running); + p_esolver->after_all_runners(ucell); ModuleBase::timer::end("Relax_Driver", "relax_driver"); return; } @@ -51,6 +52,13 @@ void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver, // After relax_step() proposes a new geometry, it is not evaluated until // the next esolve() call; if we exit the loop early, force/stress are stale. bool geometry_evaluated = false; + RelaxStopReason stop_reason = RelaxStopReason::Finished; + bool after_all_runners_called = false; + + if (inp.relax_nmax == 0 && (inp.calculation == "relax" || inp.calculation == "cell-relax")) + { + stop_reason = RelaxStopReason::DryRun; + } while (steps[0] < inp.relax_nmax) { @@ -58,29 +66,59 @@ void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver, this->esolve(steps[0], p_esolver, ucell, inp, force, stress, etot); geometry_evaluated = true; this->stru_out(steps[0], ucell, inp, etot, stress, force); - bool converged = this->relax_step(steps, p_esolver, ucell, inp, force, stress, etot, ofs_running); - if (!converged) + this->json_out(p_esolver, ucell, inp, force, stress); + + if (ModuleIO::read_exit_file(my_rank_, "EXIT", ofs_running)) { - geometry_evaluated = false; + stop_reason = RelaxStopReason::UserExit; + break; } - this->json_out(p_esolver, ucell, inp, force, stress); + + const bool laststep = + (inp.calculation == "relax" || inp.calculation == "cell-relax") + && steps[0] + 1 == inp.relax_nmax; + if (laststep) + { + p_esolver->after_all_runners(ucell); + after_all_runners_called = true; + } + + const bool converged = this->relax_step(steps, p_esolver, ucell, inp, force, stress, etot, ofs_running); // Check stop conditions + if (ucell.ionic_position_updated || ucell.cell_parameter_updated) + { + geometry_evaluated = false; + } + if (converged) { - // Relaxation converged, exit loop immediately + if (inp.calculation == "relax" || inp.calculation == "cell-relax") + { + stop_reason = RelaxStopReason::Converged; + } + else + { + stop_reason = RelaxStopReason::Finished; + } break; } - else if (ModuleIO::read_exit_file(my_rank_, "EXIT", ofs_running)) + + if (steps[0] + 1 >= inp.relax_nmax) { - // EXIT file detected, exit loop + stop_reason = RelaxStopReason::MaxSteps; break; } ++steps[0]; } - this->final_out(steps[0], ucell, inp, etot, stress, force, geometry_evaluated, ofs_running); + this->final_out(steps[0], stop_reason, ucell, inp, etot, stress, force, geometry_evaluated, ofs_running); + + if (!after_all_runners_called) + { + p_esolver->after_all_runners(ucell); + } ModuleBase::timer::end("Relax_Driver", "relax_driver"); return; @@ -260,6 +298,7 @@ void Relax_Driver::json_out(ModuleESolver::ESolver* p_esolver, } void Relax_Driver::final_out(const int istep, + const RelaxStopReason stop_reason, UnitCell& ucell, const Input_para& inp, const double etot, @@ -268,6 +307,12 @@ void Relax_Driver::final_out(const int istep, const bool geometry_evaluated, std::ofstream& ofs_running) { + int ionic_steps = istep + 1; + if (stop_reason == RelaxStopReason::DryRun) + { + ionic_steps = 0; + } + // Structure final output is effective for scf/nscf/relax/cell-relax; // relax-specific screen messages remain guarded below. const bool is_relax = (inp.calculation == "relax" || inp.calculation == "cell-relax"); @@ -294,37 +339,37 @@ void Relax_Driver::final_out(const int istep, write_force); } - // relax_nmax == 0 is a valid dry-run mode: no relaxation step was ever - // taken, so it must not be reported as a failed relaxation. - if (istep == inp.relax_nmax && inp.relax_nmax > 0) + if (stop_reason == RelaxStopReason::MaxSteps) { - if (is_relax) - { - print_not_converged_summary(inp, ofs_running); + print_not_converged_summary(ionic_steps, inp, ofs_running); - std::cout << "\n ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; - std::cout << " Geometry relaxation stops here due to reaching the maximum " << std::endl; - std::cout << " relaxation steps. More steps are needed to converge the results " << std::endl; - std::cout << " ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; - } + std::cout << "\n ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; + std::cout << " Geometry relaxation stops here due to reaching the maximum " << std::endl; + std::cout << " relaxation steps. More steps are needed to converge the results " << std::endl; + std::cout << " ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; } - else if (inp.relax_nmax > 0) + else if (stop_reason == RelaxStopReason::Converged) { - if (is_relax) - { - std::cout << "\n ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; - std::cout << " Geometry relaxation thresholds are reached within " << istep << " steps." << std::endl; - std::cout << " ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; - } + std::cout << "\n ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; + std::cout << " Geometry relaxation thresholds are reached within " << ionic_steps << " steps." << std::endl; + std::cout << " ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; } - - if (is_relax && inp.relax_nmax == 0) + else if (stop_reason == RelaxStopReason::UserExit) + { + ofs_running << " Geometry relaxation stopped by EXIT after " + << ionic_steps << " ionic step(s)." << std::endl; + std::cout << " Geometry relaxation stopped by EXIT after " + << ionic_steps << " ionic step(s)." << std::endl; + } + else if (stop_reason == RelaxStopReason::DryRun) { print_dry_run_message(); } } -void Relax_Driver::print_not_converged_summary(const Input_para& inp, std::ofstream& ofs_running) const +void Relax_Driver::print_not_converged_summary(const int ionic_steps, + const Input_para& inp, + std::ofstream& ofs_running) const { // Unified not-converged summary for both relaxation paths so ASE and // users can read the method, the actual number of steps taken, and @@ -333,7 +378,6 @@ void Relax_Driver::print_not_converged_summary(const Input_para& inp, std::ofstr = inp.uses_simultaneous_relaxation() ? rl.get_max_force_history() : rl_old.get_max_force_history(); const std::vector& stress_hist = inp.uses_simultaneous_relaxation() ? rl.get_max_stress_history() : rl_old.get_max_stress_history(); - const int ionic_steps = static_cast(force_hist.size()); const std::string method = inp.relax_method.empty() ? "unknown" : inp.relax_method; ofs_running << " Relaxation method: " << method << std::endl; ofs_running << " Relaxation stopped after " << ionic_steps << " ionic step(s) (relax_nmax = " << inp.relax_nmax diff --git a/source/source_relax/relax_driver.h b/source/source_relax/relax_driver.h index 7050eba29b..34540f00b7 100644 --- a/source/source_relax/relax_driver.h +++ b/source/source_relax/relax_driver.h @@ -43,6 +43,15 @@ class Relax_Driver std::ofstream& ofs_running); private: + enum class RelaxStopReason + { + Converged, + MaxSteps, + UserExit, + DryRun, + Finished + }; + /// New relaxation optimizer (Relax class) Relax rl; /// Old relaxation optimizer (IonCellOptimizer class) @@ -148,6 +157,7 @@ class Relax_Driver * @brief Output final results after relaxation. * * @param istep Final iteration step. + * @param stop_reason Reason why the calculation stopped. * @param ucell Reference to the unit cell. * @param inp Input parameters for the calculation. * @param etot Total energy of the final step. @@ -158,6 +168,7 @@ class Relax_Driver * omitted / marked N/A. */ void final_out(const int istep, + const RelaxStopReason stop_reason, UnitCell& ucell, const Input_para& inp, const double etot, @@ -173,10 +184,13 @@ class Relax_Driver * per-step largest force / stress history, so ASE and users can diagnose * why the relaxation stopped at relax_nmax. * + * @param ionic_steps Number of evaluated ionic steps. * @param inp Input parameters (selects the active optimizer and relax_nmax). * @param ofs_running Output stream for the running log. */ - void print_not_converged_summary(const Input_para& inp, std::ofstream& ofs_running) const; + void print_not_converged_summary(const int ionic_steps, + const Input_para& inp, + std::ofstream& ofs_running) const; /** * @brief Print the dry-run success message to screen. diff --git a/source/source_relax/relax_nsync.cpp b/source/source_relax/relax_nsync.cpp index 0bd4732b29..345e800b65 100644 --- a/source/source_relax/relax_nsync.cpp +++ b/source/source_relax/relax_nsync.cpp @@ -54,14 +54,13 @@ void IonCellOptimizer::init_relax(const int& natom, const Input_para& inp) * @brief Perform one step of relaxation (atomic and/or cell). * * Main relaxation loop that coordinates atomic and cell relaxation: - * 1. Check for maximum iteration limit - * 2. Determine calculation mode (relax vs cell-relax) - * 3. Perform atomic relaxation if needed and atoms can move - * 4. If in cell-relax mode and atomic relaxation converged, perform cell relaxation + * 1. Determine calculation mode (relax vs cell-relax) + * 2. Perform atomic relaxation if needed and atoms can move + * 3. If in cell-relax mode and atomic relaxation converged, perform cell relaxation * * Convergence behavior: * - Returns false if relaxation is still in progress - * - Returns true if relaxation has converged or maximum iterations reached + * - Returns true if the current geometry has converged * * @param istep Current total iteration step * @param energy Total energy of the system @@ -95,27 +94,6 @@ bool IonCellOptimizer::relax_step(const int& istep, return true; } - // Check if we've reached the maximum number of iterations. - if (istep == inp_->relax_nmax) - { - // This step never ran cal_movement, so no force was recorded for it; - // the actual number of ionic steps taken is the history length. - const int ionic_steps = static_cast(max_force_history_.size()); - ofs_running << " Relaxation method: " << inp_->relax_method << std::endl; - ofs_running << " Relaxation stopped after " << ionic_steps - << " ionic step(s) (relax_nmax = " << inp_->relax_nmax << " reached)." << std::endl; - if (!max_force_history_.empty()) - { - ofs_running << " Largest force per step (eV/Angstrom):" << format_relax_history(max_force_history_); - } - if (!max_stress_history_.empty()) - { - ofs_running << " Largest stress per step (kbar):" << format_relax_history(max_stress_history_); - } - ofs_running << " Relaxation is not converged after reaching relax_nmax!" << std::endl; - return true; - } - // Determine calculation mode const bool is_cell_relax = (inp_->calculation == "cell-relax"); const bool is_relax = (inp_->calculation == "relax"); diff --git a/source/source_relax/relax_sync.cpp b/source/source_relax/relax_sync.cpp index 34df2b21cc..ba8deed576 100644 --- a/source/source_relax/relax_sync.cpp +++ b/source/source_relax/relax_sync.cpp @@ -68,6 +68,9 @@ bool Relax::relax_step(UnitCell& ucell, { ModuleBase::TITLE("Relax", "relax_step"); + ucell.ionic_position_updated = false; + ucell.cell_parameter_updated = false; + etot = etot_in * ModuleBase::Ry_to_eV; // convert to eV if (istep == 0) { diff --git a/source/source_relax/relax_sync.h b/source/source_relax/relax_sync.h index 8cdf70a27e..b43e5495bb 100644 --- a/source/source_relax/relax_sync.h +++ b/source/source_relax/relax_sync.h @@ -26,7 +26,6 @@ class Relax const ModuleBase::matrix& stress, const double etot_in, std::ofstream& ofs_running); - private: int istep = 0; // count ionic step @@ -170,4 +169,4 @@ class Relax } }; -#endif \ No newline at end of file +#endif diff --git a/source/source_relax/unittests/CMakeLists.txt b/source/source_relax/unittests/CMakeLists.txt index 30eaa2b1f7..e0490708bb 100644 --- a/source/source_relax/unittests/CMakeLists.txt +++ b/source/source_relax/unittests/CMakeLists.txt @@ -313,3 +313,32 @@ addtest( mock_remake_cell.cpp ../../source_cell/bcast_cell.cpp ../../source_cell/print_cell.cpp) + +addtest( + TARGET + MODULE_RELAX_relax_driver_test + SOURCES + test_relax_driver.cpp + ../relax_driver.cpp + ../relax_nsync.cpp + ../relax_sync.cpp + ../relax_data.cpp + ../line_search.cpp + ../ions_move_methods.cpp + ../ions_move_bfgs.cpp + ../ions_move_bfgs2.cpp + ../ions_move_cg.cpp + ../ions_move_sd.cpp + ../ions_move_lbfgs.cpp + ../ions_move_basic.cpp + ../bfgs_basic.cpp + ../cg_base.cpp + ../matrix_methods.cpp + ../lat_change_method.cpp + ../lattice_change_basic.cpp + ../lattice_change_cg.cpp + ${socket_json_sources} + LIBS + parameter + base + device) diff --git a/source/source_relax/unittests/test_ions_move_bfgs2.cpp b/source/source_relax/unittests/test_ions_move_bfgs2.cpp index c4f0022b37..4401215b37 100644 --- a/source/source_relax/unittests/test_ions_move_bfgs2.cpp +++ b/source/source_relax/unittests/test_ions_move_bfgs2.cpp @@ -4,6 +4,7 @@ #include "source_relax/ions_move_basic.h" // for Ions_Move_Basic static members #include "source_relax/ions_move_bfgs2.h" #include "source_relax/relax_data.h" +#include "source_relax/relax_criteria.h" #include "gmock/gmock.h" #include "gtest/gtest.h" @@ -183,10 +184,11 @@ TEST_F(BFGSTest, RelaxStepAutoInitialize) force(1, 2) = 0.0; std::ofstream ofs_running; + const Relax_Criteria criteria; // Before relax_step, is_initialized should be false EXPECT_FALSE(get_is_initialized()); - bfgs.relax_step(force, ucell, ofs_running); + bfgs.relax_step(force, ucell, ofs_running, criteria); // After relax_step, is_initialized should be true EXPECT_TRUE(get_is_initialized()); } @@ -387,8 +389,9 @@ TEST_F(BFGSTest, RelaxStepBasic) bfgs.allocate(ucell.nat); std::ofstream ofs_running; + const Relax_Criteria criteria; - bfgs.relax_step(force, ucell, ofs_running); + bfgs.relax_step(force, ucell, ofs_running, criteria); // Check that ionic_position_updated is true EXPECT_TRUE(ucell.ionic_position_updated); // Check that force values are set (converted units) @@ -398,3 +401,38 @@ TEST_F(BFGSTest, RelaxStepBasic) EXPECT_NEAR(get_pos()[0][0], 0.0, 1e-12); EXPECT_NE(get_pos()[1][0], 1.0); } + +TEST_F(BFGSTest, RelaxStepZeroForceReturnsBeforeStepConstruction) +{ + UnitCell ucell; + ucell.ntype = 1; + ucell.nat = 2; + ucell.lat0 = 1.0; + ucell.atoms = new Atom[ucell.ntype]; + ucell.atoms[0].na = 2; + ucell.atoms[0].tau = std::vector>(2); + ucell.atoms[0].taud = std::vector>(2); + ucell.atoms[0].mbl = std::vector>(2, {1, 1, 1}); + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); + ucell.iat2it[0] = 0; + ucell.iat2it[1] = 0; + ucell.iat2ia[0] = 0; + ucell.iat2ia[1] = 1; + ucell.atoms[0].tau[0].x = 0.25; + ucell.atoms[0].tau[1].x = 0.75; + + const double initial_tau0_x = ucell.atoms[0].tau[0].x; + const double initial_tau1_x = ucell.atoms[0].tau[1].x; + ModuleBase::matrix force(2, 3); + std::ofstream ofs_running; + Relax_Criteria criteria; + criteria.force_thr_ev = 1.0; + + const bool converged = bfgs.relax_step(force, ucell, ofs_running, criteria); + + EXPECT_TRUE(converged); + EXPECT_FALSE(ucell.ionic_position_updated); + EXPECT_NEAR(ucell.atoms[0].tau[0].x, initial_tau0_x, 1e-12); + EXPECT_NEAR(ucell.atoms[0].tau[1].x, initial_tau1_x, 1e-12); +} diff --git a/source/source_relax/unittests/test_ions_move_lbfgs.cpp b/source/source_relax/unittests/test_ions_move_lbfgs.cpp index 416f2d18a2..57c4ca0d2c 100644 --- a/source/source_relax/unittests/test_ions_move_lbfgs.cpp +++ b/source/source_relax/unittests/test_ions_move_lbfgs.cpp @@ -126,6 +126,7 @@ class IonsMoveLBFGSTest : public ::testing::Test TEST_F(IonsMoveLBFGSTest, RelaxStepConverged) { UnitCell ucell; + const double initial_tau_x = ucell.atoms[0].tau[0].x; ModuleBase::matrix force(natom, 3); force(0, 0) = 1.0e-4; // Ry/Bohr, ~0.005 eV/Angstrom const double etot = 0.0; @@ -140,6 +141,8 @@ TEST_F(IonsMoveLBFGSTest, RelaxStepConverged) EXPECT_TRUE(converged); EXPECT_NEAR(Ions_Move_Basic::largest_grad, 1.0e-4, 1e-12); // Ry/Bohr, lat0 = 1 EXPECT_LT(Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A, criteria.force_thr_ev); + EXPECT_NEAR(ucell.atoms[0].tau[0].x, initial_tau_x, 1e-12); + EXPECT_FALSE(ucell.ionic_position_updated); } // A force above force_thr_ev must be reported as not converged. diff --git a/source/source_relax/unittests/test_relax_driver.cpp b/source/source_relax/unittests/test_relax_driver.cpp new file mode 100644 index 0000000000..8aa08a6417 --- /dev/null +++ b/source/source_relax/unittests/test_relax_driver.cpp @@ -0,0 +1,414 @@ +#include "source_cell/cell_tools.h" +#include "source_cell/print_cell.h" +#include "source_cell/update_cell.h" +#include "source_esolver/esolver.h" +#include "for_test.h" +#include "source_relax/relax_driver.h" +#include "source_relax/relax_stru_io.h" +#include "source_relax/socket_driver.h" + +#include "gmock/gmock.h" +#include "gtest/gtest.h" + +#include +#include +#include +#include + +namespace unitcell +{ +bool if_atoms_can_move(const Atom*, const int) +{ + return true; +} + +bool if_cell_can_change(const std::vector&) +{ + return false; +} + +void update_pos_tau(const Lattice& lat, const double* pos, const int ntype, const int nat, Atom* atoms) +{ + int iat = 0; + for (int it = 0; it < ntype; ++it) + { + for (int ia = 0; ia < atoms[it].na; ++ia) + { + for (int ik = 0; ik < 3; ++ik) + { + atoms[it].tau[ia][ik] = pos[3 * iat + ik] / lat.lat0; + } + ++iat; + } + } + EXPECT_EQ(iat, nat); +} + +void update_pos_taud(const Lattice&, const double*, const int, const int, Atom*) +{ +} + +void setup_cell_after_vc(UnitCell&, std::ofstream&, const int) +{ +} + +void remake_cell(Lattice&) +{ +} + +void print_tau(Atom*, const std::string&, const int, const double, std::ofstream&) +{ +} +} // namespace unitcell + +namespace ModuleIO +{ +void print_screen(const int&, const int&, const int&) +{ +} + +void output_after_relax(const bool, const bool, std::ofstream&) +{ +} +} // namespace ModuleIO + +namespace relax_stru_io +{ +std::string build_stru_header(const int, + const double, + const ModuleBase::matrix&, + const Input_para&, + const bool, + const bool) +{ + return std::string(); +} + +bool need_orbital(const Input_para&) +{ + return false; +} + +void write_stru(UnitCell&, + const Input_para&, + const std::string&, + const std::string&, + const ModuleBase::matrix&, + const bool, + const bool, + const int, + const bool) +{ +} +} // namespace relax_stru_io + +void Socket_Driver::socket_driver(ModuleESolver::ESolver*, UnitCell&, const Input_para&, std::ofstream&) +{ +} + +namespace +{ + +class DriverTestESolver : public ModuleESolver::ESolver +{ + public: + explicit DriverTestESolver(const double force_x) : force_x_(force_x) + { + } + + void before_all_runners(BaseCell&, const Input_para&) override + { + } + + void runner(BaseCell&, const int) override + { + ++runner_count; + } + + void after_all_runners(BaseCell&) override + { + ++after_all_count; + } + + double cal_energy() override + { + return 0.0; + } + + void cal_force(BaseCell&, ModuleBase::matrix& force) override + { + force.zero_out(); + force(0, 0) = force_x_; + } + + void cal_stress(BaseCell&, ModuleBase::matrix& stress) override + { + stress.zero_out(); + } + + int runner_count = 0; + int after_all_count = 0; + + private: + double force_x_; +}; + +class RelaxDriverTest : public ::testing::Test +{ + protected: + void SetUp() override + { + log_file_ = testing::TempDir() + "running.log"; + } + + static void make_cell(UnitCell& ucell) + { + ucell.ntype = 1; + ucell.nat = 1; + ucell.lat0 = 1.0; + ucell.omega = 1.0; + ucell.lat.lat0 = 1.0; + ucell.lat.latvec.Identity(); + ucell.latvec.Identity(); + ucell.atoms = new Atom[1]; + ucell.atoms[0].na = 1; + ucell.atoms[0].tau.resize(1); + ucell.atoms[0].taud.resize(1); + ucell.atoms[0].mbl.resize(1, {1, 1, 1}); + ucell.iat2it.resize(1); + ucell.iat2ia.resize(1); + ucell.iat2it[0] = 0; + ucell.iat2ia[0] = 0; + } + + std::string read_running_log() const + { + std::ifstream ifs(log_file_); + return std::string((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); + } + + std::string log_file_; +}; + +TEST_F(RelaxDriverTest, RelaxNmaxZeroDoesNotRunESolver) +{ + Input_para inp; + inp.calculation = "relax"; + inp.relax_method = "lbfgs"; + inp.relax_nmax = 0; + inp.out_level = "m"; + + UnitCell ucell; + make_cell(ucell); + DriverTestESolver esolver(1.0); + std::ofstream ofs(log_file_); + Relax_Driver driver; + testing::internal::CaptureStdout(); + driver.relax_driver(&esolver, ucell, inp, ofs); + const std::string output = testing::internal::GetCapturedStdout(); + ofs.close(); + + EXPECT_EQ(esolver.runner_count, 0); + EXPECT_EQ(esolver.after_all_count, 1); + EXPECT_THAT(output, testing::HasSubstr("relax_nmax = 0, DRY RUN TEST SUCCEEDS")); +} + +TEST_F(RelaxDriverTest, FirstStepConvergesOnNonSimultaneousPath) +{ + Input_para inp; + inp.calculation = "relax"; + inp.relax_method = "lbfgs"; + inp.relax_nmax = 3; + inp.force_thr = 0.1; + inp.force_thr_ev = 0.1 * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A; + inp.cal_force = 1; + inp.out_level = "m"; + + UnitCell ucell; + make_cell(ucell); + DriverTestESolver esolver(1.0e-4); + std::ofstream ofs(log_file_); + Relax_Driver driver; + testing::internal::CaptureStdout(); + driver.relax_driver(&esolver, ucell, inp, ofs); + const std::string output = testing::internal::GetCapturedStdout(); + ofs.close(); + + EXPECT_EQ(esolver.runner_count, 1); + EXPECT_EQ(esolver.after_all_count, 1); + EXPECT_THAT(read_running_log(), testing::HasSubstr("Relaxation is converged!")); + EXPECT_THAT(output, testing::HasSubstr("Geometry relaxation thresholds are reached within 1 steps.")); +} + +TEST_F(RelaxDriverTest, FirstStepConvergesOnSimultaneousPath) +{ + Input_para inp; + inp.calculation = "relax"; + inp.relax_method = "cg"; + inp.relax_nmax = 3; + inp.force_thr = 0.1; + inp.force_thr_ev = 0.1 * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A; + inp.cal_force = 1; + inp.out_level = "m"; + + UnitCell ucell; + make_cell(ucell); + DriverTestESolver esolver(1.0e-4); + std::ofstream ofs(log_file_); + Relax_Driver driver; + testing::internal::CaptureStdout(); + driver.relax_driver(&esolver, ucell, inp, ofs); + const std::string output = testing::internal::GetCapturedStdout(); + ofs.close(); + + EXPECT_EQ(esolver.runner_count, 1); + EXPECT_EQ(esolver.after_all_count, 1); + EXPECT_THAT(read_running_log(), testing::HasSubstr("Relaxation is converged!")); + EXPECT_THAT(output, testing::HasSubstr("Geometry relaxation thresholds are reached within 1 steps.")); +} + +TEST_F(RelaxDriverTest, MaxStepsMovesFinalNonSimultaneousGeometry) +{ + Input_para inp; + inp.calculation = "relax"; + inp.relax_method = "lbfgs"; + inp.relax_nmax = 1; + inp.force_thr = 0.0; + inp.force_thr_ev = 0.0; + inp.cal_force = 1; + inp.relax_bfgs_rmax = 0.2; + inp.out_level = "m"; + + UnitCell ucell; + make_cell(ucell); + DriverTestESolver esolver(1.0); + std::ofstream ofs(log_file_); + Relax_Driver driver; + driver.relax_driver(&esolver, ucell, inp, ofs); + ofs.close(); + + EXPECT_EQ(esolver.runner_count, 1); + EXPECT_EQ(esolver.after_all_count, 1); + EXPECT_THAT(read_running_log(), testing::HasSubstr("reaching relax_nmax")); + EXPECT_THAT(read_running_log(), testing::HasSubstr("stopped after 1 ionic step(s) (relax_nmax = 1 reached)")); +} + +TEST_F(RelaxDriverTest, MaxStepsStopsOnSimultaneousPath) +{ + Input_para inp; + inp.calculation = "relax"; + inp.relax_method = "cg"; + inp.relax_nmax = 1; + inp.force_thr = -1.0; + inp.force_thr_ev = -1.0; + inp.cal_force = 1; + inp.out_level = "m"; + + UnitCell ucell; + make_cell(ucell); + DriverTestESolver esolver(1.0); + std::ofstream ofs(log_file_); + Relax_Driver driver; + driver.relax_driver(&esolver, ucell, inp, ofs); + ofs.close(); + + EXPECT_EQ(esolver.runner_count, 1); + EXPECT_EQ(esolver.after_all_count, 1); + EXPECT_THAT(read_running_log(), testing::HasSubstr("reaching relax_nmax")); + EXPECT_THAT(read_running_log(), testing::HasSubstr("stopped after 1 ionic step(s) (relax_nmax = 1 reached)")); +} + +TEST_F(RelaxDriverTest, FinalAllowedStepConvergesAndCallsAfterAllOnce) +{ + Input_para inp; + inp.calculation = "relax"; + inp.relax_method = "lbfgs"; + inp.relax_nmax = 1; + inp.force_thr = 0.1; + inp.force_thr_ev = 0.1 * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A; + inp.cal_force = 1; + inp.out_level = "m"; + + UnitCell ucell; + make_cell(ucell); + DriverTestESolver esolver(1.0e-4); + std::ofstream ofs(log_file_); + Relax_Driver driver; + driver.relax_driver(&esolver, ucell, inp, ofs); + ofs.close(); + + EXPECT_EQ(esolver.runner_count, 1); + EXPECT_EQ(esolver.after_all_count, 1); + EXPECT_THAT(read_running_log(), testing::HasSubstr("Relaxation is converged!")); +} + +TEST_F(RelaxDriverTest, ExitStopsBeforeNonSimultaneousRelaxStep) +{ + char old_cwd[4096]; + ASSERT_NE(::getcwd(old_cwd, sizeof(old_cwd)), nullptr); + ASSERT_EQ(::chdir(testing::TempDir().c_str()), 0); + { + std::ofstream exit_file("EXIT"); + exit_file << "stop_ion true\n"; + } + + Input_para inp; + inp.calculation = "relax"; + inp.relax_method = "lbfgs"; + inp.relax_nmax = 3; + inp.force_thr = 0.0; + inp.force_thr_ev = 0.0; + inp.cal_force = 1; + inp.out_level = "m"; + + UnitCell ucell; + make_cell(ucell); + DriverTestESolver esolver(1.0); + std::ofstream ofs(log_file_); + Relax_Driver driver; + driver.relax_driver(&esolver, ucell, inp, ofs); + ofs.close(); + + ASSERT_EQ(::chdir(old_cwd), 0); + EXPECT_EQ(esolver.runner_count, 1); + EXPECT_EQ(esolver.after_all_count, 1); + EXPECT_DOUBLE_EQ(ucell.atoms[0].tau[0][0], 0.0); + EXPECT_THAT(read_running_log(), testing::HasSubstr("stopped by EXIT")); + EXPECT_THAT(read_running_log(), testing::HasSubstr("stopped by EXIT after 1 ionic step(s)")); +} + +TEST_F(RelaxDriverTest, ExitStopsBeforeSimultaneousRelaxStep) +{ + char old_cwd[4096]; + ASSERT_NE(::getcwd(old_cwd, sizeof(old_cwd)), nullptr); + ASSERT_EQ(::chdir(testing::TempDir().c_str()), 0); + { + std::ofstream exit_file("EXIT"); + exit_file << "stop_ion true\n"; + } + + Input_para inp; + inp.calculation = "relax"; + inp.relax_method = "cg"; + inp.relax_nmax = 3; + inp.force_thr = -1.0; + inp.force_thr_ev = -1.0; + inp.cal_force = 1; + inp.out_level = "m"; + + UnitCell ucell; + make_cell(ucell); + DriverTestESolver esolver(1.0); + std::ofstream ofs(log_file_); + Relax_Driver driver; + driver.relax_driver(&esolver, ucell, inp, ofs); + ofs.close(); + + ASSERT_EQ(::chdir(old_cwd), 0); + EXPECT_EQ(esolver.runner_count, 1); + EXPECT_EQ(esolver.after_all_count, 1); + EXPECT_THAT(read_running_log(), testing::HasSubstr("stopped by EXIT")); + EXPECT_THAT(read_running_log(), testing::HasSubstr("stopped by EXIT after 1 ionic step(s)")); +} + +} // namespace diff --git a/source/source_relax/unittests/test_relax_nsync.cpp b/source/source_relax/unittests/test_relax_nsync.cpp index 7351c5eccf..9005d2f66c 100644 --- a/source/source_relax/unittests/test_relax_nsync.cpp +++ b/source/source_relax/unittests/test_relax_nsync.cpp @@ -1,6 +1,7 @@ #include "source_base/global_variable.h" #include "source_cell/unitcell.h" #include "source_io/module_parameter/parameter.h" +#include "source_relax/ions_move_bfgs2.h" #include "source_relax/lattice_change_basic.h" #include "source_relax/relax_nsync.h" @@ -23,7 +24,6 @@ * Largest stress is ... kbar while threshold is ... kbar * Relaxation is converged! * Relaxation is not converged yet! - * Relaxation is not converged after reaching relax_nmax! * The markers are what the ASE interface greps for (issue #6546). */ @@ -227,6 +227,27 @@ TEST_F(IonCellOptimizerTest, RelaxConverged) EXPECT_THAT(log, testing::HasSubstr(" Largest force per step (eV/Angstrom):")); EXPECT_THAT(log, testing::HasSubstr("\n Relaxation is converged!")); EXPECT_THAT(log, testing::Not(testing::HasSubstr("not converged"))); + EXPECT_FALSE(ucell.ionic_position_updated); +} + +TEST_F(IonCellOptimizerTest, BfgsTradConvergedBeforeMove) +{ + inp.relax_method = "bfgs_trad"; + inp.relax_nmax = 2; + UnitCell ucell; + setup_ucell(ucell); + optimizer.init_relax(natom, inp); + ModuleBase::matrix force(natom, 3); + ModuleBase::matrix stress(3, 3); + int force_step = 1; + int stress_step = 1; + std::ofstream ofs(log_file); + + const bool done = optimizer.relax_step(1, etot, ucell, force, stress, force_step, stress_step, ofs); + ofs.close(); + + EXPECT_TRUE(done); + EXPECT_FALSE(ucell.ionic_position_updated); } // relax + large force: not converged, and no ASE converged marker may appear. @@ -252,13 +273,16 @@ TEST_F(IonCellOptimizerTest, RelaxNotConverged) EXPECT_THAT(log, testing::Not(testing::HasSubstr("Relaxation is converged!"))); } -// Reaching relax_nmax must terminate with an explicit not-converged marker. -TEST_F(IonCellOptimizerTest, RelaxNmaxReached) +// The last allowed ionic step still performs its optimizer move before the +// driver reports that relax_nmax was reached. +TEST_F(IonCellOptimizerTest, RelaxNmaxLastStepConverged) { UnitCell ucell; setup_ucell(ucell); + inp.relax_nmax = 1; optimizer.init_relax(natom, inp); ModuleBase::matrix force(natom, 3); + force(0, 0) = 1.0e-4; ModuleBase::matrix stress(3, 3); int force_step = 1; int stress_step = 1; @@ -268,15 +292,35 @@ TEST_F(IonCellOptimizerTest, RelaxNmaxReached) ofs.close(); EXPECT_TRUE(done); - // Key invariant: when relax_nmax terminates the loop, cal_movement is - // skipped and the geometry is NOT touched. The force/stress captured by - // the driver before relax_step are still consistent with ucell, so the - // driver's geometry_evaluated=true is honest. - EXPECT_FALSE(ucell.ionic_position_updated); const std::string log = read_log(); EXPECT_THAT(log, testing::HasSubstr(" Relaxation method: lbfgs")); - EXPECT_THAT(log, testing::HasSubstr(" ionic step(s) (relax_nmax = 50 reached).")); - EXPECT_THAT(log, testing::HasSubstr(" Relaxation is not converged after reaching relax_nmax!")); + EXPECT_THAT(log, testing::HasSubstr(" Relaxation converged in 1 step(s).")); + EXPECT_THAT(log, testing::HasSubstr("\n Relaxation is converged!")); +} + +// A non-converged last allowed step must return not-converged to the driver; +// the driver, rather than IonCellOptimizer, owns the relax_nmax summary. +TEST_F(IonCellOptimizerTest, RelaxNmaxLastStepNotConverged) +{ + UnitCell ucell; + setup_ucell(ucell); + inp.relax_nmax = 1; + optimizer.init_relax(natom, inp); + ModuleBase::matrix force(natom, 3); + force(0, 0) = 1.0; + ModuleBase::matrix stress(3, 3); + int force_step = 1; + int stress_step = 1; + std::ofstream ofs(log_file); + + const bool done = optimizer.relax_step(inp.relax_nmax, etot, ucell, force, stress, force_step, stress_step, ofs); + ofs.close(); + + EXPECT_FALSE(done); + EXPECT_TRUE(ucell.ionic_position_updated); + const std::string log = read_log(); + EXPECT_THAT(log, testing::HasSubstr(" Relaxation is not converged yet!")); + EXPECT_THAT(log, testing::Not(testing::HasSubstr(" relax_nmax = 1 reached."))); EXPECT_THAT(log, testing::Not(testing::HasSubstr("Relaxation is converged!"))); } @@ -315,6 +359,7 @@ TEST_F(IonCellOptimizerTest, CellRelaxConverged) TEST_F(IonCellOptimizerTest, CellRelaxNotConverged) { inp.calculation = "cell-relax"; + inp.relax_nmax = 1; UnitCell ucell; setup_ucell(ucell); optimizer.init_relax(natom, inp); @@ -330,6 +375,8 @@ TEST_F(IonCellOptimizerTest, CellRelaxNotConverged) ofs.close(); EXPECT_FALSE(done); + EXPECT_FALSE(ucell.ionic_position_updated); + EXPECT_TRUE(ucell.cell_parameter_updated); const std::string log = read_log(); EXPECT_THAT(log, testing::HasSubstr(" Largest stress is ")); EXPECT_THAT(log, testing::HasSubstr("\n Relaxation is not converged yet!")); diff --git a/source/source_relax/unittests/test_relax_sync.cpp b/source/source_relax/unittests/test_relax_sync.cpp index dc27dd0157..65cb65c91f 100644 --- a/source/source_relax/unittests/test_relax_sync.cpp +++ b/source/source_relax/unittests/test_relax_sync.cpp @@ -5,6 +5,7 @@ #include "gmock/gmock.h" #include "gtest/gtest.h" +#include #include #include @@ -780,8 +781,7 @@ TEST(RelaxSyncCellMove, FixedShapeIsotropicScale) // --------------------------------------------------------------------------- // Max-step termination on the Relax (simultaneous) path. // -// Unlike IonCellOptimizer::relax_step (which checks istep == relax_nmax and -// returns true without moving), Relax::relax_step has no max-step guard. +// Relax::relax_step has no max-step guard. // When the driver loop hits relax_nmax, the last relax_step still calls // move_cell_ions, so the geometry IS updated and the force/stress in the // driver become stale. The driver then sets geometry_evaluated=false. @@ -841,6 +841,62 @@ TEST(RelaxSyncMaxStep, GeometryMovedAndNotConverged) delete[] ucell.atoms; } +TEST(RelaxSyncConvergence, ClearsPreviousGeometryUpdate) +{ + const int nat = 1; + Input_para inp; + inp.calculation = "relax"; + inp.relax_method = "cg"; + inp.force_thr = 0.01; + inp.force_thr_ev = 1.0; + inp.fixed_axes = "None"; + inp.fixed_ibrav = false; + + UnitCell ucell; + ucell.ntype = 1; + ucell.nat = nat; + ucell.atoms = new Atom[1]; + ucell.atoms[0].na = nat; + ucell.atoms[0].label = "Si"; + ucell.omega = 1.0; + ucell.lat0 = 1.0; + ucell.iat2it.resize(nat); + ucell.iat2ia.resize(nat); + ucell.iat2it[0] = 0; + ucell.iat2ia[0] = 0; + ucell.atoms[0].mbl.resize(nat); + ucell.atoms[0].taud.resize(nat); + ucell.atoms[0].tau.resize(nat); + ucell.atoms[0].dis.resize(nat); + ucell.atoms[0].mag.resize(nat); + ucell.atoms[0].vel.resize(nat); + ucell.atoms[0].mbl[0] = {1, 1, 1}; + ucell.latvec.Identity(); + + ModuleBase::matrix moving_force(nat, 3); + moving_force(0, 0) = 0.1; + ModuleBase::matrix converged_force(nat, 3); + ModuleBase::matrix stress(3, 3); + + Relax rl; + rl.init_relax(nat, inp); + const std::string log_file = testing::TempDir() + "running_relax_sync_flags.log"; + std::ofstream ofs(log_file); + const bool first_done = rl.relax_step(ucell, moving_force, stress, 0.0, ofs); + EXPECT_FALSE(first_done); + EXPECT_TRUE(ucell.ionic_position_updated); + + const bool second_done = rl.relax_step(ucell, converged_force, stress, 0.0, ofs); + ofs.close(); + std::remove(log_file.c_str()); + + EXPECT_TRUE(second_done); + EXPECT_FALSE(ucell.ionic_position_updated); + EXPECT_FALSE(ucell.cell_parameter_updated); + + delete[] ucell.atoms; +} + // fixed_axes="volume": the cell may change shape but its volume is re-scaled // back to the original omega after each move (relax_sync.cpp Step 1). TEST(RelaxSyncCellMove, FixedVolumePreservesOmega) @@ -973,4 +1029,4 @@ TEST(RelaxSyncCellMove, MixedMblMovesOnlyFreeComponents) EXPECT_NE(ucell.atoms[1].taud[1].z, 0.0); free_two_type_cell(ucell); -} \ No newline at end of file +}