Skip to content

Commit af4ea04

Browse files
19hello19hello
committed
bfgs_trad
1 parent d8dfacb commit af4ea04

File tree

10 files changed

+109
-121
lines changed

10 files changed

+109
-121
lines changed

‎source/module_relax/CMakeLists.txt‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@ add_library(
66

77
relax_new/relax.cpp
88
relax_new/line_search.cpp
9+
relax_new/bfgs.cpp
910

1011
relax_old/relax_old.cpp
1112
relax_old/bfgs_basic.cpp
@@ -17,7 +18,7 @@ add_library(
1718
relax_old/lattice_change_basic.cpp
1819
relax_old/lattice_change_cg.cpp
1920
relax_old/lattice_change_methods.cpp
20-
bfgs_traditional/bfgs.cpp
21+
2122
)
2223

2324
if(ENABLE_COVERAGE)

‎source/module_relax/relax_driver.cpp‎

Lines changed: 3 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@ void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver)
1818
{
1919
if(PARAM.inp.relax_method == "bfgs_trad")
2020
{
21-
mybfgs.initialize(GlobalC::ucell.nat);
21+
bfgs_trad.init_relax(GlobalC::ucell.nat,GlobalC::ucell);
2222
}
2323
else if (!PARAM.inp.relax_new)
2424
{
@@ -83,26 +83,14 @@ void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver)
8383
if (PARAM.inp.calculation == "relax" || PARAM.inp.calculation == "cell-relax")
8484
{
8585
if(PARAM.inp.relax_method == "bfgs_trad")
86-
{
87-
std::vector<std::vector<double>> _force(force.nr,std::vector<double>(force.nc,0));
88-
for(int i = 0; i < force.nr; i++)
89-
{
90-
91-
for(int j=0;j<force.nc;j++)
92-
{
93-
_force[i][j]=force(i,j)*13.605693009/ModuleBase::BOHR_TO_A;
94-
//std::cout<<_force[i][j]<<' ';
95-
}
96-
//std::cout<<std::endl;
97-
}
98-
stop=mybfgs.Step(_force,GlobalC::ucell);
86+
{
87+
stop=bfgs_trad.relax_step(force,GlobalC::ucell);
9988
}
10089

10190
else if (PARAM.inp.relax_new)
10291
{
10392
stop = rl.relax_step(force, stress, this->etot);
10493
}
105-
10694
else
10795
{
10896
stop = rl_old.relax_step(istep,

‎source/module_relax/relax_driver.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@
55
#include "module_esolver/esolver_ks.h"
66
#include "relax_new/relax.h"
77
#include "relax_old/relax_old.h"
8-
#include "bfgs_traditional/bfgs.h"
8+
#include "relax_new/bfgs.h"
99
class Relax_Driver
1010
{
1111

@@ -27,7 +27,7 @@ class Relax_Driver
2727
// old relaxation method
2828
Relax_old rl_old;
2929

30-
bfgs mybfgs;
30+
BFGS bfgs_trad;
3131
};
3232

3333
#endif

‎source/module_relax/bfgs_traditional/bfgs.cpp‎ renamed to ‎source/module_relax/relax_new/bfgs.cpp‎

Lines changed: 57 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -3,10 +3,10 @@
33
#include "module_base/matrix3.h"
44

55

6-
void bfgs::initialize(int _size) // initialize H0、H、pos0、force0、force
6+
void BFGS::init_relax(const int _size,UnitCell& ucell) // initialize H0、H、pos0、force0、force
77
{
88
alpha=70;//relax_scale_force
9-
maxstep=100; //relax_nmax
9+
maxstep=100;
1010
size=_size;
1111
sign =true;
1212
H = std::vector<std::vector<double>>(3*size, std::vector<double>(3*size, 0.0));
@@ -19,36 +19,45 @@ void bfgs::initialize(int _size) // initialize H0、H、pos0、force0、force
1919
force0 = std::vector<double>(3*size, 0.0);
2020
force = std::vector<std::vector<double>>(size, std::vector<double>(3, 0.0));
2121
steplength = std::vector<double>(size, 0.0);
22-
GetPos();
22+
this->GetPos(ucell,pos);
2323

2424
}
2525

26-
bool bfgs::Step(std::vector<std::vector<double>> _force,UnitCell& ucell)
26+
bool BFGS::relax_step(ModuleBase::matrix _force,UnitCell& ucell)
2727
{
2828
std::cout<<"enter Step"<<std::endl;
29-
GlobalC::ucell.ionic_position_updated = true;
30-
force = _force;
29+
ucell.ionic_position_updated = true;
30+
for(int i = 0; i < _force.nr; i++)
31+
{
32+
33+
for(int j=0;j<_force.nc;j++)
34+
{
35+
force[i][j]=_force(i,j)*ModuleBase::Ry_to_eV/ModuleBase::BOHR_TO_A;
36+
std::cout<<force[i][j]<<' ';
37+
}
38+
std::cout<<std::endl;
39+
}
3140
//std::cout<<"enter Step0"<<std::endl;
3241
//std::cout<<size<<std::endl;
3342
//GetPos(ucell);
34-
PrepareStep();
43+
this->PrepareStep(force,pos,H,pos0,force0,steplength);
3544
//std::cout<<"enter Step1"<<std::endl;
36-
DetermineStep();
37-
//std::cout<<"enter Step2"<<std::endl;
38-
/*for(int i=0;i<size;i++)
45+
this->DetermineStep(steplength,dpos,maxstep);
46+
std::cout<<"enter Step2"<<std::endl;
47+
for(int i=0;i<size;i++)
3948
{
4049
for(int j=0;j<3;j++)
4150
{
4251
std::cout<<dpos[i][j]<<' ';
4352
}
4453
std::cout<<std::endl;
45-
}*/
46-
UpdatePos();
54+
}
55+
this->UpdatePos(ucell);
4756
//std::cout<<"enter Step3"<<std::endl;
48-
return IsRestrain();
57+
return this->IsRestrain(dpos);
4958
}
5059

51-
void bfgs::GetPos()
60+
void BFGS::GetPos(UnitCell& ucell,std::vector<std::vector<double>>& pos)
5261
{
5362
int k=0;
5463
for(int i=0;i<GlobalC::ucell.ntype;i++)
@@ -64,22 +73,22 @@ void bfgs::GetPos()
6473
}
6574

6675

67-
void bfgs::PrepareStep()
76+
void BFGS::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)
6877
{
6978
//std::cout<<"enter prepareStep0"<<std::endl;
70-
std::vector<double> changedforce = ReshapeMToV(force);
71-
std::vector<double> changedpos = ReshapeMToV(pos);
79+
std::vector<double> changedforce = this->ReshapeMToV(force);
80+
std::vector<double> changedpos = this->ReshapeMToV(pos);
7281
//std::cout<<"enter prepareStep1"<<std::endl;
73-
Update(changedpos, changedforce);
82+
this->Update(changedpos, changedforce,H);
7483
//std::cout<<"enter prepareStep2"<<std::endl;
75-
/*for(int i = 0; i < 3*size; i++)
84+
for(int i = 0; i < 3*size; i++)
7685
{
7786
for(int j = 0; j < 3*size; j++)
7887
{
7988
std::cout<<H[i][j]<<' ';
8089
}
8190
std::cout<<std::endl;
82-
}*/
91+
}
8392
//call dysev
8493
//std::cout<<size<<std::endl;
8594
std::vector<double> omega(3*size);
@@ -112,16 +121,16 @@ void bfgs::PrepareStep()
112121
std::cout<<omega[i]<<' ';
113122
}
114123
std::cout<<std::endl;*/
115-
std::vector<double> a=DotInMAndV2(V, changedforce);
124+
std::vector<double> a=this->DotInMAndV2(V, changedforce);
116125
for(int i = 0; i < a.size(); i++)
117126
{
118127
a[i] /= abs(omega[i]);
119128
}
120129

121-
std::vector<double> tmpdpos = DotInMAndV1(V, a);
130+
std::vector<double> tmpdpos = this->DotInMAndV1(V, a);
122131

123132

124-
dpos = ReshapeVToM(tmpdpos);
133+
dpos = this->ReshapeVToM(tmpdpos);
125134
for(int i = 0; i < size; i++)
126135
{
127136
double k = 0;
@@ -131,11 +140,11 @@ void bfgs::PrepareStep()
131140
}
132141
steplength[i] = sqrt(k);
133142
}
134-
pos0 = ReshapeMToV(pos);
135-
force0 = ReshapeMToV(force);
143+
pos0 = this->ReshapeMToV(pos);
144+
force0 = this->ReshapeMToV(force);
136145
}
137146

138-
void bfgs::Update(std::vector<double> pos, std::vector<double> force)
147+
void BFGS::Update(std::vector<double> pos, std::vector<double> force,std::vector<std::vector<double>>& H)
139148
{
140149
if(sign)
141150
{
@@ -154,20 +163,20 @@ void bfgs::Update(std::vector<double> pos, std::vector<double> force)
154163
std::cout<<pos0[i]<<' ';
155164
}
156165
std::cout<<std::endl;*/
157-
std::vector<double> dpos = VSubV(pos, pos0);
166+
std::vector<double> dpos = this->VSubV(pos, pos0);
158167
if(*max_element(dpos.begin(), dpos.end()) < 1e-7)
159168
{
160169
return;
161170
}
162-
std::vector<double> dforce = VSubV(force, force0);
163-
double a = DotInVAndV(dpos, dforce);
164-
std::vector<double> dg = DotInMAndV1(H, dpos);
165-
double b = DotInVAndV(dpos, dg);
166-
H = MSubM(H, MPlus(OuterVAndV(dforce, dforce), a));
167-
H = MSubM(H, MPlus(OuterVAndV(dg, dg), b));
171+
std::vector<double> dforce = this->VSubV(force, force0);
172+
double a = this->DotInVAndV(dpos, dforce);
173+
std::vector<double> dg = this->DotInMAndV1(H, dpos);
174+
double b = this->DotInVAndV(dpos, dg);
175+
H = this->MSubM(H, this->MPlus(this->OuterVAndV(dforce, dforce), a));
176+
H = this->MSubM(H, this->MPlus(this->OuterVAndV(dg, dg), b));
168177
}
169178

170-
void bfgs::DetermineStep()
179+
void BFGS::DetermineStep(std::vector<double> steplength,std::vector<std::vector<double>>& dpos,int maxstep)
171180
{
172181
auto maxsteplength = max_element(steplength.begin(), steplength.end());
173182
double a = *maxsteplength;
@@ -184,7 +193,7 @@ void bfgs::DetermineStep()
184193
}
185194
}
186195

187-
void bfgs::UpdatePos()
196+
void BFGS::UpdatePos(UnitCell& ucell)
188197
{
189198
double move_ion[3*size];
190199
ModuleBase::zeros(move_ion, size*3);
@@ -225,10 +234,10 @@ void bfgs::UpdatePos()
225234
}
226235
}
227236
GlobalC::ucell.update_pos_taud(move_ion);
228-
pos = MAddM(pos, dpos);
237+
pos = this->MAddM(pos, dpos);
229238
}
230239

231-
bool bfgs::IsRestrain()
240+
bool BFGS::IsRestrain(std::vector<std::vector<double>>& dpos)
232241
{
233242
double a=0;
234243
for(int i=0;i<size;i++)
@@ -255,7 +264,7 @@ bool bfgs::IsRestrain()
255264
}
256265
// matrix methods
257266

258-
std::vector<double> bfgs::ReshapeMToV(std::vector<std::vector<double>> matrix)
267+
std::vector<double> BFGS::ReshapeMToV(std::vector<std::vector<double>> matrix)
259268
{
260269
int size = matrix.size();
261270
std::vector<double> result;
@@ -266,7 +275,7 @@ std::vector<double> bfgs::ReshapeMToV(std::vector<std::vector<double>> matrix)
266275
return result;
267276
}
268277

269-
std::vector<std::vector<double>> bfgs::MAddM(std::vector<std::vector<double>> a, std::vector<std::vector<double>> b)
278+
std::vector<std::vector<double>> BFGS::MAddM(std::vector<std::vector<double>> a, std::vector<std::vector<double>> b)
270279
{
271280
std::vector<std::vector<double>> result = std::vector<std::vector<double>>(a.size(), std::vector<double>(a[0].size(), 0.0));
272281
for(int i = 0; i < a.size(); i++)
@@ -279,7 +288,7 @@ std::vector<std::vector<double>> bfgs::MAddM(std::vector<std::vector<double>> a,
279288
return result;
280289
}
281290

282-
std::vector<double> bfgs::VSubV(std::vector<double> a, std::vector<double> b)
291+
std::vector<double> BFGS::VSubV(std::vector<double> a, std::vector<double> b)
283292
{
284293
std::vector<double> result = std::vector<double>(a.size(), 0.0);
285294
for(int i = 0; i < a.size(); i++)
@@ -289,7 +298,7 @@ std::vector<double> bfgs::VSubV(std::vector<double> a, std::vector<double> b)
289298
return result;
290299
}
291300

292-
std::vector<std::vector<double>> bfgs::ReshapeVToM(std::vector<double> matrix)
301+
std::vector<std::vector<double>> BFGS::ReshapeVToM(std::vector<double> matrix)
293302
{
294303
std::vector<std::vector<double>> result = std::vector<std::vector<double>>(matrix.size() / 3, std::vector<double>(3));
295304
for(int i = 0; i < result.size(); i++)
@@ -302,7 +311,7 @@ std::vector<std::vector<double>> bfgs::ReshapeVToM(std::vector<double> matrix)
302311
return result;
303312
}
304313

305-
std::vector<double> bfgs::DotInMAndV1(std::vector<std::vector<double>> matrix, std::vector<double> vec)
314+
std::vector<double> BFGS::DotInMAndV1(std::vector<std::vector<double>> matrix, std::vector<double> vec)
306315
{
307316
std::vector<double> result(matrix.size(), 0.0);
308317
for(int i = 0; i < result.size(); i++)
@@ -314,7 +323,7 @@ std::vector<double> bfgs::DotInMAndV1(std::vector<std::vector<double>> matrix, s
314323
}
315324
return result;
316325
}
317-
std::vector<double> bfgs::DotInMAndV2(std::vector<std::vector<double>> matrix, std::vector<double> vec)
326+
std::vector<double> BFGS::DotInMAndV2(std::vector<std::vector<double>> matrix, std::vector<double> vec)
318327
{
319328
std::vector<double> result(matrix.size(), 0.0);
320329
for(int i = 0; i < result.size(); i++)
@@ -327,7 +336,7 @@ std::vector<double> bfgs::DotInMAndV2(std::vector<std::vector<double>> matrix, s
327336
return result;
328337
}
329338

330-
double bfgs::DotInVAndV(std::vector<double> vec1, std::vector<double> vec2)
339+
double BFGS::DotInVAndV(std::vector<double> vec1, std::vector<double> vec2)
331340
{
332341
double result = 0.0;
333342
for(int i = 0; i < vec1.size(); i++)
@@ -337,7 +346,7 @@ double bfgs::DotInVAndV(std::vector<double> vec1, std::vector<double> vec2)
337346
return result;
338347
}
339348

340-
std::vector<std::vector<double>> bfgs::OuterVAndV(std::vector<double> a, std::vector<double> b)
349+
std::vector<std::vector<double>> BFGS::OuterVAndV(std::vector<double> a, std::vector<double> b)
341350
{
342351
std::vector<std::vector<double>> result = std::vector<std::vector<double>>(a.size(), std::vector<double>(b.size(), 0.0));
343352
for(int i = 0; i < a.size(); i++)
@@ -350,7 +359,7 @@ std::vector<std::vector<double>> bfgs::OuterVAndV(std::vector<double> a, std::ve
350359
return result;
351360
}
352361

353-
std::vector<std::vector<double>> bfgs::MPlus(std::vector<std::vector<double>> a, double b)
362+
std::vector<std::vector<double>> BFGS::MPlus(std::vector<std::vector<double>> a, double b)
354363
{
355364
std::vector<std::vector<double>> result = std::vector<std::vector<double>>(a.size(), std::vector<double>(a[0].size(), 0.0));
356365
for(int i = 0; i < a.size(); i++)
@@ -363,7 +372,7 @@ std::vector<std::vector<double>> bfgs::MPlus(std::vector<std::vector<double>> a,
363372
return result;
364373
}
365374

366-
std::vector<std::vector<double>> bfgs::MSubM(std::vector<std::vector<double>> a, std::vector<std::vector<double>> b)
375+
std::vector<std::vector<double>> BFGS::MSubM(std::vector<std::vector<double>> a, std::vector<std::vector<double>> b)
367376
{
368377
std::vector<std::vector<double>> result = std::vector<std::vector<double>>(a.size(), std::vector<double>(a[0].size(), 0.0));
369378
for(int i = 0; i < a.size(); i++)
@@ -376,7 +385,7 @@ std::vector<std::vector<double>> bfgs::MSubM(std::vector<std::vector<double>> a,
376385
return result;
377386
}
378387

379-
std::tuple<std::vector<double>, std::vector<std::vector<double>>> bfgs::GetEigenvalueAndEigenVector(std::vector<std::vector<double>> matrix)
388+
std::tuple<std::vector<double>, std::vector<std::vector<double>>> BFGS::GetEigenvalueAndEigenVector(std::vector<std::vector<double>> matrix)
380389
{
381390
std::vector<double> omega;
382391
std::vector<std::vector<double>> V;

0 commit comments

Comments
 (0)