Skip to content

Commit c66790b

Browse files
19helloFei Yang
19hello
and
Fei Yang
authored
Modify bfgs test code (deepmodeling#6575)
* change BFGS name and make lattice_change_cg and ions_move_cg shorter * change input parameters * change BFGS name and make lattice_change_cg and ions_move_cg shorter * change input parameters * change input parameters * change input parameters * change input parameters * change input parameters * fix INPUT problem * fix INPUT problem * fix INPUT problem * fix INPUT problem * fix INPUT problem * fix INPUT problem * fix INPUT problem * fix INPUT problem * fix INPUT problem * code optimazation * fix INPUT problem * code optimazation * modify bfgstest code * fixcode --------- Co-authored-by: Fei Yang <2501213217@stu.pku.edu.cn>
1 parent e7ee599 commit c66790b

File tree

5 files changed

+248
-75
lines changed

5 files changed

+248
-75
lines changed

‎source/source_relax/bfgs.cpp‎

Lines changed: 22 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -10,10 +10,12 @@
1010
//! initialize H0、H、pos0、force0、force
1111
void BFGS::allocate(const int _size)
1212
{
13+
assert(_size > 0);
1314
alpha=70;//default value in ase is 70
1415
maxstep=PARAM.inp.relax_bfgs_rmax;
1516
size=_size;
16-
sign =true;
17+
largest_grad=0.0;
18+
sign=true;
1719
H = std::vector<std::vector<double>>(3*size, std::vector<double>(3*size, 0.0));
1820

1921
for (int i = 0; i < 3*size; ++i)
@@ -37,6 +39,7 @@ void BFGS::relax_step(const ModuleBase::matrix& _force,UnitCell& ucell)
3739
GetPos(ucell,pos);
3840
GetPostaud(ucell,pos_taud);
3941
ucell.ionic_position_updated = true;
42+
assert(_force.nr == force.size() && _force.nc == force[0].size());
4043
for(int i = 0; i < _force.nr; i++)
4144
{
4245
for(int j=0;j<_force.nc;j++)
@@ -65,7 +68,7 @@ void BFGS::relax_step(const ModuleBase::matrix& _force,UnitCell& ucell)
6568
k+=ucell.atoms[i].na;
6669
}
6770

68-
this->PrepareStep(force,pos,H,pos0,force0,steplength,dpos,ucell);
71+
this->PrepareStep(force,pos,H,pos0,force0,steplength,dpos,size,ucell);
6972
this->DetermineStep(steplength,dpos,maxstep);
7073
this->UpdatePos(ucell);
7174
this->CalculateLargestGrad(_force,ucell);
@@ -76,6 +79,7 @@ void BFGS::relax_step(const ModuleBase::matrix& _force,UnitCell& ucell)
7679

7780
void BFGS::GetPos(UnitCell& ucell,std::vector<std::vector<double>>& pos)
7881
{
82+
assert(pos.size() == ucell.nat);
7983
int k=0;
8084
for(int i=0;i<ucell.ntype;i++)
8185
{
@@ -92,6 +96,7 @@ void BFGS::GetPos(UnitCell& ucell,std::vector<std::vector<double>>& pos)
9296
void BFGS::GetPostaud(UnitCell& ucell,
9397
std::vector<std::vector<double>>& pos_taud)
9498
{
99+
assert(pos_taud.size() == ucell.nat);
95100
int k=0;
96101
for(int i=0;i<ucell.ntype;i++)
97102
{
@@ -112,6 +117,7 @@ void BFGS::PrepareStep(std::vector<std::vector<double>>& force,
112117
std::vector<double>& force0,
113118
std::vector<double>& steplength,
114119
std::vector<std::vector<double>>& dpos,
120+
int& size,
115121
UnitCell& ucell)
116122
{
117123
std::vector<double> changedforce = ReshapeMToV(force);
@@ -142,8 +148,10 @@ void BFGS::PrepareStep(std::vector<std::vector<double>>& force,
142148
}
143149
}
144150
std::vector<double> a=DotInMAndV2(V, changedforce);
151+
double threshold=1e-8;
145152
for(int i = 0; i < a.size(); i++)
146153
{
154+
assert(std::abs(omega[i]) > threshold);
147155
a[i]/=std::abs(omega[i]);
148156
}
149157
std::vector<double> tmpdpos = DotInMAndV1(V, a);
@@ -221,7 +229,8 @@ void BFGS::Update(std::vector<double>& pos,
221229
dpos[iat * 3 + 2] = move_ion_dr.z ;
222230
}
223231
}
224-
if(*max_element(dpos.begin(), dpos.end()) < 1e-7)
232+
double threshold=1e-7;
233+
if(*max_element(dpos.begin(), dpos.end()) < threshold)
225234
{
226235
return;
227236
}
@@ -243,6 +252,8 @@ void BFGS::DetermineStep(std::vector<double>& steplength,
243252
{
244253
std::vector<double>::iterator maxsteplength = max_element(steplength.begin(), steplength.end());
245254
double a = *maxsteplength;
255+
double threshold=1e-10;
256+
assert(a > threshold);
246257
if(a >= maxstep)
247258
{
248259
double scale = maxstep / a;
@@ -311,8 +322,8 @@ void BFGS::UpdatePos(UnitCell& ucell)
311322

312323
void BFGS::IsRestrain(std::vector<std::vector<double>>& dpos)
313324
{
314-
Ions_Move_Basic::converged = Ions_Move_Basic::largest_grad
315-
* ModuleBase::Ry_to_eV / 0.529177<PARAM.inp.force_thr_ev;
325+
Ions_Move_Basic::converged = largest_grad
326+
* ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A<PARAM.inp.force_thr_ev;
316327
}
317328

318329
void BFGS::CalculateLargestGrad(const ModuleBase::matrix& _force,UnitCell& ucell)
@@ -334,20 +345,20 @@ void BFGS::CalculateLargestGrad(const ModuleBase::matrix& _force,UnitCell& ucell
334345
++iat;
335346
}
336347
}
337-
Ions_Move_Basic::largest_grad = 0.0;
348+
largest_grad = 0.0;
338349
for (int i = 0; i < 3*size; i++)
339350
{
340-
if (Ions_Move_Basic::largest_grad < std::abs(grad[i]))
351+
if (largest_grad < std::abs(grad[i]))
341352
{
342-
Ions_Move_Basic::largest_grad = std::abs(grad[i]);
353+
largest_grad = std::abs(grad[i]);
343354
}
344355
}
356+
assert(ucell.lat0 != 0);
345357
Ions_Move_Basic::largest_grad /= ucell.lat0;
346358
if (PARAM.inp.out_level == "ie")
347359
{
348-
std::cout << " LARGEST GRAD (eV/Angstrom) : " << Ions_Move_Basic::largest_grad
349-
* ModuleBase::Ry_to_eV / 0.5291772109
360+
std::cout << " LARGEST GRAD (eV/Angstrom) : " << largest_grad
361+
* ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A
350362
<< std::endl;
351363
}
352-
353364
}

‎source/source_relax/bfgs.h‎

Lines changed: 12 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,17 @@
1313
class BFGS
1414
{
1515
public:
16+
void allocate(const int _size);//initialize parameters
17+
void relax_step(const ModuleBase::matrix& _force,UnitCell& ucell);//a full iteration step
18+
19+
20+
private:
21+
bool sign;//check if this is the first iteration
22+
double alpha;//initialize H,diagonal element is alpha
23+
double maxstep;//every movement smaller than maxstep
24+
double largest_grad;
25+
int size;//number of atoms
26+
1627
std::vector<double> steplength;//the length of atoms displacement
1728
std::vector<std::vector<double>> H;//Hessian matrix
1829
std::vector<double> force0;//force in previous step
@@ -22,17 +33,8 @@ class BFGS
2233
std::vector<double> pos_taud0;//atom pos in previous step(relative coordinates)
2334
std::vector<std::vector<double>> pos_taud;
2435
std::vector<std::vector<double>> dpos;
25-
26-
void allocate(const int _size);//initialize parameters
27-
void relax_step(const ModuleBase::matrix& _force,UnitCell& ucell);//a full iteration step
28-
void PrepareStep(std::vector<std::vector<double>>& force,std::vector<std::vector<double>>& pos,std::vector<std::vector<double>>& H,std::vector<double>& pos0,std::vector<double>& force0,std::vector<double>& steplength,std::vector<std::vector<double>>& dpos,UnitCell& ucell);//calculate the atomic displacement in one iteration step
29-
30-
private:
31-
bool sign;//check if this is the first iteration
32-
double alpha;//initialize H,diagonal element is alpha
33-
double maxstep;//every movement smaller than maxstep
34-
int size;//number of atoms
3536

37+
void PrepareStep(std::vector<std::vector<double>>& force,std::vector<std::vector<double>>& pos,std::vector<std::vector<double>>& H,std::vector<double>& pos0,std::vector<double>& force0,std::vector<double>& steplength,std::vector<std::vector<double>>& dpos,int& size,UnitCell& ucell);//calculate the atomic displacement in one iteration step
3638
void IsRestrain(std::vector<std::vector<double>>& dpos);//check if converged
3739
void CalculateLargestGrad(const ModuleBase::matrix& _force,UnitCell& ucell);
3840
void GetPos(UnitCell& ucell,std::vector<std::vector<double>>& pos);

‎source/source_relax/matrix_methods.cpp‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,8 @@
44

55
std::vector<double> ReshapeMToV(std::vector<std::vector<double>>& matrix)
66
{
7+
assert(!matrix.empty());
8+
assert(matrix[0].size() == 3);
79
int size = matrix.size();
810
std::vector<double> result;
911
result.reserve(3*size);
@@ -16,6 +18,8 @@ std::vector<double> ReshapeMToV(std::vector<std::vector<double>>& matrix)
1618
std::vector<std::vector<double>> MAddM(std::vector<std::vector<double>>& a,
1719
std::vector<std::vector<double>>& b)
1820
{
21+
assert(!a.empty() && !b.empty());
22+
assert(a.size() == b.size() && a[0].size() == b[0].size());
1923
std::vector<std::vector<double>> result = std::vector<std::vector<double>>(a.size(), std::vector<double>(a[0].size(), 0.0));
2024
for(int i = 0; i < a.size(); i++)
2125
{
@@ -29,6 +33,7 @@ std::vector<std::vector<double>> MAddM(std::vector<std::vector<double>>& a,
2933

3034
std::vector<double> VSubV(std::vector<double>& a, std::vector<double>& b)
3135
{
36+
assert(a.size() == b.size());
3237
std::vector<double> result = std::vector<double>(a.size(), 0.0);
3338
for(int i = 0; i < a.size(); i++)
3439
{
@@ -39,6 +44,7 @@ std::vector<double> VSubV(std::vector<double>& a, std::vector<double>& b)
3944

4045
std::vector<std::vector<double>> ReshapeVToM(std::vector<double>& matrix)
4146
{
47+
assert(matrix.size() % 3 == 0);
4248
std::vector<std::vector<double>> result = std::vector<std::vector<double>>(matrix.size() / 3, std::vector<double>(3));
4349
for(int i = 0; i < result.size(); i++)
4450
{
@@ -52,6 +58,8 @@ std::vector<std::vector<double>> ReshapeVToM(std::vector<double>& matrix)
5258

5359
std::vector<double> DotInMAndV1(std::vector<std::vector<double>>& matrix, std::vector<double>& vec)
5460
{
61+
assert(!matrix.empty());
62+
assert(matrix[0].size() == vec.size());
5563
std::vector<double> result(matrix.size(), 0.0);
5664
for(int i = 0; i < result.size(); i++)
5765
{
@@ -64,6 +72,8 @@ std::vector<double> DotInMAndV1(std::vector<std::vector<double>>& matrix, std::v
6472
}
6573
std::vector<double> DotInMAndV2(std::vector<std::vector<double>>& matrix, std::vector<double>& vec)
6674
{
75+
assert(!matrix.empty());
76+
assert(matrix.size() == vec.size());
6777
std::vector<double> result(matrix.size(), 0.0);
6878
for(int i = 0; i < result.size(); i++)
6979
{
@@ -77,6 +87,7 @@ std::vector<double> DotInMAndV2(std::vector<std::vector<double>>& matrix, std::v
7787

7888
double DotInVAndV(std::vector<double>& vec1, std::vector<double>& vec2)
7989
{
90+
assert(vec1.size() == vec2.size());
8091
double result = 0.0;
8192
for(int i = 0; i < vec1.size(); i++)
8293
{
@@ -87,6 +98,7 @@ double DotInVAndV(std::vector<double>& vec1, std::vector<double>& vec2)
8798

8899
std::vector<std::vector<double>> OuterVAndV(std::vector<double>& a, std::vector<double>& b)
89100
{
101+
assert(a.size() == b.size());
90102
std::vector<std::vector<double>> result = std::vector<std::vector<double>>(a.size(), std::vector<double>(b.size(), 0.0));
91103
for(int i = 0; i < a.size(); i++)
92104
{
@@ -100,6 +112,8 @@ std::vector<std::vector<double>> OuterVAndV(std::vector<double>& a, std::vector<
100112

101113
std::vector<std::vector<double>> MPlus(std::vector<std::vector<double>>& a, double b)
102114
{
115+
assert(!a.empty());
116+
assert(b != 0);
103117
std::vector<std::vector<double>> result = std::vector<std::vector<double>>(a.size(), std::vector<double>(a[0].size(), 0.0));
104118
for(int i = 0; i < a.size(); i++)
105119
{
@@ -113,6 +127,8 @@ std::vector<std::vector<double>> MPlus(std::vector<std::vector<double>>& a, doub
113127

114128
std::vector<std::vector<double>> MSubM(std::vector<std::vector<double>>& a, std::vector<std::vector<double>>& b)
115129
{
130+
assert(!a.empty() && !b.empty());
131+
assert(a.size() == b.size() && a[0].size() == b[0].size());
116132
std::vector<std::vector<double>> result = std::vector<std::vector<double>>(a.size(), std::vector<double>(a[0].size(), 0.0));
117133
for(int i = 0; i < a.size(); i++)
118134
{
@@ -126,6 +142,7 @@ std::vector<std::vector<double>> MSubM(std::vector<std::vector<double>>& a, std:
126142

127143
std::vector<double> DotInVAndFloat(std::vector<double>& vec, double b)
128144
{
145+
assert(b != 0);
129146
std::vector<double> result(vec.size(), 0.0);
130147
for(int i = 0; i < vec.size(); i++)
131148
{
@@ -136,6 +153,7 @@ std::vector<double> DotInVAndFloat(std::vector<double>& vec, double b)
136153

137154
std::vector<double> VAddV(std::vector<double>& a, std::vector<double>& b)
138155
{
156+
assert(a.size() == b.size());
139157
std::vector<double> result = std::vector<double>(a.size(), 0.0);
140158
for(int i = 0; i < a.size(); i++)
141159
{

‎source/source_relax/matrix_methods.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@
22
#define MATRIX_METHODS
33

44
#include <vector>
5+
#include <cassert>
56

67

78

0 commit comments

Comments
 (0)