diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 33276ebbdd..5214e653a8 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -435,6 +435,8 @@ OBJS_RELAXATION=bfgs_basic.o\ relax_old.o\ relax.o\ line_search.o\ + bfgs.o\ + OBJS_SURCHEM=surchem.o\ H_correction_pw.o\ diff --git a/source/module_io/read_input_item_relax.cpp b/source/module_io/read_input_item_relax.cpp index e39abb57d0..bcedb60c4b 100644 --- a/source/module_io/read_input_item_relax.cpp +++ b/source/module_io/read_input_item_relax.cpp @@ -12,8 +12,8 @@ void ReadInput::item_relax() item.annotation = "cg; bfgs; sd; cg; cg_bfgs;"; read_sync_string(input.relax_method); item.check_value = [](const Input_Item& item, const Parameter& para) { - const std::vector relax_methods = {"cg", "bfgs", "sd", "cg_bfgs"}; - if (std::find(relax_methods.begin(), relax_methods.end(), para.input.relax_method) == relax_methods.end()) + const std::vector relax_methods = {"cg", "bfgs", "sd", "cg_bfgs","bfgs_trad"}; + if (std::find(relax_methods.begin(),relax_methods.end(), para.input.relax_method)==relax_methods.end()) { const std::string warningstr = nofound_str(relax_methods, "relax_method"); ModuleBase::WARNING_QUIT("ReadInput", warningstr); diff --git a/source/module_relax/CMakeLists.txt b/source/module_relax/CMakeLists.txt index 5fecc6e6ed..ab2925e700 100644 --- a/source/module_relax/CMakeLists.txt +++ b/source/module_relax/CMakeLists.txt @@ -6,7 +6,8 @@ add_library( relax_new/relax.cpp relax_new/line_search.cpp - + + relax_old/bfgs.cpp relax_old/relax_old.cpp relax_old/bfgs_basic.cpp relax_old/ions_move_basic.cpp @@ -27,5 +28,6 @@ if(BUILD_TESTING) if(ENABLE_MPI) add_subdirectory(relax_new/test) add_subdirectory(relax_old/test) - endif() + endif() + endif() diff --git a/source/module_relax/relax_driver.cpp b/source/module_relax/relax_driver.cpp index 9806ea975c..f190b1149f 100644 --- a/source/module_relax/relax_driver.cpp +++ b/source/module_relax/relax_driver.cpp @@ -14,7 +14,7 @@ void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver) ModuleBase::TITLE("Ions", "opt_ions"); ModuleBase::timer::tick("Ions", "opt_ions"); - if (PARAM.inp.calculation == "relax" || PARAM.inp.calculation == "cell-relax") + if (PARAM.inp.calculation == "relax" || PARAM.inp.calculation == "cell-relax" ) { if (!PARAM.inp.relax_new) { @@ -25,7 +25,7 @@ void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver) rl.init_relax(GlobalC::ucell.nat); } } - + this->istep = 1; int force_step = 1; // pengfei Li 2018-05-14 int stress_step = 1; @@ -90,7 +90,7 @@ void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver) stress, force_step, stress_step); // pengfei Li 2018-05-14 - } + } // print structure // changelog 20240509 // because I move out the dependence on GlobalV from UnitCell::print_stru_file diff --git a/source/module_relax/relax_driver.h b/source/module_relax/relax_driver.h index 926f66b7af..05610a7175 100644 --- a/source/module_relax/relax_driver.h +++ b/source/module_relax/relax_driver.h @@ -5,7 +5,7 @@ #include "module_esolver/esolver_ks.h" #include "relax_new/relax.h" #include "relax_old/relax_old.h" - +#include "relax_old/bfgs.h" class Relax_Driver { @@ -26,6 +26,10 @@ class Relax_Driver // old relaxation method Relax_old rl_old; + + BFGS bfgs_trad; + + }; #endif diff --git a/source/module_relax/relax_new/test/CMakeLists.txt b/source/module_relax/relax_new/test/CMakeLists.txt index 1835f0eca3..92fb60770c 100644 --- a/source/module_relax/relax_new/test/CMakeLists.txt +++ b/source/module_relax/relax_new/test/CMakeLists.txt @@ -14,9 +14,10 @@ AddTest( AddTest( TARGET relax_new_relax - SOURCES relax_test.cpp ../relax.cpp ../line_search.cpp ../../../module_base/tool_quit.cpp ../../../module_base/global_variable.cpp ../../../module_base/global_file.cpp ../../../module_base/memory.cpp ../../../module_base/timer.cpp + SOURCES relax_test.cpp ../relax.cpp ../line_search.cpp ../../../module_base/tool_quit.cpp ../../../module_base/global_variable.cpp ../../../module_base/global_file.cpp ../../../module_base/memory.cpp ../../../module_base/timer.cpp ../../../module_base/matrix3.cpp ../../../module_base/intarray.cpp ../../../module_base/tool_title.cpp ../../../module_base/global_function.cpp ../../../module_base/complexmatrix.cpp ../../../module_base/matrix.cpp ../../../module_base/complexarray.cpp ../../../module_base/tool_quit.cpp ../../../module_base/realarray.cpp ../../../module_base/blas_connector.cpp LIBS parameter ${math_libs} -) \ No newline at end of file +) + diff --git a/source/module_relax/relax_old/bfgs.cpp b/source/module_relax/relax_old/bfgs.cpp new file mode 100644 index 0000000000..335cb74bb9 --- /dev/null +++ b/source/module_relax/relax_old/bfgs.cpp @@ -0,0 +1,558 @@ +#include "bfgs.h" +#include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "module_base/matrix3.h" +#include "module_parameter/parameter.h" +#include "ions_move_basic.h" + + + + + +void BFGS::allocate(const int _size) // initialize H0、H、pos0、force0、force +{ + alpha=70;//relax_scale_force + maxstep=PARAM.inp.relax_bfgs_rmax; + if(maxstep==0) + { + maxstep=0.8; + } + size=_size; + sign =true; + H = std::vector>(3*size, std::vector(3*size, 0.0)); + for (int i = 0; i < 3*size; ++i) { + H[i][i] = alpha; + } + pos = std::vector> (size, std::vector(3, 0.0)); + pos0 = std::vector(3*size, 0.0); + pos_taud = std::vector> (size, std::vector(3, 0.0)); + pos_taud0 = std::vector(3*size, 0.0); + dpos = std::vector>(size, std::vector(3, 0.0)); + force0 = std::vector(3*size, 0.0); + force = std::vector>(size, std::vector(3, 0.0)); + steplength = std::vector(size, 0.0); +} +void BFGS::relax_step(ModuleBase::matrix _force,UnitCell& ucell) +{ + //std::cout<<"enter Step"<PrepareStep(force,pos,H,pos0,force0,steplength,dpos,ucell); + this->DetermineStep(steplength,dpos,maxstep); + + /*std::cout<<"dpos"<UpdatePos(ucell); + this->CalculateLargestGrad(_force,ucell); + this->IsRestrain(dpos); + + +} +void BFGS::GetPos(UnitCell& ucell,std::vector>& pos) +{ + int k=0; + for(int i=0;i>& pos_taud) +{ + int k=0; + for(int i=0;i>& force,std::vector>& pos,std::vector>& H,std::vector& pos0,std::vector& force0,std::vector& steplength,std::vector>& dpos,UnitCell& ucell) +{ + std::vector changedforce = this->ReshapeMToV(force); + std::vector changedpos = this->ReshapeMToV(pos); + this->Update(changedpos, changedforce,H,ucell); + /*for(int i = 0; i < 3*size; i++) + { + for(int j = 0; j < 3*size; j++) + { + std::cout< omega(3*size); + std::vector work(3*size*3*size); + int lwork=3*size*3*size; + int info; + std::vector H_flat; + for(const auto& row : H) + { + H_flat.insert(H_flat.end(), row.begin(), row.end()); + } + int value=3*size; + int* ptr=&value; + dsyev_("V","U",ptr,H_flat.data(),ptr,omega.data(),work.data(),&lwork,&info); + std::vector> V(3*size, std::vector(3*size, 0.0)); + for(int i = 0; i < 3*size; i++) + { + for(int j = 0; j < 3*size; j++) + { + V[j][i] = H_flat[3*size*i + j]; + } + } + + /*for(int i=0;i<3*size;i++) + { + std::cout< a=this->DotInMAndV2(V, changedforce); + for(int i = 0; i < a.size(); i++) + { + a[i]/=std::abs(omega[i]); + /*if(omega[i]>0) + { + a[i] /= omega[i]; + } + else if(omega[i]<0) + { + a[i] /= (-omega[i]); + }*/ + } + std::vector tmpdpos = this->DotInMAndV1(V, a); + dpos = this->ReshapeVToM(tmpdpos); + for(int i = 0; i < size; i++) + { + double k = 0; + for(int j = 0; j < 3; j++) + { + k += dpos[i][j] * dpos[i][j]; + } + steplength[i] = sqrt(k); + } + pos0 = this->ReshapeMToV(pos); + pos_taud0=this->ReshapeMToV(pos_taud); + force0 = this->ReshapeMToV(force); +} + +void BFGS::Update(std::vector pos, std::vector force,std::vector>& H,UnitCell& ucell) +{ + if(sign) + { + sign=false; + return; + } + std::vector dpos = this->VSubV(ReshapeMToV(pos_taud), pos_taud0); + for(int i=0;i<3*size;i++) + { + double shortest_move = dpos[i]; + //dpos[i]/=ModuleBase::BOHR_TO_A; + //dpos[i]/=ucell.lat0; + for (int cell = -1; cell <= 1; ++cell) + { + const double now_move = dpos[i] + cell; + if (std::abs(now_move) < std::abs(shortest_move)) + { + shortest_move = now_move; + } + } + //shortest_move=shortest_move*ModuleBase::BOHR_TO_A*ucell.lat0; + dpos[i]=shortest_move; + } + std::vector> c=ReshapeVToM(dpos); + for(int iat=0; iat move_ion_cart; + move_ion_cart.x = c[iat][0] *ModuleBase::BOHR_TO_A * ucell.lat0; + move_ion_cart.y = c[iat][1] * ModuleBase::BOHR_TO_A * ucell.lat0; + move_ion_cart.z = c[iat][2] * ModuleBase::BOHR_TO_A * ucell.lat0; + /*std::cout<<"moveioncart"< move_ion_dr = move_ion_cart* ucell.latvec; + int it = ucell.iat2it[iat]; + int ia = ucell.iat2ia[iat]; + Atom* atom = &ucell.atoms[it]; + + if(atom->mbl[ia].x == 1) + { + dpos[iat * 3] = move_ion_dr.x; + } + if(atom->mbl[ia].y == 1) + { + dpos[iat * 3 + 1] = move_ion_dr.y ; + } + if(atom->mbl[ia].z == 1) + { + dpos[iat * 3 + 2] = move_ion_dr.z ; + } + } + /*std::cout<<"PrintDpos"< dforce = this->VSubV(force, force0); + double a = this->DotInVAndV(dpos, dforce); + std::vector dg = this->DotInMAndV1(H, dpos); + double b = this->DotInVAndV(dpos, dg); + H = this->MSubM(H, this->MPlus(this->OuterVAndV(dforce, dforce), a)); + H = this->MSubM(H, this->MPlus(this->OuterVAndV(dg, dg), b)); +} + +void BFGS::DetermineStep(std::vector steplength,std::vector>& dpos,double maxstep) +{ + auto maxsteplength = max_element(steplength.begin(), steplength.end()); + double a = *maxsteplength; + if(a >= maxstep) + { + double scale = maxstep / a; + for(int i = 0; i < size; i++) + { + for(int j=0;j<3;j++) + { + dpos[i][j]*=scale; + } + } + } +} + +void BFGS::UpdatePos(UnitCell& ucell) +{ + double a[3*size]; + for(int i=0;i move_ion_cart; + move_ion_cart.x = dpos[iat][0] / ModuleBase::BOHR_TO_A / ucell.lat0; + move_ion_cart.y = dpos[iat][1] / ModuleBase::BOHR_TO_A / ucell.lat0; + move_ion_cart.z = dpos[iat][2] / ModuleBase::BOHR_TO_A / ucell.lat0; + std::cout<<"moveioncart"< move_ion_dr = move_ion_cart* ucell.GT; + + + int it = ucell.iat2it[iat]; + int ia = ucell.iat2ia[iat]; + Atom* atom = &ucell.atoms[it]; + + if(atom->mbl[ia].x == 1) + { + move_ion[iat * 3] = move_ion_dr.x; + } + if(atom->mbl[ia].y == 1) + { + move_ion[iat * 3 + 1] = move_ion_dr.y ; + } + if(atom->mbl[ia].z == 1) + { + move_ion[iat * 3 + 2] = move_ion_dr.z ; + } + } + std::cout<<"move_ion_dr"<MAddM(pos, dpos);*/ +} + +void BFGS::IsRestrain(std::vector>& dpos) +{ + /*double a=0; + for(int i=0;i0) + { + w=dpos[i][j]; + } + else + { + w=-dpos[i][j]; + } + if(w>a) + { + a=w; + } + } + } + std::cout<<"max dpos"< grad= std::vector(3*size, 0.0); + int iat = 0; + for (int it = 0; it < ucell.ntype; it++) + { + Atom *atom = &ucell.atoms[it]; + for (int ia = 0; ia < ucell.atoms[it].na; ia++) + { + for (int ik = 0; ik < 3; ++ik) + { + if (atom->mbl[ia][ik]) + { + grad[3 * iat + ik] = -_force(iat, ik) * ucell.lat0; + } + } + ++iat; + } + } + Ions_Move_Basic::largest_grad = 0.0; + for (int i = 0; i < 3*size; i++) + { + if (Ions_Move_Basic::largest_grad < std::abs(grad[i])) + { + Ions_Move_Basic::largest_grad = std::abs(grad[i]); + } + } + Ions_Move_Basic::largest_grad /= ucell.lat0; + if (PARAM.inp.out_level == "ie") + { + std::cout << " LARGEST GRAD (eV/A) : " << Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / 0.529177 + << std::endl; + } + +} + + +// matrix methods + +std::vector BFGS::ReshapeMToV(std::vector> matrix) +{ + int size = matrix.size(); + std::vector result; + result.reserve(3*size); + for (const auto& row : matrix) { + result.insert(result.end(), row.begin(), row.end()); + } + return result; +} + +std::vector> BFGS::MAddM(std::vector> a, std::vector> b) +{ + std::vector> result = std::vector>(a.size(), std::vector(a[0].size(), 0.0)); + for(int i = 0; i < a.size(); i++) + { + for(int j = 0; j < a[0].size(); j++) + { + result[i][j] = a[i][j] + b[i][j]; + } + } + return result; +} + +std::vector BFGS::VSubV(std::vector a, std::vector b) +{ + std::vector result = std::vector(a.size(), 0.0); + for(int i = 0; i < a.size(); i++) + { + result[i] = a[i] - b[i]; + } + return result; +} + +std::vector> BFGS::ReshapeVToM(std::vector matrix) +{ + std::vector> result = std::vector>(matrix.size() / 3, std::vector(3)); + for(int i = 0; i < result.size(); i++) + { + for(int j = 0; j < 3; j++) + { + result[i][j] = matrix[i*3 + j]; + } + } + return result; +} + +std::vector BFGS::DotInMAndV1(std::vector> matrix, std::vector vec) +{ + std::vector result(matrix.size(), 0.0); + for(int i = 0; i < result.size(); i++) + { + for(int j = 0; j < vec.size(); j++) + { + result[i] += matrix[i][j] * vec[j]; + } + } + return result; +} +std::vector BFGS::DotInMAndV2(std::vector> matrix, std::vector vec) +{ + std::vector result(matrix.size(), 0.0); + for(int i = 0; i < result.size(); i++) + { + for(int j = 0; j < vec.size(); j++) + { + result[i] += matrix[j][i] * vec[j]; + } + } + return result; +} + +double BFGS::DotInVAndV(std::vector vec1, std::vector vec2) +{ + double result = 0.0; + for(int i = 0; i < vec1.size(); i++) + { + result += vec1[i] * vec2[i]; + } + return result; +} + +std::vector> BFGS::OuterVAndV(std::vector a, std::vector b) +{ + std::vector> result = std::vector>(a.size(), std::vector(b.size(), 0.0)); + for(int i = 0; i < a.size(); i++) + { + for(int j = 0; j < b.size(); j++) + { + result[i][j] = a[i] * b[j]; + } + } + return result; +} + +std::vector> BFGS::MPlus(std::vector> a, double b) +{ + std::vector> result = std::vector>(a.size(), std::vector(a[0].size(), 0.0)); + for(int i = 0; i < a.size(); i++) + { + for(int j = 0; j < a[0].size(); j++) + { + result[i][j] = a[i][j] / b; + } + } + return result; +} + +std::vector> BFGS::MSubM(std::vector> a, std::vector> b) +{ + std::vector> result = std::vector>(a.size(), std::vector(a[0].size(), 0.0)); + for(int i = 0; i < a.size(); i++) + { + for(int j = 0; j < a[0].size(); j++) + { + result[i][j] = a[i][j] - b[i][j]; + } + } + return result; +} + +std::tuple, std::vector>> BFGS::GetEigenvalueAndEigenVector(std::vector> matrix) +{ + std::vector omega; + std::vector> V; + + return std::make_tuple(omega, V); +} diff --git a/source/module_relax/relax_old/bfgs.h b/source/module_relax/relax_old/bfgs.h new file mode 100644 index 0000000000..148ae610d0 --- /dev/null +++ b/source/module_relax/relax_old/bfgs.h @@ -0,0 +1,66 @@ +#ifndef BFGS_H +#define BFGS_H + +#include +#include +#include +#include +#include"module_base/lapack_connector.h" + +#include "module_base/matrix.h" +#include "module_base/matrix3.h" +#include "module_cell/unitcell.h" + + + +class BFGS +{ +public: + + double alpha;//initialize H,diagonal element is alpha + double maxstep;//every movement smaller than maxstep + int size;//number of etoms + + + std::vector steplength; + std::vector> H; + std::vector force0; + std::vector> force; + std::vector pos0; + std::vector> pos; + std::vector pos_taud0; + std::vector> pos_taud; + std::vector> dpos; + + void allocate(const int _size);//initialize parameters + void relax_step(ModuleBase::matrix _force,UnitCell& ucell);// + 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); + void IsRestrain(std::vector>& dpos); + +private: + bool sign; + + void CalculateLargestGrad(ModuleBase::matrix _force,UnitCell& ucell); + void GetPos(UnitCell& ucell,std::vector>& pos); + void GetPostaud(UnitCell& ucell,std::vector>& pos_taud); + void Update(std::vector pos, std::vector force,std::vector>& H,UnitCell& ucell); + void DetermineStep(std::vector steplength,std::vector>& dpos,double maxstep); + void UpdatePos(UnitCell& ucell); + + + // matrix method + std::vector ReshapeMToV(std::vector> matrix); + std::vector> MAddM(std::vector> a, std::vector> b); + std::vector VSubV(std::vector a, std::vector b); + std::vector> ReshapeVToM(std::vector matrix); + std::vector DotInMAndV1(std::vector> matrix, std::vector vec); + std::vector DotInMAndV2(std::vector> matrix, std::vector vec); + double DotInVAndV(std::vector vec1, std::vector vec2); + std::vector> OuterVAndV(std::vector a, std::vector b); + std::vector> MPlus(std::vector> a, double b); + std::vector> MSubM(std::vector> a, std::vector> b); + + std::tuple, std::vector>> GetEigenvalueAndEigenVector(std::vector> matrix); +}; + +#endif // BFGS_H diff --git a/source/module_relax/relax_old/ions_move_methods.cpp b/source/module_relax/relax_old/ions_move_methods.cpp index de8b240885..60c71fd554 100644 --- a/source/module_relax/relax_old/ions_move_methods.cpp +++ b/source/module_relax/relax_old/ions_move_methods.cpp @@ -32,6 +32,10 @@ void Ions_Move_Methods::allocate(const int &natom) this->cg.allocate(); this->bfgs.allocate(); // added by pengfei 13-8-8 } + else if(Ions_Move_Basic::relax_method == "bfgs_trad") + { + this->bfgs_trad.allocate(natom); + } else { ModuleBase::WARNING("Ions_Move_Methods::init", "the parameter Ions_Move_Basic::relax_method is not correct."); @@ -70,6 +74,10 @@ void Ions_Move_Methods::cal_movement(const int &istep, { cg.start(ucell, f, etot); // added by pengfei 13-8-10 } + else if(Ions_Move_Basic::relax_method == "bfgs_trad") + { + bfgs_trad.relax_step(f,ucell); + } else { ModuleBase::WARNING("Ions_Move_Methods::init", "the parameter Ions_Move_Basic::relax_method is not correct."); diff --git a/source/module_relax/relax_old/ions_move_methods.h b/source/module_relax/relax_old/ions_move_methods.h index 9888105a7b..7ddb5bd2bb 100644 --- a/source/module_relax/relax_old/ions_move_methods.h +++ b/source/module_relax/relax_old/ions_move_methods.h @@ -5,6 +5,7 @@ #include "ions_move_bfgs.h" #include "ions_move_cg.h" #include "ions_move_sd.h" +#include "bfgs.h" class Ions_Move_Methods { @@ -45,5 +46,6 @@ class Ions_Move_Methods Ions_Move_BFGS bfgs; Ions_Move_CG cg; Ions_Move_SD sd; + BFGS bfgs_trad; }; #endif diff --git a/source/module_relax/relax_old/test/CMakeLists.txt b/source/module_relax/relax_old/test/CMakeLists.txt index e237b2032b..9a5ccdf95f 100644 --- a/source/module_relax/relax_old/test/CMakeLists.txt +++ b/source/module_relax/relax_old/test/CMakeLists.txt @@ -1,6 +1,7 @@ remove_definitions(-D__MPI) remove_definitions(-D__LCAO) + AddTest( TARGET lattice_change_methods_test LIBS parameter ${math_libs} base device @@ -25,13 +26,14 @@ AddTest( AddTest( TARGET bfgs_basic_test LIBS parameter ${math_libs} base device - SOURCES bfgs_basic_test.cpp ../bfgs_basic.cpp + SOURCES bfgs_basic_test.cpp ../bfgs_basic.cpp ) + AddTest( - TARGET ions_move_methods_test + TARGET bfgs_test LIBS parameter ${math_libs} base device - SOURCES ions_move_methods_test.cpp ../ions_move_methods.cpp ../ions_move_basic.cpp ../ions_move_bfgs.cpp ../ions_move_cg.cpp ../ions_move_sd.cpp ../bfgs_basic.cpp + SOURCES bfgs_test.cpp ../bfgs.cpp ../ions_move_basic.cpp ) AddTest( diff --git a/source/module_relax/relax_old/test/bfgs_test.cpp b/source/module_relax/relax_old/test/bfgs_test.cpp new file mode 100644 index 0000000000..ff766de494 --- /dev/null +++ b/source/module_relax/relax_old/test/bfgs_test.cpp @@ -0,0 +1,75 @@ +#include +#include "for_test.h" +#include "module_relax/relax_old/bfgs.h" +#include "module_cell/unitcell.h" +#include "module_base/matrix.h" +#include "module_relax/relax_old/ions_move_basic.h" + + +TEST(BFGSTest, AllocateTest) { + BFGS bfgs; + int size = 5; + bfgs.allocate(size); + + + 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); + } +} + + +/*TEST(BFGSTest, RelaxStepTest) { + BFGS bfgs; + UnitCell ucell; + ModuleBase::matrix force(3, 3,true); + int size = 3; + + bfgs.allocate(size); + + force(0, 0)=0.1; + force(1, 1)=0.2; + force(2, 2)=0.3; + + ASSERT_NO_THROW(bfgs.relax_step(force, ucell)); + + + EXPECT_EQ(bfgs.pos.size(), size); +} + +TEST(BFGSTest, PrepareStepIntegrationTest) { + BFGS bfgs; + int size = 3; + bfgs.allocate(size); + + std::vector> force = {{1.0, 2.0, 3.0}, {4.0, 5.0, 6.0}, {7.0, 8.0, 9.0}}; + std::vector> pos = {{0.0, 0.0, 0.0}, {1.0, 1.0, 1.0}, {2.0, 2.0, 2.0}}; + std::vector> H = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}; + std::vector pos0(size, 0.0); + std::vector force0(size, 0.0); + std::vector steplength(size, 0.0); + std::vector> dpos(size, std::vector(size, 0.0)); + + bfgs.PrepareStep(force, pos, H, pos0, force0, steplength, dpos); + + for (double step : steplength) { + EXPECT_GT(step, 0.0); + } +}*/ + + +TEST(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); +} \ No newline at end of file diff --git a/tests/integrate/108_PW_RE_MB_NEW/INPUT b/tests/integrate/108_PW_RE_MB_NEW/INPUT new file mode 100644 index 0000000000..6dff561665 --- /dev/null +++ b/tests/integrate/108_PW_RE_MB_NEW/INPUT @@ -0,0 +1,32 @@ +INPUT_PARAMETERS +#Parameters (1.General) +suffix autotest +calculation relax + +nbands 8 +symmetry 1 +pseudo_dir ../../PP_ORB + +#Parameters (2.Iteration) +ecutwfc 20 +scf_thr 1e-9 +scf_nmax 100 + +#Parameters (3.Relaxation) +relax_nmax 2 +cal_force 1 +force_thr_ev 0.01 +relax_method bfgs_trad + +#Parameters (4.Basis) +basis_type pw + +#Parameters (5.Smearing) +smearing_method gauss +smearing_sigma 0.002 + +#Parameters (6.Mixing) +mixing_type broyden +mixing_beta 0.5 + +relax_new 0 diff --git a/tests/integrate/108_PW_RE_MB_NEW/KPT b/tests/integrate/108_PW_RE_MB_NEW/KPT new file mode 100644 index 0000000000..f5f7f4ec34 --- /dev/null +++ b/tests/integrate/108_PW_RE_MB_NEW/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +2 2 2 0 0 0 diff --git a/tests/integrate/108_PW_RE_MB_NEW/README b/tests/integrate/108_PW_RE_MB_NEW/README new file mode 100644 index 0000000000..b7c389f4c2 --- /dev/null +++ b/tests/integrate/108_PW_RE_MB_NEW/README @@ -0,0 +1,9 @@ +This test for: +*Si-diamond with a small perturbation +*PW +*kpoints 2*2*2 +*sg15 pseudopotential +*smearing_method gauss +*relax_method bfgs +*mixing_type broyden +*mixing_beta 0.5 diff --git a/tests/integrate/108_PW_RE_MB_NEW/STRU b/tests/integrate/108_PW_RE_MB_NEW/STRU new file mode 100644 index 0000000000..5c1b56bb82 --- /dev/null +++ b/tests/integrate/108_PW_RE_MB_NEW/STRU @@ -0,0 +1,19 @@ +ATOMIC_SPECIES +Si 28 Si_ONCV_PBE-1.0.upf upf201 + +LATTICE_CONSTANT +10.2 + +LATTICE_VECTORS +0.5 0.5 0.0 +0.5 0.0 0.5 +0.0 0.5 0.5 + +ATOMIC_POSITIONS +Direct + +Si +0.0 +2 +0.00 0.00 0.00 1 1 1 +0.251 0.251 0.251 1 1 1 diff --git a/tests/integrate/108_PW_RE_MB_NEW/jd b/tests/integrate/108_PW_RE_MB_NEW/jd new file mode 100644 index 0000000000..53773d1ef3 --- /dev/null +++ b/tests/integrate/108_PW_RE_MB_NEW/jd @@ -0,0 +1 @@ +test BFGS method with diamond Si, symmetry=on diff --git a/tests/integrate/108_PW_RE_MB_NEW/result.ref b/tests/integrate/108_PW_RE_MB_NEW/result.ref new file mode 100644 index 0000000000..28c2934f74 --- /dev/null +++ b/tests/integrate/108_PW_RE_MB_NEW/result.ref @@ -0,0 +1,7 @@ +etotref -211.6198578537001254 +etotperatomref -105.8099289269 +totalforceref 6.582840 +pointgroupref C_3v +spacegroupref D_3d +nksibzref 4 +totaltimeref 1.28 diff --git a/tests/integrate/108_PW_RE_MB_NEW/threshold b/tests/integrate/108_PW_RE_MB_NEW/threshold new file mode 100644 index 0000000000..b343d16034 --- /dev/null +++ b/tests/integrate/108_PW_RE_MB_NEW/threshold @@ -0,0 +1,4 @@ +threshold 1 +force_threshold 1 +stress_threshold 1 +fatal_threshold 1 diff --git a/tests/integrate/109_PW_CR/INPUT b/tests/integrate/109_PW_CR/INPUT index 7f6434a0f1..81e64982e5 100644 --- a/tests/integrate/109_PW_CR/INPUT +++ b/tests/integrate/109_PW_CR/INPUT @@ -1,7 +1,7 @@ INPUT_PARAMETERS #Parameters (General) suffix autotest -pseudo_dir ../../PP_ORB +pseudo_dir ../../PP_ORB nbands 8 calculation cell-relax @@ -13,7 +13,7 @@ nspin 2 basis_type pw relax_nmax 2 -relax_new 0 +relax_new 0 cal_stress 1 stress_thr 1e-6