Add initialization checking in BFGS, and replace the old BFGS with a new one in CG-BFGS.

source/source_relax/bfgs.cpp

@@ -30,11 +30,16 @@ void BFGS::allocate(const int _size)
     force0 = std::vector<double>(3*size, 0.0);
     force = std::vector<ModuleBase::Vector3<double>>(size, ModuleBase::Vector3<double>(0.0, 0.0, 0.0));
     steplength = std::vector<double>(size, 0.0);  
+    is_initialized = true;
 }
 
 
 void BFGS::relax_step(const ModuleBase::matrix& _force,UnitCell& ucell) 
 {
+    if(!is_initialized)
+    {
+        allocate(ucell.nat);
+    }
     GetPos(ucell,pos);  
     GetPostaud(ucell,pos_taud);
     ucell.ionic_position_updated = true;

source/source_relax/bfgs.h

@@ -24,6 +24,7 @@ class BFGS
     double maxstep;//every movement smaller than maxstep
     double largest_grad;
     int size;//number of atoms
+    bool is_initialized=false;
 
     std::vector<double> steplength;//the length of atoms displacement 
     std::vector<std::vector<double>> H;//Hessian matrix

source/source_relax/ions_move_cg.cpp

@@ -2,6 +2,7 @@
 #include "ions_move_basic.h"
 #include "source_base/global_function.h"
 #include "source_base/global_variable.h"
+
 using namespace Ions_Move_Basic;
 
 double Ions_Move_CG::RELAX_CG_THR = -1.0; // default is 0.5
@@ -169,7 +170,7 @@ void Ions_Move_CG::start(UnitCell &ucell, const ModuleBase::matrix &force, const
                     < RELAX_CG_THR) // cg to bfgs  by pengfei 13-8-8
                 {
                     Ions_Move_Basic::relax_method[0] = "bfgs";
-                    Ions_Move_Basic::relax_method[1] = "2";
+                    Ions_Move_Basic::relax_method[1] = "1";
                 }
                 Ions_Move_Basic::best_xxx = steplength;
             }

source/source_relax/test/bfgs_test.cpp

@@ -54,6 +54,34 @@ TEST_F(BFGSTest, TestAllocate)
     EXPECT_EQ(bfgs.largest_grad,0.0);
 }
 
+// Test that relax_step will auto-initialize if not already initialized
+TEST_F(BFGSTest, RelaxStepAutoInitialize)
+{
+    bfgs.is_initialized = false;
+
+    UnitCell ucell;
+    ucell.nat = 2;
+    ucell.ntype = 1;
+    ucell.atoms = new Atom[ucell.ntype];
+    ucell.atoms[0].na = 2;
+    ucell.atoms[0].tau = std::vector<ModuleBase::Vector3<double>>(2);
+    ucell.atoms[0].taud = std::vector<ModuleBase::Vector3<double>>(2);
+    ucell.atoms[0].mbl = std::vector<ModuleBase::Vector3<int>>(2, {1, 1, 1});
+    ucell.atoms[0].tau[0].x = 0.0; ucell.atoms[0].tau[0].y = 0.0; ucell.atoms[0].tau[0].z = 0.0;
+    ucell.atoms[0].tau[1].x = 1.0; ucell.atoms[0].tau[1].y = 0.0; ucell.atoms[0].tau[1].z = 0.0;
+    ucell.lat0 = 1.0;
+
+    ModuleBase::matrix force(2, 3);
+    force(0, 0) = 0.1; force(0, 1) = 0.0; force(0, 2) = 0.0;
+    force(1, 0) = -0.1; force(1, 1) = 0.0; force(1, 2) = 0.0;
+
+    // Before relax_step, is_initialized should be false
+    EXPECT_FALSE(bfgs.is_initialized);
+    bfgs.relax_step(force, ucell);
+    // After relax_step, is_initialized should be true
+    EXPECT_TRUE(bfgs.is_initialized);
+}
+
 // Test if a dimension less than or equal to 0 results in an assertion error
 TEST_F(BFGSTest, TestAllocateWithZeroDimension)
 {
@@ -103,10 +131,8 @@ TEST_F(BFGSTest, GetPosAndPostaud)
     ucell.atoms[0].mbl = std::vector<ModuleBase::Vector3<int>>(2, {1, 1, 1});
 
     // set coordinates
-    ucell.atoms[0].tau[0].x = 1.0; ucell.atoms[0].tau[0].y = 2.0; ucell.atoms[0].tau[0].z = 3.0;
-    ucell.atoms[0].tau[1].x = 2.0; ucell.atoms[0].tau[1].y = 3.0; ucell.atoms[0].tau[1].z = 4.0;
-    ucell.atoms[0].taud[0].x = 0.1; ucell.atoms[0].taud[0].y = 0.2; ucell.atoms[0].taud[0].z = 0.3;
-    ucell.atoms[0].taud[1].x = 0.4; ucell.atoms[0].taud[1].y = 0.5; ucell.atoms[0].taud[1].z = 0.6;
+    ucell.atoms[0].tau[0].x = 0.0; ucell.atoms[0].tau[0].y = 0.0; ucell.atoms[0].tau[0].z = 0.0;
+    ucell.atoms[0].tau[1].x = 1.0; ucell.atoms[0].tau[1].y = 0.0; ucell.atoms[0].tau[1].z = 0.0;
 
     // allocate mapping arrays 
     ucell.iat2it = new int[ucell.nat];