Skip to content

Commit c8b7066

Browse files
19hello19hello
committed
Refactor the code and add comments
1 parent ca96d15 commit c8b7066

File tree

10 files changed

+360
-1156
lines changed

10 files changed

+360
-1156
lines changed

‎source/Makefile.Objects‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -444,6 +444,7 @@ OBJS_RELAXATION=bfgs_basic.o\
444444
line_search.o\
445445
bfgs.o\
446446
lbfgs.o\
447+
matrix_methods.o\
447448

448449

449450
OBJS_SURCHEM=surchem.o\

‎source/module_relax/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,7 @@ add_library(
1919
relax_old/lattice_change_basic.cpp
2020
relax_old/lattice_change_cg.cpp
2121
relax_old/lattice_change_methods.cpp
22+
relax_old/matrix_methods.cpp
2223
)
2324

2425
if(ENABLE_COVERAGE)

‎source/module_relax/relax_old/bfgs.cpp‎

Lines changed: 21 additions & 228 deletions
Original file line numberDiff line numberDiff line change
@@ -44,16 +44,6 @@ void BFGS::relax_step(const ModuleBase::matrix& _force,UnitCell& ucell)
4444
force[i][j]=_force(i,j)*ModuleBase::Ry_to_eV/ModuleBase::BOHR_TO_A;
4545
}
4646
}
47-
/*std::cout<<"force"<<std::endl;
48-
for(int i=0;i<size;i++)
49-
{
50-
for(int j=0;j<3;j++)
51-
{
52-
std::cout<<force[i][j]<<' ';
53-
}
54-
std::cout<<std::endl;
55-
}*/
56-
5747
int k=0;
5848
for(int i=0;i<ucell.ntype;i++)
5949
{
@@ -77,43 +67,6 @@ void BFGS::relax_step(const ModuleBase::matrix& _force,UnitCell& ucell)
7767

7868
this->PrepareStep(force,pos,H,pos0,force0,steplength,dpos,ucell);
7969
this->DetermineStep(steplength,dpos,maxstep);
80-
/*std::cout<<"dpos"<<std::endl;
81-
for(int i=0;i<force.size();i++)
82-
{
83-
for(int j=0;j<3;j++)
84-
{
85-
std::cout<<dpos[i][j]<<' ';
86-
}
87-
std::cout<<std::endl;
88-
}
89-
std::cout<<"force"<<std::endl;
90-
for(int i=0;i<size;i++)
91-
{
92-
for(int j=0;j<3;j++)
93-
{
94-
std::cout<<force[i][j]<<' ';
95-
}
96-
std::cout<<std::endl;
97-
}*/
98-
/*std::cout<<"dpos"<<std::endl;
99-
for(int i=0;i<size;i++)
100-
{
101-
for(int j=0;j<3;j++)
102-
{
103-
std::cout<<dpos[i][j]<<' ';
104-
}
105-
std::cout<<std::endl;
106-
}*/
107-
/*std::cout<<"pos"<<std::endl;
108-
for(int i=0;i<size;i++)
109-
{
110-
for(int j=0;j<3;j++)
111-
{
112-
std::cout<<pos[i][j]<<' ';
113-
}
114-
std::cout<<std::endl;
115-
}*/
116-
11770
this->UpdatePos(ucell);
11871
this->CalculateLargestGrad(_force,ucell);
11972
this->IsRestrain(dpos);
@@ -159,8 +112,8 @@ void BFGS::PrepareStep(std::vector<std::vector<double>>& force,
159112
std::vector<std::vector<double>>& dpos,
160113
UnitCell& ucell)
161114
{
162-
std::vector<double> changedforce = this->ReshapeMToV(force);
163-
std::vector<double> changedpos = this->ReshapeMToV(pos);
115+
std::vector<double> changedforce = ReshapeMToV(force);
116+
std::vector<double> changedpos = ReshapeMToV(pos);
164117
this->Update(changedpos, changedforce,H,ucell);
165118

166119
//! call dysev
@@ -186,22 +139,13 @@ void BFGS::PrepareStep(std::vector<std::vector<double>>& force,
186139
V[j][i] = H_flat[3*size*i + j];
187140
}
188141
}
189-
std::vector<double> a=this->DotInMAndV2(V, changedforce);
142+
std::vector<double> a=DotInMAndV2(V, changedforce);
190143
for(int i = 0; i < a.size(); i++)
191144
{
192145
a[i]/=std::abs(omega[i]);
193146
}
194-
std::vector<double> tmpdpos = this->DotInMAndV1(V, a);
195-
dpos = this->ReshapeVToM(tmpdpos);
196-
/*std::cout<<"dpos0"<<std::endl;
197-
for(int i=0;i<size;i++)
198-
{
199-
for(int j=0;j<3;j++)
200-
{
201-
std::cout<<dpos[i][j]<<' ';
202-
}
203-
std::cout<<std::endl;
204-
}*/
147+
std::vector<double> tmpdpos = DotInMAndV1(V, a);
148+
dpos = ReshapeVToM(tmpdpos);
205149
for(int i = 0; i < size; i++)
206150
{
207151
double k = 0;
@@ -211,9 +155,9 @@ void BFGS::PrepareStep(std::vector<std::vector<double>>& force,
211155
}
212156
steplength[i] = sqrt(k);
213157
}
214-
pos0 = this->ReshapeMToV(pos);
215-
pos_taud0=this->ReshapeMToV(pos_taud);
216-
force0 = this->ReshapeMToV(force);
158+
pos0 = ReshapeMToV(pos);
159+
pos_taud0=ReshapeMToV(pos_taud);
160+
force0 = ReshapeMToV(force);
217161
}
218162

219163
void BFGS::Update(std::vector<double>& pos,
@@ -227,8 +171,8 @@ void BFGS::Update(std::vector<double>& pos,
227171
return;
228172
}
229173
//std::vector<double> dpos=this->VSubV(pos,pos0);
230-
auto term=this->ReshapeMToV(pos_taud);
231-
std::vector<double> dpos = this->VSubV(term, pos_taud0);
174+
auto term=ReshapeMToV(pos_taud);
175+
std::vector<double> dpos = VSubV(term, pos_taud0);
232176
for(int i=0;i<3*size;i++)
233177
{
234178
double shortest_move = dpos[i];
@@ -275,44 +219,20 @@ void BFGS::Update(std::vector<double>& pos,
275219
dpos[iat * 3 + 2] = move_ion_dr.z ;
276220
}
277221
}
278-
/*std::cout<<"Printpos"<<std::endl;
279-
for(int i=0;i<3*size;i++)
280-
{
281-
std::cout<<pos[i]<<' ';
282-
}
283-
std::cout<<std::endl;
284-
std::cout<<"Printpos0"<<std::endl;
285-
for(int i=0;i<3*size;i++)
286-
{
287-
std::cout<<pos0[i]<<' ';
288-
}
289-
std::cout<<std::endl;*/
290-
/*std::cout<<"PrintDpos"<<std::endl;
291-
for(int i=0;i<3*size;i++)
292-
{
293-
std::cout<<dpos[i]<<' ';
294-
}
295-
std::cout<<std::endl;*/
296222
if(*max_element(dpos.begin(), dpos.end()) < 1e-7)
297223
{
298224
return;
299-
}
300-
301-
std::vector<double> dforce = this->VSubV(force, force0);
302-
double a = this->DotInVAndV(dpos, dforce);
303-
std::vector<double> dg = this->DotInMAndV1(H, dpos);
304-
double b = this->DotInVAndV(dpos, dg);
305-
306-
/*std::cout<<"a"<<std::endl;
307-
std::cout<<a<<std::endl;
308-
std::cout<<"b"<<std::endl;
309-
std::cout<<b<<std::endl;*/
310-
auto term1=this->OuterVAndV(dforce, dforce);
311-
auto term2=this->OuterVAndV(dg, dg);
312-
auto term3=this->MPlus(term1, a);
313-
auto term4=this->MPlus(term2, b);
314-
H = this->MSubM(H, term3);
315-
H = this->MSubM(H, term4);
225+
}
226+
std::vector<double> dforce = VSubV(force, force0);
227+
double a = DotInVAndV(dpos, dforce);
228+
std::vector<double> dg = DotInMAndV1(H, dpos);
229+
double b = DotInVAndV(dpos, dg);
230+
auto term1=OuterVAndV(dforce, dforce);
231+
auto term2=OuterVAndV(dg, dg);
232+
auto term3=MPlus(term1, a);
233+
auto term4=MPlus(term2, b);
234+
H = MSubM(H, term3);
235+
H = MSubM(H, term4);
316236
}
317237

318238
void BFGS::DetermineStep(std::vector<double>& steplength,
@@ -321,10 +241,6 @@ void BFGS::DetermineStep(std::vector<double>& steplength,
321241
{
322242
auto maxsteplength = max_element(steplength.begin(), steplength.end());
323243
double a = *maxsteplength;
324-
/*std::cout<<"maxstep"<<std::endl;
325-
std::cout<<maxstep<<std::endl;
326-
std::cout<<"maxsteplength"<<std::endl;
327-
std::cout<<a<<std::endl;*/
328244
if(a >= maxstep)
329245
{
330246
double scale = maxstep / a;
@@ -433,126 +349,3 @@ void BFGS::CalculateLargestGrad(const ModuleBase::matrix& _force,UnitCell& ucell
433349
}
434350

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

0 commit comments

Comments
 (0)