Skip to content

Commit c4a9719

Browse files
19hello19hello
committed
Add determine whether an atom is movable code and modify the value of relax_bfgs_rmax to 0.2
1 parent 12d607e commit c4a9719

File tree

3 files changed

+115
-136
lines changed

3 files changed

+115
-136
lines changed

‎source/module_parameter/input_parameter.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -163,7 +163,7 @@ struct Input_para
163163
double press3 = 0;
164164
double relax_bfgs_w1 = 0.01; ///< wolfe condition 1
165165
double relax_bfgs_w2 = 0.5; ///< wolfe condition 2
166-
double relax_bfgs_rmax = 0.8; ///< trust radius max
166+
double relax_bfgs_rmax = 0.2; ///< trust radius max
167167
double relax_bfgs_rmin = 1e-05; ///< trust radius min
168168
double relax_bfgs_init = 0.5; ///< initial move
169169
std::string fixed_axes = "None"; ///< which axes are fixed

‎source/module_relax/relax_old/bfgs.cpp‎

Lines changed: 93 additions & 130 deletions
Original file line numberDiff line numberDiff line change
@@ -4,18 +4,10 @@
44
#include "module_parameter/parameter.h"
55
#include "ions_move_basic.h"
66

7-
8-
9-
10-
117
void BFGS::allocate(const int _size) // initialize H0、H、pos0、force0、force
128
{
13-
alpha=70;//relax_scale_force
9+
alpha=70;
1410
maxstep=PARAM.inp.relax_bfgs_rmax;
15-
if(maxstep==0)
16-
{
17-
maxstep=0.8;
18-
}
1911
size=_size;
2012
sign =true;
2113
H = std::vector<std::vector<double>>(3*size, std::vector<double>(3*size, 0.0));
@@ -31,41 +23,73 @@ void BFGS::allocate(const int _size) // initialize H0、H、pos0、force0、forc
3123
force = std::vector<std::vector<double>>(size, std::vector<double>(3, 0.0));
3224
steplength = std::vector<double>(size, 0.0);
3325
}
26+
3427
void BFGS::relax_step(ModuleBase::matrix _force,UnitCell& ucell)
3528
{
36-
//std::cout<<"enter Step"<<std::endl;
37-
GetPos(ucell,pos);
29+
GetPos(ucell,pos);
3830
GetPostaud(ucell,pos_taud);
3931
ucell.ionic_position_updated = true;
4032
for(int i = 0; i < _force.nr; i++)
4133
{
42-
4334
for(int j=0;j<_force.nc;j++)
4435
{
4536
force[i][j]=_force(i,j)*ModuleBase::Ry_to_eV/ModuleBase::BOHR_TO_A;
46-
//std::cout<<force[i][j]<<' ';
4737
}
48-
//std::cout<<std::endl;
38+
}
39+
int k=0;
40+
for(int i=0;i<ucell.ntype;i++)
41+
{
42+
for(int j=0;j<ucell.atoms[i].na;j++)
43+
{
44+
if(ucell.atoms[i].mbl[j].x==0)
45+
{
46+
force[k+j][0]=0;
47+
}
48+
if(ucell.atoms[i].mbl[j].y==0)
49+
{
50+
force[k+j][1]=0;
51+
}
52+
if(ucell.atoms[i].mbl[j].z==0)
53+
{
54+
force[k+j][2]=0;
55+
}
56+
}
57+
k+=ucell.atoms[i].na;
4958
}
5059
this->PrepareStep(force,pos,H,pos0,force0,steplength,dpos,ucell);
5160
this->DetermineStep(steplength,dpos,maxstep);
52-
53-
/*std::cout<<"dpos"<<std::endl;
61+
/*std::cout<<"force"<<std::endl;
62+
for(int i=0;i<size;i++)
63+
{
64+
for(int j=0;j<3;j++)
65+
{
66+
std::cout<<force[i][j]<<' ';
67+
}
68+
std::cout<<std::endl;
69+
}
70+
std::cout<<"dpos"<<std::endl;
5471
for(int i=0;i<size;i++)
5572
{
5673
for(int j=0;j<3;j++)
5774
{
5875
std::cout<<dpos[i][j]<<' ';
5976
}
6077
std::cout<<std::endl;
78+
}
79+
std::cout<<"pos"<<std::endl;
80+
for(int i=0;i<size;i++)
81+
{
82+
for(int j=0;j<3;j++)
83+
{
84+
std::cout<<pos[i][j]<<' ';
85+
}
86+
std::cout<<std::endl;
6187
}*/
62-
6388
this->UpdatePos(ucell);
6489
this->CalculateLargestGrad(_force,ucell);
65-
this->IsRestrain(dpos);
66-
67-
90+
this->IsRestrain(dpos);
6891
}
92+
6993
void BFGS::GetPos(UnitCell& ucell,std::vector<std::vector<double>>& pos)
7094
{
7195
int k=0;
@@ -96,26 +120,23 @@ void BFGS::GetPostaud(UnitCell& ucell,std::vector<std::vector<double>>& pos_taud
96120
}
97121
}
98122

99-
100-
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,std::vector<std::vector<double>>& dpos,UnitCell& ucell)
123+
void BFGS::PrepareStep(std::vector<std::vector<double>>& force,
124+
std::vector<std::vector<double>>& pos,
125+
std::vector<std::vector<double>>& H,
126+
std::vector<double>& pos0,
127+
std::vector<double>& force0,
128+
std::vector<double>& steplength,
129+
std::vector<std::vector<double>>& dpos,
130+
UnitCell& ucell)
101131
{
102132
std::vector<double> changedforce = this->ReshapeMToV(force);
103133
std::vector<double> changedpos = this->ReshapeMToV(pos);
104134
this->Update(changedpos, changedforce,H,ucell);
105-
/*for(int i = 0; i < 3*size; i++)
106-
{
107-
for(int j = 0; j < 3*size; j++)
108-
{
109-
std::cout<<H[i][j]<<' ';
110-
}
111-
std::cout<<std::endl;
112-
}*/
113-
114135
//call dysev
115136
std::vector<double> omega(3*size);
116137
std::vector<double> work(3*size*3*size);
117138
int lwork=3*size*3*size;
118-
int info;
139+
int info=0;
119140
std::vector<double> H_flat;
120141
for(const auto& row : H)
121142
{
@@ -132,28 +153,22 @@ void BFGS::PrepareStep(std::vector<std::vector<double>>& force,std::vector<std::
132153
V[j][i] = H_flat[3*size*i + j];
133154
}
134155
}
135-
136-
/*for(int i=0;i<3*size;i++)
137-
{
138-
std::cout<<omega[i]<<' ';
139-
}
140-
std::cout<<std::endl;*/
141-
142156
std::vector<double> a=this->DotInMAndV2(V, changedforce);
143157
for(int i = 0; i < a.size(); i++)
144158
{
145-
a[i]/=std::abs(omega[i]);
146-
/*if(omega[i]>0)
147-
{
148-
a[i] /= omega[i];
149-
}
150-
else if(omega[i]<0)
151-
{
152-
a[i] /= (-omega[i]);
153-
}*/
159+
a[i]/=std::abs(omega[i]);
154160
}
155161
std::vector<double> tmpdpos = this->DotInMAndV1(V, a);
156162
dpos = this->ReshapeVToM(tmpdpos);
163+
/*std::cout<<"dpos0"<<std::endl;
164+
for(int i=0;i<size;i++)
165+
{
166+
for(int j=0;j<3;j++)
167+
{
168+
std::cout<<dpos[i][j]<<' ';
169+
}
170+
std::cout<<std::endl;
171+
}*/
157172
for(int i = 0; i < size; i++)
158173
{
159174
double k = 0;
@@ -168,13 +183,14 @@ void BFGS::PrepareStep(std::vector<std::vector<double>>& force,std::vector<std::
168183
force0 = this->ReshapeMToV(force);
169184
}
170185

171-
void BFGS::Update(std::vector<double> pos, std::vector<double> force,std::vector<std::vector<double>>& H,UnitCell& ucell)
186+
void BFGS::Update(std::vector<double>& pos, std::vector<double>& force,std::vector<std::vector<double>>& H,UnitCell& ucell)
172187
{
173188
if(sign)
174189
{
175190
sign=false;
176191
return;
177192
}
193+
//std::vector<double> dpos=this->VSubV(pos,pos0);
178194
std::vector<double> dpos = this->VSubV(ReshapeMToV(pos_taud), pos_taud0);
179195
for(int i=0;i<3*size;i++)
180196
{
@@ -203,19 +219,12 @@ void BFGS::Update(std::vector<double> pos, std::vector<double> force,std::vector
203219
move_ion_cart.x = c[iat][0] *ModuleBase::BOHR_TO_A * ucell.lat0;
204220
move_ion_cart.y = c[iat][1] * ModuleBase::BOHR_TO_A * ucell.lat0;
205221
move_ion_cart.z = c[iat][2] * ModuleBase::BOHR_TO_A * ucell.lat0;
206-
/*std::cout<<"moveioncart"<<std::endl;
207-
std::cout<<move_ion_cart.x<<' ';
208-
std::cout<<move_ion_cart.y<<' ';
209-
std::cout<<move_ion_cart.z<<' ';
210-
std::cout<<std::endl;*/
211-
212222

213223
//convert pos
214224
ModuleBase::Vector3<double> move_ion_dr = move_ion_cart* ucell.latvec;
215225
int it = ucell.iat2it[iat];
216226
int ia = ucell.iat2ia[iat];
217227
Atom* atom = &ucell.atoms[it];
218-
219228
if(atom->mbl[ia].x == 1)
220229
{
221230
dpos[iat * 3] = move_ion_dr.x;
@@ -229,6 +238,18 @@ void BFGS::Update(std::vector<double> pos, std::vector<double> force,std::vector
229238
dpos[iat * 3 + 2] = move_ion_dr.z ;
230239
}
231240
}
241+
/*std::cout<<"Printpos"<<std::endl;
242+
for(int i=0;i<3*size;i++)
243+
{
244+
std::cout<<pos[i]<<' ';
245+
}
246+
std::cout<<std::endl;
247+
std::cout<<"Printpos0"<<std::endl;
248+
for(int i=0;i<3*size;i++)
249+
{
250+
std::cout<<pos0[i]<<' ';
251+
}
252+
std::cout<<std::endl;*/
232253
/*std::cout<<"PrintDpos"<<std::endl;
233254
for(int i=0;i<3*size;i++)
234255
{
@@ -243,6 +264,10 @@ void BFGS::Update(std::vector<double> pos, std::vector<double> force,std::vector
243264
double a = this->DotInVAndV(dpos, dforce);
244265
std::vector<double> dg = this->DotInMAndV1(H, dpos);
245266
double b = this->DotInVAndV(dpos, dg);
267+
/*std::cout<<"a"<<std::endl;
268+
std::cout<<a<<std::endl;
269+
std::cout<<"b"<<std::endl;
270+
std::cout<<b<<std::endl;*/
246271
H = this->MSubM(H, this->MPlus(this->OuterVAndV(dforce, dforce), a));
247272
H = this->MSubM(H, this->MPlus(this->OuterVAndV(dg, dg), b));
248273
}
@@ -251,6 +276,10 @@ void BFGS::DetermineStep(std::vector<double> steplength,std::vector<std::vector<
251276
{
252277
auto maxsteplength = max_element(steplength.begin(), steplength.end());
253278
double a = *maxsteplength;
279+
/*std::cout<<"maxstep"<<std::endl;
280+
std::cout<<maxstep<<std::endl;
281+
std::cout<<"maxsteplength"<<std::endl;
282+
std::cout<<a<<std::endl;*/
254283
if(a >= maxstep)
255284
{
256285
double scale = maxstep / a;
@@ -275,8 +304,9 @@ void BFGS::UpdatePos(UnitCell& ucell)
275304
a[i*3+j]/=ModuleBase::BOHR_TO_A;
276305
}
277306
}
307+
std::cout<<std::endl;
308+
int k=0;
278309
ucell.update_pos_tau(a);
279-
280310
/*double move_ion[3*size];
281311
ModuleBase::zeros(move_ion, size*3);
282312
@@ -290,11 +320,6 @@ void BFGS::UpdatePos(UnitCell& ucell)
290320
move_ion_cart.x = dpos[iat][0] / ModuleBase::BOHR_TO_A / ucell.lat0;
291321
move_ion_cart.y = dpos[iat][1] / ModuleBase::BOHR_TO_A / ucell.lat0;
292322
move_ion_cart.z = dpos[iat][2] / ModuleBase::BOHR_TO_A / ucell.lat0;
293-
std::cout<<"moveioncart"<<std::endl;
294-
std::cout<<move_ion_cart.x<<' ';
295-
std::cout<<move_ion_cart.y<<' ';
296-
std::cout<<move_ion_cart.z<<' ';
297-
std::cout<<std::endl;
298323
299324
//convert to Direct coordinate
300325
//note here the old GT is used
@@ -320,76 +345,16 @@ void BFGS::UpdatePos(UnitCell& ucell)
320345
move_ion[iat * 3 + 2] = move_ion_dr.z ;
321346
}
322347
}
323-
std::cout<<"move_ion_dr"<<std::endl;
324-
for(int i=0;i<3*size;i++)
325-
{
326-
std::cout<<move_ion[i]<<' ';
327-
}
328-
std::cout<<std::endl;
329-
std::cout<<"printtau"<<std::endl;
330-
for(int i=0;i<ucell.ntype;i++)
331-
{
332-
for(int j=0;j<ucell.atoms[i].na;j++)
333-
{
334-
std::cout<<ucell.atoms[i].tau[j].x<<' ';
335-
std::cout<<ucell.atoms[i].tau[j].y<<' ';
336-
std::cout<<ucell.atoms[i].tau[j].z<<' ';
337-
std::cout<<ucell.atoms[i].taud[j].x<<' ';
338-
std::cout<<ucell.atoms[i].taud[j].y<<' ';
339-
std::cout<<ucell.atoms[i].taud[j].z<<' ';
340-
}
341-
std::cout<<std::endl;
342-
}
343-
//ucell.update_pos_taud(move_ion);
344-
345-
std::cout<<"printtau"<<std::endl;
346-
for(int i=0;i<ucell.ntype;i++)
347-
{
348-
for(int j=0;j<ucell.atoms[i].na;j++)
349-
{
350-
std::cout<<ucell.atoms[i].tau[j].x<<' ';
351-
std::cout<<ucell.atoms[i].tau[j].y<<' ';
352-
std::cout<<ucell.atoms[i].tau[j].z<<' ';
353-
std::cout<<ucell.atoms[i].taud[j].x<<' ';
354-
std::cout<<ucell.atoms[i].taud[j].y<<' ';
355-
std::cout<<ucell.atoms[i].taud[j].z<<' ';
356-
}
357-
std::cout<<std::endl;
358-
}
359-
//pos = this->MAddM(pos, dpos);*/
348+
ucell.update_pos_taud(move_ion);
349+
pos = this->MAddM(pos, dpos);*/
360350
}
361351

362352
void BFGS::IsRestrain(std::vector<std::vector<double>>& dpos)
363353
{
364-
/*double a=0;
365-
for(int i=0;i<size;i++)
366-
{
367-
for(int j=0;j<3;j++)
368-
{
369-
double w;
370-
if(dpos[i][j]>0)
371-
{
372-
w=dpos[i][j];
373-
}
374-
else
375-
{
376-
w=-dpos[i][j];
377-
}
378-
if(w>a)
379-
{
380-
a=w;
381-
}
382-
}
383-
}
384-
std::cout<<"max dpos"<<std::endl;
385-
std::cout<<a<<std::endl;
386-
Ions_Move_Basic::converged = a<0.00001;
387-
std::cout<<Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / 0.529177<<std::endl;
388-
std::cout<<PARAM.inp.force_thr_ev<<std::endl;*/
389354
Ions_Move_Basic::converged = Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / 0.529177<PARAM.inp.force_thr_ev;
390355
}
391356

392-
void BFGS::CalculateLargestGrad(ModuleBase::matrix _force,UnitCell& ucell)
357+
void BFGS::CalculateLargestGrad(ModuleBase::matrix& _force,UnitCell& ucell)
393358
{
394359
std::vector<double> grad= std::vector<double>(3*size, 0.0);
395360
int iat = 0;
@@ -419,16 +384,14 @@ void BFGS::CalculateLargestGrad(ModuleBase::matrix _force,UnitCell& ucell)
419384
Ions_Move_Basic::largest_grad /= ucell.lat0;
420385
if (PARAM.inp.out_level == "ie")
421386
{
422-
std::cout << " LARGEST GRAD (eV/A) : " << Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / 0.529177
387+
std::cout << " LARGEST GRAD (eV/A) : " << Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / 0.5291772109
423388
<< std::endl;
424389
}
425390

426391
}
427-
428-
429392
// matrix methods
430393

431-
std::vector<double> BFGS::ReshapeMToV(std::vector<std::vector<double>> matrix)
394+
std::vector<double> BFGS::ReshapeMToV(std::vector<std::vector<double>>& matrix)
432395
{
433396
int size = matrix.size();
434397
std::vector<double> result;
@@ -439,7 +402,7 @@ std::vector<double> BFGS::ReshapeMToV(std::vector<std::vector<double>> matrix)
439402
return result;
440403
}
441404

442-
std::vector<std::vector<double>> BFGS::MAddM(std::vector<std::vector<double>> a, std::vector<std::vector<double>> b)
405+
std::vector<std::vector<double>> BFGS::MAddM(std::vector<std::vector<double>>& a, std::vector<std::vector<double>>& b)
443406
{
444407
std::vector<std::vector<double>> result = std::vector<std::vector<double>>(a.size(), std::vector<double>(a[0].size(), 0.0));
445408
for(int i = 0; i < a.size(); i++)

0 commit comments

Comments
 (0)