Fix bug: Reset BFGS history after cell changes

source/source_relax/ions_move_bfgs.cpp

@@ -37,6 +37,27 @@ void Ions_Move_BFGS::allocate()
     return;
 }
 
+void Ions_Move_BFGS::reset()
+{
+    if (init_done)
+    {
+        std::fill(pos.begin(), pos.end(), 0.0);
+        std::fill(pos_p.begin(), pos_p.end(), 0.0);
+        std::fill(grad.begin(), grad.end(), 0.0);
+        std::fill(grad_p.begin(), grad_p.end(), 0.0);
+        std::fill(move.begin(), move.end(), 0.0);
+        std::fill(move_p.begin(), move_p.end(), 0.0);
+
+        this->reset_hessian();
+    }
+    this->save_flag = false;
+    this->tr_min_hit = false;
+    this->first_step = true;
+
+    Ions_Move_Basic::trust_radius = 0.0;
+    Ions_Move_Basic::trust_radius_old = 0.0;
+}
+
 bool Ions_Move_BFGS::start(UnitCell& ucell, const ModuleBase::matrix& force, const double& energy_in, const int istep, int& update_iter, std::ofstream& ofs, std::vector<double>& etot_info)
 {
     ModuleBase::TITLE("Ions_Move_BFGS", "start");

source/source_relax/ions_move_bfgs.h

@@ -14,6 +14,7 @@ class Ions_Move_BFGS : public BFGS_Basic
     ~Ions_Move_BFGS();
 
     void allocate(void);
+    void reset(void);
     bool start(UnitCell& ucell, const ModuleBase::matrix& force, const double& energy_in, const int istep, int& update_iter, std::ofstream& ofs, std::vector<double>& etot_info);
 
   private:

source/source_relax/ions_move_bfgs2.cpp

@@ -33,6 +33,11 @@ void Ions_Move_BFGS2::allocate(const int _size)
     is_initialized = true;
 }
 
+void Ions_Move_BFGS2::reset()
+{
+    is_initialized = false;
+}
+
 
 bool Ions_Move_BFGS2::relax_step(const ModuleBase::matrix& _force,UnitCell& ucell, std::ofstream& ofs_running) 
 {

source/source_relax/ions_move_bfgs2.h

@@ -15,6 +15,7 @@ class Ions_Move_BFGS2
 {
 public:
     void allocate(const int _size);//initialize parameters
+    void reset();
     bool relax_step(const ModuleBase::matrix& _force,UnitCell& ucell, std::ofstream& ofs_running);//a full iteration step
 
 

source/source_relax/ions_move_methods.cpp

@@ -1,5 +1,7 @@
 #include "ions_move_methods.h"
 
+#include <algorithm>
+
 #include "ions_move_basic.h"
 #include "source_base/global_function.h"
 #include "source_base/global_variable.h"
@@ -93,3 +95,34 @@ void Ions_Move_Methods::cal_movement(const int &istep,
     }
     return;
 }
+
+void Ions_Move_Methods::reset_after_cell_change(const std::vector<std::string>& relax_method, std::ofstream& ofs)
+{
+    ModuleBase::TITLE("Ions_Move_Methods", "reset_after_cell_change");
+
+    if (relax_method.empty())
+    {
+        return;
+    }
+
+    const std::string method = relax_method[0];
+    const std::string method_arg = relax_method.size() > 1 ? relax_method[1] : "";
+    const auto reset_common_state = [this]() {
+        this->converged_ = false;
+        this->update_iter_ = 0;
+        std::fill(this->etot_info_.begin(), this->etot_info_.end(), 0.0);
+    };
+
+    if (method == "bfgs" && method_arg != "1")
+    {
+        reset_common_state();
+        this->bfgs.reset();
+        ofs << " Reset ionic BFGS history after cell change." << std::endl;
+    }
+    else if (method == "bfgs" && method_arg == "1")
+    {
+        reset_common_state();
+        this->bfgs_trad.reset();
+        ofs << " Reset traditional ionic BFGS history after cell change." << std::endl;
+    }
+}

source/source_relax/ions_move_methods.h

@@ -25,6 +25,7 @@ class Ions_Move_Methods
                       UnitCell &ucell,
                       std::ofstream& ofs,
                       std::vector<std::string>& relax_method);
+    void reset_after_cell_change(const std::vector<std::string>& relax_method, std::ofstream& ofs);
 
     bool get_converged() const
     {
@@ -56,6 +57,6 @@ class Ions_Move_Methods
     Ions_Move_LBFGS lbfgs;
     bool converged_ = false;
     int update_iter_ = 0;
-    std::vector<double> etot_info_{2, 0.0}; // [etot, etot_p]
+    std::vector<double> etot_info_{0.0, 0.0}; // [etot, etot_p]
 };
 #endif

source/source_relax/relax_nsync.cpp

@@ -128,6 +128,7 @@ bool IonCellOptimizer::relax_step(const int& istep,
             // Reset force_step counter after cell change for fresh atomic relaxation
             force_step = 1;
             stress_step++;
+            IMM.reset_after_cell_change(PARAM.inp.relax_method, ofs_running);
             ucell.cell_parameter_updated = true;
 
             // Update cell-related parameters after volume change
@@ -145,4 +146,3 @@ bool IonCellOptimizer::relax_step(const int& istep,
 
     return true;
 }
-

source/source_relax/test/CMakeLists.txt

@@ -86,6 +86,27 @@ AddTest(
           ${cell_source_files}
 )
 
+AddTest(
+  TARGET MODULE_RELAX_ions_move_methods_test
+  LIBS parameter ${math_libs} base device
+  SOURCES ions_move_methods_test.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
+          ../relax_data.cpp
+          ../../source_io/module_output/orb_io.cpp
+          ../../source_cell/bcast_cell.cpp
+          ../../source_cell/print_cell.cpp
+          ../../source_io/module_output/output.cpp
+)
+
 AddTest(
   TARGET MODULE_RELAX_ions_move_cg_test
   LIBS parameter ${math_libs} base device 
@@ -102,4 +123,4 @@ AddTest(
   TARGET MODULE_RELAX_ions_move_sd_test
   LIBS parameter ${math_libs} base device 
   SOURCES ions_move_sd_test.cpp ../ions_move_sd.cpp ../ions_move_basic.cpp ../relax_data.cpp ${cell_source_files}
-)
+)

source/source_relax/test/ions_move_methods_test.cpp

@@ -1,8 +1,16 @@
 #include "for_test.h"
 #include "gmock/gmock.h"
 #include "gtest/gtest.h"
+#include <algorithm>
+#include <cstdio>
+#include <fstream>
+#include <iterator>
+#include <string>
+#include <vector>
 #define private public
+#define protected public
 #include "source_relax/ions_move_methods.h"
+#undef protected
 #undef private
 /************************************************
  *  unit tests of class Ions_Move_Methods
@@ -62,6 +70,10 @@ namespace unitcell
         }
     }
 
+    void update_pos_tau(const Lattice&, const double*, const int, const int, Atom*)
+    {
+    }
+
     // Helper function to reset mock state
     void reset_remake_cell_mock()
     {
@@ -130,7 +142,9 @@ TEST_F(IonsMoveMethodsTest, CalMovement)
 {
     const int istep = 0;
     const int force_step = 1;
-    const ModuleBase::matrix f(3, 3);
+    ModuleBase::matrix f(natom, 3);
+    f(0, 0) = 0.1;
+    f(1, 1) = -0.1;
     const double etot = 0.0;
     UnitCell ucell;
     std::ofstream ofs;
@@ -205,7 +219,74 @@ TEST_F(IonsMoveMethodsTest, GetTrustRadius)
 // Test the get_update_iter() function
 TEST_F(IonsMoveMethodsTest, GetUpdateIter)
 {
-    Ions_Move_Basic::update_iter = 4;
+    imm.update_iter_ = 4;
 
     EXPECT_EQ(imm.get_update_iter(), 4);
-}
+}
+
+TEST_F(IonsMoveMethodsTest, ResetAfterCellChange)
+{
+    const std::string log_file = "reset_after_cell_change.log";
+    std::ofstream ofs(log_file);
+
+    imm.allocate(natom, "bfgs", "2");
+    imm.converged_ = true;
+    imm.update_iter_ = 4;
+    imm.etot_info_ = {-1.0, -2.0};
+    imm.bfgs.first_step = false;
+    imm.bfgs.save_flag = true;
+    imm.bfgs.tr_min_hit = true;
+    std::fill(imm.bfgs.pos.begin(), imm.bfgs.pos.end(), 1.0);
+    std::fill(imm.bfgs.pos_p.begin(), imm.bfgs.pos_p.end(), 2.0);
+    std::fill(imm.bfgs.grad.begin(), imm.bfgs.grad.end(), 3.0);
+    std::fill(imm.bfgs.grad_p.begin(), imm.bfgs.grad_p.end(), 4.0);
+    std::fill(imm.bfgs.move.begin(), imm.bfgs.move.end(), 5.0);
+    std::fill(imm.bfgs.move_p.begin(), imm.bfgs.move_p.end(), 6.0);
+    Ions_Move_Basic::trust_radius = 0.3;
+    Ions_Move_Basic::trust_radius_old = 0.2;
+
+    imm.reset_after_cell_change({"bfgs", "2"}, ofs);
+
+    EXPECT_FALSE(imm.converged_);
+    EXPECT_EQ(imm.update_iter_, 0);
+    EXPECT_THAT(imm.etot_info_, testing::Each(0.0));
+    EXPECT_TRUE(imm.bfgs.first_step);
+    EXPECT_FALSE(imm.bfgs.save_flag);
+    EXPECT_FALSE(imm.bfgs.tr_min_hit);
+    EXPECT_THAT(imm.bfgs.pos, testing::Each(0.0));
+    EXPECT_THAT(imm.bfgs.pos_p, testing::Each(0.0));
+    EXPECT_THAT(imm.bfgs.grad, testing::Each(0.0));
+    EXPECT_THAT(imm.bfgs.grad_p, testing::Each(0.0));
+    EXPECT_THAT(imm.bfgs.move, testing::Each(0.0));
+    EXPECT_THAT(imm.bfgs.move_p, testing::Each(0.0));
+    for (int i = 0; i < Ions_Move_Basic::dim; ++i)
+    {
+        for (int j = 0; j < Ions_Move_Basic::dim; ++j)
+        {
+            EXPECT_DOUBLE_EQ(imm.bfgs.inv_hess(i, j), i == j ? 1.0 : 0.0);
+        }
+    }
+    EXPECT_DOUBLE_EQ(Ions_Move_Basic::trust_radius, 0.0);
+    EXPECT_DOUBLE_EQ(Ions_Move_Basic::trust_radius_old, 0.0);
+
+    imm.allocate(natom, "bfgs", "1");
+    imm.converged_ = true;
+    imm.update_iter_ = 3;
+    imm.etot_info_ = {-3.0, -4.0};
+    ASSERT_TRUE(imm.bfgs_trad.is_initialized);
+
+    imm.reset_after_cell_change({"bfgs", "1"}, ofs);
+    ofs.close();
+
+    EXPECT_FALSE(imm.converged_);
+    EXPECT_EQ(imm.update_iter_, 0);
+    EXPECT_THAT(imm.etot_info_, testing::Each(0.0));
+    EXPECT_FALSE(imm.bfgs_trad.is_initialized);
+
+    std::ifstream ifs(log_file);
+    const std::string output((std::istreambuf_iterator<char>(ifs)), std::istreambuf_iterator<char>());
+    EXPECT_THAT(output, testing::HasSubstr("Reset ionic BFGS history after cell change."));
+    EXPECT_THAT(output, testing::HasSubstr("Reset traditional ionic BFGS history after cell change."));
+    ifs.close();
+    std::remove(log_file.c_str());
+}