diff --git a/source/source_relax/bfgs.cpp b/source/source_relax/bfgs.cpp index 829af3fa1c..00f5de2cc4 100644 --- a/source/source_relax/bfgs.cpp +++ b/source/source_relax/bfgs.cpp @@ -10,10 +10,12 @@ //! initialize H0、H、pos0、force0、force void BFGS::allocate(const int _size) { + assert(_size > 0); alpha=70;//default value in ase is 70 maxstep=PARAM.inp.relax_bfgs_rmax; size=_size; - sign =true; + largest_grad=0.0; + sign=true; H = std::vector>(3*size, std::vector(3*size, 0.0)); for (int i = 0; i < 3*size; ++i) @@ -37,6 +39,7 @@ void BFGS::relax_step(const ModuleBase::matrix& _force,UnitCell& ucell) GetPos(ucell,pos); GetPostaud(ucell,pos_taud); ucell.ionic_position_updated = true; + assert(_force.nr == force.size() && _force.nc == force[0].size()); for(int i = 0; i < _force.nr; i++) { for(int j=0;j<_force.nc;j++) @@ -65,7 +68,7 @@ void BFGS::relax_step(const ModuleBase::matrix& _force,UnitCell& ucell) k+=ucell.atoms[i].na; } - this->PrepareStep(force,pos,H,pos0,force0,steplength,dpos,ucell); + this->PrepareStep(force,pos,H,pos0,force0,steplength,dpos,size,ucell); this->DetermineStep(steplength,dpos,maxstep); this->UpdatePos(ucell); this->CalculateLargestGrad(_force,ucell); @@ -76,6 +79,7 @@ void BFGS::relax_step(const ModuleBase::matrix& _force,UnitCell& ucell) void BFGS::GetPos(UnitCell& ucell,std::vector>& pos) { + assert(pos.size() == ucell.nat); int k=0; for(int i=0;i>& pos) void BFGS::GetPostaud(UnitCell& ucell, std::vector>& pos_taud) { + assert(pos_taud.size() == ucell.nat); int k=0; for(int i=0;i>& force, std::vector& force0, std::vector& steplength, std::vector>& dpos, + int& size, UnitCell& ucell) { std::vector changedforce = ReshapeMToV(force); @@ -142,8 +148,10 @@ void BFGS::PrepareStep(std::vector>& force, } } std::vector a=DotInMAndV2(V, changedforce); + double threshold=1e-8; for(int i = 0; i < a.size(); i++) { + assert(std::abs(omega[i]) > threshold); a[i]/=std::abs(omega[i]); } std::vector tmpdpos = DotInMAndV1(V, a); @@ -221,7 +229,8 @@ void BFGS::Update(std::vector& pos, dpos[iat * 3 + 2] = move_ion_dr.z ; } } - if(*max_element(dpos.begin(), dpos.end()) < 1e-7) + double threshold=1e-7; + if(*max_element(dpos.begin(), dpos.end()) < threshold) { return; } @@ -243,6 +252,8 @@ void BFGS::DetermineStep(std::vector& steplength, { std::vector::iterator maxsteplength = max_element(steplength.begin(), steplength.end()); double a = *maxsteplength; + double threshold=1e-10; + assert(a > threshold); if(a >= maxstep) { double scale = maxstep / a; @@ -311,8 +322,8 @@ void BFGS::UpdatePos(UnitCell& ucell) void BFGS::IsRestrain(std::vector>& dpos) { - Ions_Move_Basic::converged = Ions_Move_Basic::largest_grad - * ModuleBase::Ry_to_eV / 0.529177 steplength;//the length of atoms displacement std::vector> H;//Hessian matrix std::vector force0;//force in previous step @@ -22,17 +33,8 @@ class BFGS std::vector pos_taud0;//atom pos in previous step(relative coordinates) std::vector> pos_taud; std::vector> dpos; - - void allocate(const int _size);//initialize parameters - void relax_step(const ModuleBase::matrix& _force,UnitCell& ucell);//a full iteration step - void PrepareStep(std::vector>& force,std::vector>& pos,std::vector>& H,std::vector& pos0,std::vector& force0,std::vector& steplength,std::vector>& dpos,UnitCell& ucell);//calculate the atomic displacement in one iteration step - -private: - bool sign;//check if this is the first iteration - double alpha;//initialize H,diagonal element is alpha - double maxstep;//every movement smaller than maxstep - int size;//number of atoms + void PrepareStep(std::vector>& force,std::vector>& pos,std::vector>& H,std::vector& pos0,std::vector& force0,std::vector& steplength,std::vector>& dpos,int& size,UnitCell& ucell);//calculate the atomic displacement in one iteration step void IsRestrain(std::vector>& dpos);//check if converged void CalculateLargestGrad(const ModuleBase::matrix& _force,UnitCell& ucell); void GetPos(UnitCell& ucell,std::vector>& pos); diff --git a/source/source_relax/matrix_methods.cpp b/source/source_relax/matrix_methods.cpp index 99019f0668..27f106ae56 100644 --- a/source/source_relax/matrix_methods.cpp +++ b/source/source_relax/matrix_methods.cpp @@ -4,6 +4,8 @@ std::vector ReshapeMToV(std::vector>& matrix) { + assert(!matrix.empty()); + assert(matrix[0].size() == 3); int size = matrix.size(); std::vector result; result.reserve(3*size); @@ -16,6 +18,8 @@ std::vector ReshapeMToV(std::vector>& matrix) std::vector> MAddM(std::vector>& a, std::vector>& b) { + assert(!a.empty() && !b.empty()); + assert(a.size() == b.size() && a[0].size() == b[0].size()); std::vector> result = std::vector>(a.size(), std::vector(a[0].size(), 0.0)); for(int i = 0; i < a.size(); i++) { @@ -29,6 +33,7 @@ std::vector> MAddM(std::vector>& a, std::vector VSubV(std::vector& a, std::vector& b) { + assert(a.size() == b.size()); std::vector result = std::vector(a.size(), 0.0); for(int i = 0; i < a.size(); i++) { @@ -39,6 +44,7 @@ std::vector VSubV(std::vector& a, std::vector& b) std::vector> ReshapeVToM(std::vector& matrix) { + assert(matrix.size() % 3 == 0); std::vector> result = std::vector>(matrix.size() / 3, std::vector(3)); for(int i = 0; i < result.size(); i++) { @@ -52,6 +58,8 @@ std::vector> ReshapeVToM(std::vector& matrix) std::vector DotInMAndV1(std::vector>& matrix, std::vector& vec) { + assert(!matrix.empty()); + assert(matrix[0].size() == vec.size()); std::vector result(matrix.size(), 0.0); for(int i = 0; i < result.size(); i++) { @@ -64,6 +72,8 @@ std::vector DotInMAndV1(std::vector>& matrix, std::v } std::vector DotInMAndV2(std::vector>& matrix, std::vector& vec) { + assert(!matrix.empty()); + assert(matrix.size() == vec.size()); std::vector result(matrix.size(), 0.0); for(int i = 0; i < result.size(); i++) { @@ -77,6 +87,7 @@ std::vector DotInMAndV2(std::vector>& matrix, std::v double DotInVAndV(std::vector& vec1, std::vector& vec2) { + assert(vec1.size() == vec2.size()); double result = 0.0; for(int i = 0; i < vec1.size(); i++) { @@ -87,6 +98,7 @@ double DotInVAndV(std::vector& vec1, std::vector& vec2) std::vector> OuterVAndV(std::vector& a, std::vector& b) { + assert(a.size() == b.size()); std::vector> result = std::vector>(a.size(), std::vector(b.size(), 0.0)); for(int i = 0; i < a.size(); i++) { @@ -100,6 +112,8 @@ std::vector> OuterVAndV(std::vector& a, std::vector< std::vector> MPlus(std::vector>& a, double b) { + assert(!a.empty()); + assert(b != 0); std::vector> result = std::vector>(a.size(), std::vector(a[0].size(), 0.0)); for(int i = 0; i < a.size(); i++) { @@ -113,6 +127,8 @@ std::vector> MPlus(std::vector>& a, doub std::vector> MSubM(std::vector>& a, std::vector>& b) { + assert(!a.empty() && !b.empty()); + assert(a.size() == b.size() && a[0].size() == b[0].size()); std::vector> result = std::vector>(a.size(), std::vector(a[0].size(), 0.0)); for(int i = 0; i < a.size(); i++) { @@ -126,6 +142,7 @@ std::vector> MSubM(std::vector>& a, std: std::vector DotInVAndFloat(std::vector& vec, double b) { + assert(b != 0); std::vector result(vec.size(), 0.0); for(int i = 0; i < vec.size(); i++) { @@ -136,6 +153,7 @@ std::vector DotInVAndFloat(std::vector& vec, double b) std::vector VAddV(std::vector& a, std::vector& b) { + assert(a.size() == b.size()); std::vector result = std::vector(a.size(), 0.0); for(int i = 0; i < a.size(); i++) { diff --git a/source/source_relax/matrix_methods.h b/source/source_relax/matrix_methods.h index 7f576c53af..c87979d52e 100644 --- a/source/source_relax/matrix_methods.h +++ b/source/source_relax/matrix_methods.h @@ -2,6 +2,7 @@ #define MATRIX_METHODS #include +#include diff --git a/source/source_relax/test/bfgs_test.cpp b/source/source_relax/test/bfgs_test.cpp index dcd9bdaa3c..a9cc3efcf8 100644 --- a/source/source_relax/test/bfgs_test.cpp +++ b/source/source_relax/test/bfgs_test.cpp @@ -1,74 +1,215 @@ -#include +// source/source_relax/test/bfgs_test.cpp +#include "gtest/gtest.h" +#include "gmock/gmock.h" #include "for_test.h" + +#define private public #include "source_relax/bfgs.h" -#include "source_cell/unitcell.h" -#include "source_base/matrix.h" -#include "source_relax/ions_move_basic.h" -#include "source_relax/matrix_methods.h" +#undef private + +#define private public +#include "source_io/module_parameter/parameter.h" +#undef private + +#include "source_relax/ions_move_basic.h" // for Ions_Move_Basic static members + +/************************************************ + * unit tests for BFGS (no MockUnitCell) + ***********************************************/ class BFGSTest : public ::testing::Test { protected: BFGS bfgs; - UnitCell ucell; - std::vector> force; - - void SetUp() override { - int size = 10; - bfgs.allocate(size); - - ucell.ntype = 2; - ucell.lat0 = 1.0; - ucell.nat = 10; - ucell.atoms = new Atom[ucell.ntype]; - for (int i = 0; i < ucell.ntype; i++) { - ucell.atoms[i].na = 5; - ucell.atoms[i].tau = std::vector>(5); - ucell.atoms[i].taud = std::vector>(5); - ucell.atoms[i].mbl = std::vector>(5, {1, 1, 1}); - } + void SetUp() override + { + // Initialize variables before each test + } - force = std::vector>(size, std::vector(3, 0.0)); - for (int i = 0; i < force.size(); ++i) { - for (int j = 0; j < 3; ++j) { - force[i][j] = -0.1 * (i + 1); - } - } + void TearDown() override + { + // nothing global to clean here } }; -TEST_F(BFGSTest, PrepareStep) { - bfgs.PrepareStep(force, bfgs.pos, bfgs.H, bfgs.pos0, bfgs.force0, bfgs.steplength, bfgs.dpos, ucell); - EXPECT_EQ(bfgs.steplength.size(), 10); - for (int i = 0; i < 10; ++i) { - EXPECT_GT(bfgs.steplength[i], 0); - } +// Test whether the allocate() function can correctly allocate memory space +TEST_F(BFGSTest, TestAllocate) +{ + int size = 2; + bfgs.allocate(size); + + // Check if allocated arrays are not empty + EXPECT_FALSE(bfgs.H.empty()); + EXPECT_FALSE(bfgs.pos.empty()); + EXPECT_FALSE(bfgs.pos0.empty()); + EXPECT_FALSE(bfgs.pos_taud.empty()); + EXPECT_FALSE(bfgs.pos_taud0.empty()); + EXPECT_FALSE(bfgs.force.empty()); + EXPECT_FALSE(bfgs.force0.empty()); + EXPECT_FALSE(bfgs.steplength.empty()); + EXPECT_FALSE(bfgs.dpos.empty()); + EXPECT_EQ(bfgs.size, size); + EXPECT_EQ(bfgs.alpha,70); + EXPECT_EQ(bfgs.maxstep,PARAM.inp.relax_bfgs_rmax); + EXPECT_TRUE(bfgs.sign); + EXPECT_EQ(bfgs.largest_grad,0.0); } +// Test if a dimension less than or equal to 0 results in an assertion error +TEST_F(BFGSTest, TestAllocateWithZeroDimension) +{ + int size = 0; + ASSERT_DEATH(bfgs.allocate(size), ""); +} -TEST_F(BFGSTest, AllocateTest) { - BFGS bfgs; - int size = 5; +// Test DetermineStep scaling +TEST_F(BFGSTest, DetermineStepScaling) +{ + int size = 2; bfgs.allocate(size); + std::vector steplength = {1.0, 0.1}; + std::vector> dpos = { + {1.0, 1.0, 1.0}, + {0.1, 0.1, 0.1} + }; + double maxstep = 0.5; + bfgs.DetermineStep(steplength, dpos, maxstep); + + // first atom scaled down to maxstep + EXPECT_NEAR(dpos[0][0], 0.5, 1e-12); + EXPECT_NEAR(dpos[0][1], 0.5, 1e-12); + EXPECT_NEAR(dpos[0][2], 0.5, 1e-12); + + // second atom unchanged (small) + EXPECT_NEAR(dpos[1][0], 0.05, 1e-12); + EXPECT_NEAR(dpos[1][1], 0.05, 1e-12); + EXPECT_NEAR(dpos[1][2], 0.05, 1e-12); +} + +// Test GetPos and GetPostaud without creating extra helper class +TEST_F(BFGSTest, GetPosAndPostaud) +{ + // prepare UnitCell with 1 type and 2 atoms + UnitCell ucell; + ucell.ntype = 1; + ucell.nat = 2; + ucell.lat0 = 2.0; + + // allocate atoms array + 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}); + + // 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; + + // allocate mapping arrays + ucell.iat2it = new int[ucell.nat]; + ucell.iat2ia = new int[ucell.nat]; + int k = 0; + for (int it = 0; it < ucell.ntype; ++it) { + for (int ia = 0; ia < ucell.atoms[it].na; ++ia) { + ucell.iat2it[k] = it; + ucell.iat2ia[k] = ia; + ++k; + } + } + + // allocate bfgs arrays and call getters + bfgs.allocate(ucell.nat); + bfgs.GetPos(ucell, bfgs.pos); + bfgs.GetPostaud(ucell, bfgs.pos_taud); + + // pos is tau * BOHR_TO_A * lat0 + EXPECT_DOUBLE_EQ(bfgs.pos[0][0], ucell.atoms[0].tau[0].x * ModuleBase::BOHR_TO_A * ucell.lat0); + EXPECT_DOUBLE_EQ(bfgs.pos_taud[1][2], ucell.atoms[0].taud[1].z); +} + +// Test CalculateLargestGrad (uses ModuleBase::matrix) +TEST_F(BFGSTest, CalculateLargestGrad) +{ + // UnitCell with 1 type and 2 atoms + UnitCell ucell; + ucell.ntype = 1; + ucell.nat = 2; + ucell.lat0 = 2.0; + + ucell.atoms = new Atom[ucell.ntype]; + ucell.atoms[0].na = 2; + ucell.atoms[0].mbl = std::vector>(2, {1, 1, 1}); - EXPECT_EQ(bfgs.steplength.size(), size); - EXPECT_EQ(bfgs.force0.size(), 3*size); - EXPECT_EQ(bfgs.H.size(), 3*size); - for (const auto& row : bfgs.H) { - EXPECT_EQ(row.size(), 3*size); + // mapping arrays + ucell.iat2it = new int[ucell.nat]; + ucell.iat2ia = new int[ucell.nat]; + int k = 0; + for (int it = 0; it < ucell.ntype; ++it) { + for (int ia = 0; ia < ucell.atoms[it].na; ++ia) { + ucell.iat2it[k] = it; + ucell.iat2ia[k] = ia; + ++k; + } } + + // build force matrix: 2 atoms x 3 components + ModuleBase::matrix force(2, 3); + force(0, 0) = -2.0; // this yields grad component = -(-2.0)*lat0 = 4.0 -> divided by lat0 => 2.0 + force(0, 1) = 0.0; + force(0, 2) = 1.0; + force(1, 0) = 3.0; // this yields abs = 6.0 -> divided by lat0 => 3.0 (this should be largest) + force(1, 1) = -1.0; + force(1, 2) = 0.0; + + bfgs.allocate(ucell.nat); + bfgs.CalculateLargestGrad(force, ucell); + + // expected largest_grad = 3.0 (see calculation above) + EXPECT_NEAR(bfgs.largest_grad, 6.0, 1e-12); } -TEST_F(BFGSTest, FullStepTest) -{ - BFGS bfgs; - UnitCell ucell; - ModuleBase::matrix force(3, 3); - int size = 3; - bfgs.allocate(size); - force(0, 0)=-0.5; - force(1, 1)=-0.3; - force(2, 2)=0.1; - EXPECT_EQ(bfgs.force.size(), size); - EXPECT_EQ(bfgs.pos.size(), size); +// Test relax_step basic functionality +TEST_F(BFGSTest, RelaxStepBasic) +{ + // Setup UnitCell with 1 type, 2 atoms + 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 = new int[ucell.nat]; + ucell.iat2ia = new int[ucell.nat]; + int k = 0; + for (int it = 0; it < ucell.ntype; ++it) { + for (int ia = 0; ia < ucell.atoms[it].na; ++ia) { + ucell.iat2it[k] = it; + ucell.iat2ia[k] = ia; + ++k; + } + } + // Set initial positions + 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; + // Setup force matrix + 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; + // Allocate and call relax_step + bfgs.allocate(ucell.nat); + bfgs.relax_step(force, ucell); + // Check that ionic_position_updated is true + EXPECT_TRUE(ucell.ionic_position_updated); + // Check that force values are set (converted units) + EXPECT_NEAR(bfgs.force[0][0], 0.1 * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A, 1e-12); + EXPECT_NEAR(bfgs.force[1][0], -0.1 * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A, 1e-12); + // Check that positions are updated (not equal to initial) + EXPECT_NEAR(bfgs.pos[0][0], 0.0, 1e-12); + EXPECT_NE(bfgs.pos[1][0], 1.0); } \ No newline at end of file