Skip to content

Commit 5d79cc5

Browse files
19helloFei Yang
andauthored
fix bfgs bug (#7483)
* fix bfgs bug * fix bfgs bug * fix bfgs bug --------- Co-authored-by: Fei Yang <2501213217@stu.pku.edu.cn>
1 parent 660e4d1 commit 5d79cc5

2 files changed

Lines changed: 43 additions & 2 deletions

File tree

source/source_relax/ions_move_bfgs2.cpp

Lines changed: 13 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -80,6 +80,13 @@ bool Ions_Move_BFGS2::relax_step(const ModuleBase::matrix& _force,UnitCell& ucel
8080
this->UpdatePos(ucell);
8181
this->CalculateLargestGrad(_force,ucell);
8282
bool converged = this->IsRestrain();
83+
if (converged)
84+
{
85+
ofs_running << "\n Largest force is " << largest_grad * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A
86+
<< " eV/Angstrom while threshold is "
87+
<< PARAM.inp.force_thr_ev << " eV/Angstrom" << std::endl;
88+
ofs_running << "\n Ion relaxation is converged!" << std::endl;
89+
}
8390
// print out geometry information during bfgs_trad relax
8491
unitcell::print_tau(ucell.atoms,ucell.Coordinate,ucell.ntype,ucell.lat0,ofs_running);
8592
return converged;
@@ -239,7 +246,12 @@ void Ions_Move_BFGS2::Update(std::vector<double>& pos,
239246
}
240247
}
241248
double threshold=1e-7;
242-
if(*max_element(dpos.begin(), dpos.end()) < threshold)
249+
double max_abs_dpos = 0.0;
250+
for (const double value : dpos)
251+
{
252+
max_abs_dpos = std::max(max_abs_dpos, std::abs(value));
253+
}
254+
if (max_abs_dpos < threshold)
243255
{
244256
return;
245257
}

source/source_relax/test/bfgs_test.cpp

Lines changed: 30 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -118,6 +118,35 @@ TEST_F(BFGSTest, DetermineStepScaling)
118118
EXPECT_NEAR(dpos[1][2], 0.05, 1e-12);
119119
}
120120

121+
TEST_F(BFGSTest, UpdateUsesAbsoluteDisplacementThreshold)
122+
{
123+
UnitCell ucell;
124+
ucell.ntype = 1;
125+
ucell.nat = 1;
126+
ucell.lat0 = 1.0;
127+
ucell.latvec.Identity();
128+
ucell.atoms = new Atom[ucell.ntype];
129+
ucell.atoms[0].na = 1;
130+
ucell.atoms[0].mbl = std::vector<ModuleBase::Vector3<int>>(1, {1, 1, 1});
131+
ucell.iat2it = new int[ucell.nat];
132+
ucell.iat2ia = new int[ucell.nat];
133+
ucell.iat2it[0] = 0;
134+
ucell.iat2ia[0] = 0;
135+
136+
bfgs.allocate(ucell.nat);
137+
bfgs.sign = false;
138+
bfgs.pos_taud[0].x = -1.0e-6;
139+
bfgs.force0 = {0.0, 0.0, 0.0};
140+
141+
std::vector<double> pos = {0.0, 0.0, 0.0};
142+
std::vector<double> force = {1.0, 0.0, 0.0};
143+
const double initial_h00 = bfgs.H[0][0];
144+
145+
bfgs.Update(pos, force, bfgs.H, ucell);
146+
147+
EXPECT_NE(bfgs.H[0][0], initial_h00);
148+
}
149+
121150
// Test GetPos and GetPostaud without creating extra helper class
122151
TEST_F(BFGSTest, GetPosAndPostaud)
123152
{
@@ -245,4 +274,4 @@ TEST_F(BFGSTest, RelaxStepBasic)
245274
// Check that positions are updated (not equal to initial)
246275
EXPECT_NEAR(bfgs.pos[0][0], 0.0, 1e-12);
247276
EXPECT_NE(bfgs.pos[1][0], 1.0);
248-
}
277+
}

0 commit comments

Comments
 (0)