Skip to content

Commit 2b23c18

Browse files
19hello19hello
authored
add LBFGS method (#6005)
* update bfgs method * update bfgs method and modify the parameter force in relax_step to be passed by reference * update bfgs method and modify the parameter force in relax_step to be passed by reference * Introduce the differences between the two BFGS methods * Add & in relax_step input parameters * add lbfgs method * add lbfgs method * add lbfgs * lbfgs0 * addlbfgs * add lbfgs method * add lbfgs method * Refactor the code and add comments * Modify comments * remove auto * remove auto * remove auto * remove auto * remove auto
1 parent eef959d commit 2b23c18

16 files changed

+728
-322
lines changed

‎source/Makefile.Objects

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -443,8 +443,10 @@ OBJS_RELAXATION=bfgs_basic.o\
443443
lattice_change_methods.o\
444444
relax_old.o\
445445
relax.o\
446-
line_search.o\
447446
bfgs.o\
447+
lbfgs.o\
448+
matrix_methods.o\
449+
line_search.o\
448450

449451

450452
OBJS_SURCHEM=surchem.o\

‎source/module_io/read_input_item_relax.cpp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@ void ReadInput::item_relax()
1212
item.annotation = "cg; bfgs; sd; cg; cg_bfgs;";
1313
read_sync_string(input.relax_method);
1414
item.check_value = [](const Input_Item& item, const Parameter& para) {
15-
const std::vector<std::string> relax_methods = {"cg", "bfgs", "sd", "cg_bfgs","bfgs_trad"};
15+
const std::vector<std::string> relax_methods = {"cg", "bfgs", "sd", "cg_bfgs","bfgs_trad","lbfgs"};
1616
if (std::find(relax_methods.begin(),relax_methods.end(), para.input.relax_method)==relax_methods.end())
1717
{
1818
const std::string warningstr = nofound_str(relax_methods, "relax_method");

‎source/module_relax/CMakeLists.txt

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@ add_library(
88
relax_new/line_search.cpp
99

1010
relax_old/bfgs.cpp
11+
relax_old/lbfgs.cpp
1112
relax_old/relax_old.cpp
1213
relax_old/bfgs_basic.cpp
1314
relax_old/ions_move_basic.cpp
@@ -18,6 +19,7 @@ add_library(
1819
relax_old/lattice_change_basic.cpp
1920
relax_old/lattice_change_cg.cpp
2021
relax_old/lattice_change_methods.cpp
22+
relax_old/matrix_methods.cpp
2123
)
2224

2325
if(ENABLE_COVERAGE)

‎source/module_relax/relax_driver.cpp

Lines changed: 3 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -12,10 +12,9 @@
1212
#include "module_cell/print_cell.h"
1313

1414
void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver, UnitCell& ucell)
15-
{
15+
{
1616
ModuleBase::TITLE("Ions", "opt_ions");
1717
ModuleBase::timer::tick("Ions", "opt_ions");
18-
1918
if (PARAM.inp.calculation == "relax" || PARAM.inp.calculation == "cell-relax" )
2019
{
2120
if (!PARAM.inp.relax_new)
@@ -27,7 +26,6 @@ void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver, UnitCell& uce
2726
rl.init_relax(ucell.nat);
2827
}
2928
}
30-
3129
this->istep = 1;
3230
int force_step = 1; // pengfei Li 2018-05-14
3331
int stress_step = 1;
@@ -90,7 +88,8 @@ void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver, UnitCell& uce
9088
force,
9189
stress,
9290
force_step,
93-
stress_step); // pengfei Li 2018-05-14
91+
stress_step,
92+
p_esolver); // pengfei Li 2018-05-14
9493
}
9594
// print structure
9695
// changelog 20240509

‎source/module_relax/relax_old/bfgs.cpp

Lines changed: 22 additions & 216 deletions
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,6 @@ void BFGS::allocate(const int _size)
1212
maxstep=PARAM.inp.relax_bfgs_rmax;
1313
size=_size;
1414
sign =true;
15-
1615
H = std::vector<std::vector<double>>(3*size, std::vector<double>(3*size, 0.0));
1716

1817
for (int i = 0; i < 3*size; ++i)
@@ -32,7 +31,6 @@ void BFGS::allocate(const int _size)
3231

3332

3433
void BFGS::relax_step(const ModuleBase::matrix& _force,UnitCell& ucell)
35-
3634
{
3735
GetPos(ucell,pos);
3836
GetPostaud(ucell,pos_taud);
@@ -44,7 +42,6 @@ void BFGS::relax_step(const ModuleBase::matrix& _force,UnitCell& ucell)
4442
force[i][j]=_force(i,j)*ModuleBase::Ry_to_eV/ModuleBase::BOHR_TO_A;
4543
}
4644
}
47-
4845
int k=0;
4946
for(int i=0;i<ucell.ntype;i++)
5047
{
@@ -68,35 +65,6 @@ void BFGS::relax_step(const ModuleBase::matrix& _force,UnitCell& ucell)
6865

6966
this->PrepareStep(force,pos,H,pos0,force0,steplength,dpos,ucell);
7067
this->DetermineStep(steplength,dpos,maxstep);
71-
72-
/*std::cout<<"force"<<std::endl;
73-
for(int i=0;i<size;i++)
74-
{
75-
for(int j=0;j<3;j++)
76-
{
77-
std::cout<<force[i][j]<<' ';
78-
}
79-
std::cout<<std::endl;
80-
}
81-
std::cout<<"dpos"<<std::endl;
82-
for(int i=0;i<size;i++)
83-
{
84-
for(int j=0;j<3;j++)
85-
{
86-
std::cout<<dpos[i][j]<<' ';
87-
}
88-
std::cout<<std::endl;
89-
}
90-
std::cout<<"pos"<<std::endl;
91-
for(int i=0;i<size;i++)
92-
{
93-
for(int j=0;j<3;j++)
94-
{
95-
std::cout<<pos[i][j]<<' ';
96-
}
97-
std::cout<<std::endl;
98-
}*/
99-
10068
this->UpdatePos(ucell);
10169
this->CalculateLargestGrad(_force,ucell);
10270
this->IsRestrain(dpos);
@@ -142,8 +110,8 @@ void BFGS::PrepareStep(std::vector<std::vector<double>>& force,
142110
std::vector<std::vector<double>>& dpos,
143111
UnitCell& ucell)
144112
{
145-
std::vector<double> changedforce = this->ReshapeMToV(force);
146-
std::vector<double> changedpos = this->ReshapeMToV(pos);
113+
std::vector<double> changedforce = ReshapeMToV(force);
114+
std::vector<double> changedpos = ReshapeMToV(pos);
147115
this->Update(changedpos, changedforce,H,ucell);
148116

149117
//! call dysev
@@ -169,22 +137,13 @@ void BFGS::PrepareStep(std::vector<std::vector<double>>& force,
169137
V[j][i] = H_flat[3*size*i + j];
170138
}
171139
}
172-
std::vector<double> a=this->DotInMAndV2(V, changedforce);
140+
std::vector<double> a=DotInMAndV2(V, changedforce);
173141
for(int i = 0; i < a.size(); i++)
174142
{
175143
a[i]/=std::abs(omega[i]);
176144
}
177-
std::vector<double> tmpdpos = this->DotInMAndV1(V, a);
178-
dpos = this->ReshapeVToM(tmpdpos);
179-
/*std::cout<<"dpos0"<<std::endl;
180-
for(int i=0;i<size;i++)
181-
{
182-
for(int j=0;j<3;j++)
183-
{
184-
std::cout<<dpos[i][j]<<' ';
185-
}
186-
std::cout<<std::endl;
187-
}*/
145+
std::vector<double> tmpdpos = DotInMAndV1(V, a);
146+
dpos = ReshapeVToM(tmpdpos);
188147
for(int i = 0; i < size; i++)
189148
{
190149
double k = 0;
@@ -194,9 +153,9 @@ void BFGS::PrepareStep(std::vector<std::vector<double>>& force,
194153
}
195154
steplength[i] = sqrt(k);
196155
}
197-
pos0 = this->ReshapeMToV(pos);
198-
pos_taud0=this->ReshapeMToV(pos_taud);
199-
force0 = this->ReshapeMToV(force);
156+
pos0 = ReshapeMToV(pos);
157+
pos_taud0=ReshapeMToV(pos_taud);
158+
force0 = ReshapeMToV(force);
200159
}
201160

202161
void BFGS::Update(std::vector<double>& pos,
@@ -210,8 +169,8 @@ void BFGS::Update(std::vector<double>& pos,
210169
return;
211170
}
212171
//std::vector<double> dpos=this->VSubV(pos,pos0);
213-
auto term=this->ReshapeMToV(pos_taud);
214-
std::vector<double> dpos = this->VSubV(term, pos_taud0);
172+
std::vector<double> term=ReshapeMToV(pos_taud);
173+
std::vector<double> dpos = VSubV(term, pos_taud0);
215174
for(int i=0;i<3*size;i++)
216175
{
217176
double shortest_move = dpos[i];
@@ -258,56 +217,28 @@ void BFGS::Update(std::vector<double>& pos,
258217
dpos[iat * 3 + 2] = move_ion_dr.z ;
259218
}
260219
}
261-
/*std::cout<<"Printpos"<<std::endl;
262-
for(int i=0;i<3*size;i++)
263-
{
264-
std::cout<<pos[i]<<' ';
265-
}
266-
std::cout<<std::endl;
267-
std::cout<<"Printpos0"<<std::endl;
268-
for(int i=0;i<3*size;i++)
269-
{
270-
std::cout<<pos0[i]<<' ';
271-
}
272-
std::cout<<std::endl;*/
273-
/*std::cout<<"PrintDpos"<<std::endl;
274-
for(int i=0;i<3*size;i++)
275-
{
276-
std::cout<<dpos[i]<<' ';
277-
}
278-
std::cout<<std::endl;*/
279220
if(*max_element(dpos.begin(), dpos.end()) < 1e-7)
280221
{
281222
return;
282-
}
283-
284-
std::vector<double> dforce = this->VSubV(force, force0);
285-
double a = this->DotInVAndV(dpos, dforce);
286-
std::vector<double> dg = this->DotInMAndV1(H, dpos);
287-
double b = this->DotInVAndV(dpos, dg);
288-
289-
/*std::cout<<"a"<<std::endl;
290-
std::cout<<a<<std::endl;
291-
std::cout<<"b"<<std::endl;
292-
std::cout<<b<<std::endl;*/
293-
auto term1=this->OuterVAndV(dforce, dforce);
294-
auto term2=this->OuterVAndV(dg, dg);
295-
auto term3=this->MPlus(term1, a);
296-
auto term4=this->MPlus(term2, b);
297-
H = this->MSubM(H, term3);
298-
H = this->MSubM(H, term4);
223+
}
224+
std::vector<double> dforce = VSubV(force, force0);
225+
double a = DotInVAndV(dpos, dforce);
226+
std::vector<double> dg = DotInMAndV1(H, dpos);
227+
double b = DotInVAndV(dpos, dg);
228+
std::vector<std::vector<double>> term1=OuterVAndV(dforce, dforce);
229+
std::vector<std::vector<double>> term2=OuterVAndV(dg, dg);
230+
std::vector<std::vector<double>> term3=MPlus(term1, a);
231+
std::vector<std::vector<double>> term4=MPlus(term2, b);
232+
H = MSubM(H, term3);
233+
H = MSubM(H, term4);
299234
}
300235

301236
void BFGS::DetermineStep(std::vector<double>& steplength,
302237
std::vector<std::vector<double>>& dpos,
303238
double& maxstep)
304239
{
305-
auto maxsteplength = max_element(steplength.begin(), steplength.end());
240+
std::vector<double>::iterator maxsteplength = max_element(steplength.begin(), steplength.end());
306241
double a = *maxsteplength;
307-
/*std::cout<<"maxstep"<<std::endl;
308-
std::cout<<maxstep<<std::endl;
309-
std::cout<<"maxsteplength"<<std::endl;
310-
std::cout<<a<<std::endl;*/
311242
if(a >= maxstep)
312243
{
313244
double scale = maxstep / a;
@@ -332,8 +263,6 @@ void BFGS::UpdatePos(UnitCell& ucell)
332263
a[i*3+j]/=ModuleBase::BOHR_TO_A;
333264
}
334265
}
335-
std::cout<<std::endl;
336-
int k=0;
337266
unitcell::update_pos_tau(ucell.lat,a,ucell.ntype,ucell.nat,ucell.atoms);
338267
/*double move_ion[3*size];
339268
ModuleBase::zeros(move_ion, size*3);
@@ -418,126 +347,3 @@ void BFGS::CalculateLargestGrad(const ModuleBase::matrix& _force,UnitCell& ucell
418347
}
419348

420349
}
421-
// matrix methods
422-
423-
std::vector<double> BFGS::ReshapeMToV(std::vector<std::vector<double>>& matrix)
424-
{
425-
int size = matrix.size();
426-
std::vector<double> result;
427-
result.reserve(3*size);
428-
for (const auto& row : matrix) {
429-
result.insert(result.end(), row.begin(), row.end());
430-
}
431-
return result;
432-
}
433-
434-
std::vector<std::vector<double>> BFGS::MAddM(std::vector<std::vector<double>>& a,
435-
std::vector<std::vector<double>>& b)
436-
{
437-
std::vector<std::vector<double>> result = std::vector<std::vector<double>>(a.size(), std::vector<double>(a[0].size(), 0.0));
438-
for(int i = 0; i < a.size(); i++)
439-
{
440-
for(int j = 0; j < a[0].size(); j++)
441-
{
442-
result[i][j] = a[i][j] + b[i][j];
443-
}
444-
}
445-
return result;
446-
}
447-
448-
std::vector<double> BFGS::VSubV(std::vector<double>& a, std::vector<double>& b)
449-
{
450-
std::vector<double> result = std::vector<double>(a.size(), 0.0);
451-
for(int i = 0; i < a.size(); i++)
452-
{
453-
result[i] = a[i] - b[i];
454-
}
455-
return result;
456-
}
457-
458-
std::vector<std::vector<double>> BFGS::ReshapeVToM(std::vector<double>& matrix)
459-
{
460-
std::vector<std::vector<double>> result = std::vector<std::vector<double>>(matrix.size() / 3, std::vector<double>(3));
461-
for(int i = 0; i < result.size(); i++)
462-
{
463-
for(int j = 0; j < 3; j++)
464-
{
465-
result[i][j] = matrix[i*3 + j];
466-
}
467-
}
468-
return result;
469-
}
470-
471-
std::vector<double> BFGS::DotInMAndV1(std::vector<std::vector<double>>& matrix, std::vector<double>& vec)
472-
{
473-
std::vector<double> result(matrix.size(), 0.0);
474-
for(int i = 0; i < result.size(); i++)
475-
{
476-
for(int j = 0; j < vec.size(); j++)
477-
{
478-
result[i] += matrix[i][j] * vec[j];
479-
}
480-
}
481-
return result;
482-
}
483-
std::vector<double> BFGS::DotInMAndV2(std::vector<std::vector<double>>& matrix, std::vector<double>& vec)
484-
{
485-
std::vector<double> result(matrix.size(), 0.0);
486-
for(int i = 0; i < result.size(); i++)
487-
{
488-
for(int j = 0; j < vec.size(); j++)
489-
{
490-
result[i] += matrix[j][i] * vec[j];
491-
}
492-
}
493-
return result;
494-
}
495-
496-
double BFGS::DotInVAndV(std::vector<double>& vec1, std::vector<double>& vec2)
497-
{
498-
double result = 0.0;
499-
for(int i = 0; i < vec1.size(); i++)
500-
{
501-
result += vec1[i] * vec2[i];
502-
}
503-
return result;
504-
}
505-
506-
std::vector<std::vector<double>> BFGS::OuterVAndV(std::vector<double>& a, std::vector<double>& b)
507-
{
508-
std::vector<std::vector<double>> result = std::vector<std::vector<double>>(a.size(), std::vector<double>(b.size(), 0.0));
509-
for(int i = 0; i < a.size(); i++)
510-
{
511-
for(int j = 0; j < b.size(); j++)
512-
{
513-
result[i][j] = a[i] * b[j];
514-
}
515-
}
516-
return result;
517-
}
518-
519-
std::vector<std::vector<double>> BFGS::MPlus(std::vector<std::vector<double>>& a, double& b)
520-
{
521-
std::vector<std::vector<double>> result = std::vector<std::vector<double>>(a.size(), std::vector<double>(a[0].size(), 0.0));
522-
for(int i = 0; i < a.size(); i++)
523-
{
524-
for(int j = 0; j < a[0].size(); j++)
525-
{
526-
result[i][j] = a[i][j] / b;
527-
}
528-
}
529-
return result;
530-
}
531-
532-
std::vector<std::vector<double>> BFGS::MSubM(std::vector<std::vector<double>>& a, std::vector<std::vector<double>>& b)
533-
{
534-
std::vector<std::vector<double>> result = std::vector<std::vector<double>>(a.size(), std::vector<double>(a[0].size(), 0.0));
535-
for(int i = 0; i < a.size(); i++)
536-
{
537-
for(int j = 0; j < a[0].size(); j++)
538-
{
539-
result[i][j] = a[i][j] - b[i][j];
540-
}
541-
}
542-
return result;
543-
}

0 commit comments

Comments
 (0)