Skip to content

Commit d408f32

Browse files
author
Fei Yang
committed
code optimazation
1 parent faee0f0 commit d408f32

File tree

4 files changed

+59
-20
lines changed

4 files changed

+59
-20
lines changed

‎source/source_relax/bfgs.cpp‎

Lines changed: 28 additions & 10 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,12 @@ 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+
int total_atoms = 0;
83+
for(int i = 0; i < ucell.ntype; i++)
84+
{
85+
total_atoms += ucell.atoms[i].na;
86+
}
87+
assert(pos.size() == total_atoms);
7988
int k=0;
8089
for(int i=0;i<ucell.ntype;i++)
8190
{
@@ -92,6 +101,12 @@ void BFGS::GetPos(UnitCell& ucell,std::vector<std::vector<double>>& pos)
92101
void BFGS::GetPostaud(UnitCell& ucell,
93102
std::vector<std::vector<double>>& pos_taud)
94103
{
104+
int total_atoms = 0;
105+
for(int i = 0; i < ucell.ntype; i++)
106+
{
107+
total_atoms += ucell.atoms[i].na;
108+
}
109+
assert(pos_taud.size() == total_atoms);
95110
int k=0;
96111
for(int i=0;i<ucell.ntype;i++)
97112
{
@@ -112,6 +127,7 @@ void BFGS::PrepareStep(std::vector<std::vector<double>>& force,
112127
std::vector<double>& force0,
113128
std::vector<double>& steplength,
114129
std::vector<std::vector<double>>& dpos,
130+
int& size,
115131
UnitCell& ucell)
116132
{
117133
std::vector<double> changedforce = ReshapeMToV(force);
@@ -144,6 +160,7 @@ void BFGS::PrepareStep(std::vector<std::vector<double>>& force,
144160
std::vector<double> a=DotInMAndV2(V, changedforce);
145161
for(int i = 0; i < a.size(); i++)
146162
{
163+
assert(std::abs(omega[i]) > 1e-8);
147164
a[i]/=std::abs(omega[i]);
148165
}
149166
std::vector<double> tmpdpos = DotInMAndV1(V, a);
@@ -243,6 +260,7 @@ void BFGS::DetermineStep(std::vector<double>& steplength,
243260
{
244261
std::vector<double>::iterator maxsteplength = max_element(steplength.begin(), steplength.end());
245262
double a = *maxsteplength;
263+
assert(a > 1e-10);
246264
if(a >= maxstep)
247265
{
248266
double scale = maxstep / a;
@@ -311,8 +329,8 @@ void BFGS::UpdatePos(UnitCell& ucell)
311329

312330
void BFGS::IsRestrain(std::vector<std::vector<double>>& dpos)
313331
{
314-
Ions_Move_Basic::converged = Ions_Move_Basic::largest_grad
315-
* ModuleBase::Ry_to_eV / 0.529177<PARAM.inp.force_thr_ev;
332+
Ions_Move_Basic::converged = largest_grad
333+
* ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A<PARAM.inp.force_thr_ev;
316334
}
317335

318336
void BFGS::CalculateLargestGrad(const ModuleBase::matrix& _force,UnitCell& ucell)
@@ -334,20 +352,20 @@ void BFGS::CalculateLargestGrad(const ModuleBase::matrix& _force,UnitCell& ucell
334352
++iat;
335353
}
336354
}
337-
Ions_Move_Basic::largest_grad = 0.0;
355+
largest_grad = 0.0;
338356
for (int i = 0; i < 3*size; i++)
339357
{
340-
if (Ions_Move_Basic::largest_grad < std::abs(grad[i]))
358+
if (largest_grad < std::abs(grad[i]))
341359
{
342-
Ions_Move_Basic::largest_grad = std::abs(grad[i]);
360+
largest_grad = std::abs(grad[i]);
343361
}
344362
}
363+
assert(ucell.lat0 != 0);
345364
Ions_Move_Basic::largest_grad /= ucell.lat0;
346365
if (PARAM.inp.out_level == "ie")
347366
{
348-
std::cout << " LARGEST GRAD (eV/Angstrom) : " << Ions_Move_Basic::largest_grad
349-
* ModuleBase::Ry_to_eV / 0.5291772109
367+
std::cout << " LARGEST GRAD (eV/Angstrom) : " << largest_grad
368+
* ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A
350369
<< std::endl;
351370
}
352-
353371
}

‎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)