Skip to content

Commit 9c3c405

Browse files
19hello19hello
committed
bfgs_trad
1 parent 8e9f2d8 commit 9c3c405

File tree

12 files changed

+352
-94
lines changed

12 files changed

+352
-94
lines changed

‎source/module_relax/relax_old/bfgs.cpp‎

Lines changed: 186 additions & 51 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,13 @@
11
#include "bfgs.h"
22
#include "module_hamilt_pw/hamilt_pwdft/global.h"
33
#include "module_base/matrix3.h"
4-
#include "module_relax/relax_old/ions_move_basic.h"
54
#include "module_parameter/parameter.h"
65
#include "ions_move_basic.h"
76

87

8+
9+
10+
911
void BFGS::allocate(const int _size) // initialize H0、H、pos0、force0、force
1012
{
1113
alpha=70;//relax_scale_force
@@ -22,18 +24,18 @@ void BFGS::allocate(const int _size) // initialize H0、H、pos0、force0、forc
2224
}
2325
pos = std::vector<std::vector<double>> (size, std::vector<double>(3, 0.0));
2426
pos0 = std::vector<double>(3*size, 0.0);
27+
pos_taud = std::vector<std::vector<double>> (size, std::vector<double>(3, 0.0));
28+
pos_taud0 = std::vector<double>(3*size, 0.0);
2529
dpos = std::vector<std::vector<double>>(size, std::vector<double>(3, 0.0));
2630
force0 = std::vector<double>(3*size, 0.0);
2731
force = std::vector<std::vector<double>>(size, std::vector<double>(3, 0.0));
28-
steplength = std::vector<double>(size, 0.0);
29-
this->GetPos(GlobalC::ucell,pos);
30-
31-
32+
steplength = std::vector<double>(size, 0.0);
3233
}
33-
3434
void BFGS::relax_step(ModuleBase::matrix _force,UnitCell& ucell)
3535
{
3636
//std::cout<<"enter Step"<<std::endl;
37+
GetPos(ucell,pos);
38+
GetPostaud(ucell,pos_taud);
3739
ucell.ionic_position_updated = true;
3840
for(int i = 0; i < _force.nr; i++)
3941
{
@@ -45,27 +47,24 @@ void BFGS::relax_step(ModuleBase::matrix _force,UnitCell& ucell)
4547
}
4648
//std::cout<<std::endl;
4749
}
48-
//std::cout<<"enter Step0"<<std::endl;
49-
//std::cout<<size<<std::endl;
50-
//GetPos(ucell);
51-
this->PrepareStep(force,pos,H,pos0,force0,steplength,dpos);
52-
//std::cout<<"enter Step1"<<std::endl;
50+
this->PrepareStep(force,pos,H,pos0,force0,steplength,dpos,ucell);
5351
this->DetermineStep(steplength,dpos,maxstep);
5452

55-
//std::cout<<"enter Step2"<<std::endl;
56-
/*for(int i=0;i<size;i++)
53+
/*std::cout<<"dpos"<<std::endl;
54+
for(int i=0;i<size;i++)
5755
{
5856
for(int j=0;j<3;j++)
5957
{
6058
std::cout<<dpos[i][j]<<' ';
6159
}
6260
std::cout<<std::endl;
6361
}*/
62+
6463
this->UpdatePos(ucell);
65-
//std::cout<<"enter Step3"<<std::endl;
6664
this->IsRestrain(dpos);
65+
this->CalculateLargestGrad(_force,ucell);
66+
6767
}
68-
6968
void BFGS::GetPos(UnitCell& ucell,std::vector<std::vector<double>>& pos)
7069
{
7170
int k=0;
@@ -81,15 +80,27 @@ void BFGS::GetPos(UnitCell& ucell,std::vector<std::vector<double>>& pos)
8180
}
8281
}
8382

83+
void BFGS::GetPostaud(UnitCell& ucell,std::vector<std::vector<double>>& pos_taud)
84+
{
85+
int k=0;
86+
for(int i=0;i<ucell.ntype;i++)
87+
{
88+
for(int j=0;j<ucell.atoms[i].na;j++)
89+
{
90+
pos_taud[k+j][0]=ucell.atoms[i].taud[j].x;
91+
pos_taud[k+j][1]=ucell.atoms[i].taud[j].y;
92+
pos_taud[k+j][2]=ucell.atoms[i].taud[j].z;
93+
}
94+
k+=ucell.atoms[i].na;
95+
}
96+
}
97+
8498

85-
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)
99+
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)
86100
{
87-
//std::cout<<"enter prepareStep0"<<std::endl;
88101
std::vector<double> changedforce = this->ReshapeMToV(force);
89102
std::vector<double> changedpos = this->ReshapeMToV(pos);
90-
//std::cout<<"enter prepareStep1"<<std::endl;
91-
this->Update(changedpos, changedforce,H);
92-
//std::cout<<"enter prepareStep2"<<std::endl;
103+
this->Update(changedpos, changedforce,H,ucell);
93104
/*for(int i = 0; i < 3*size; i++)
94105
{
95106
for(int j = 0; j < 3*size; j++)
@@ -98,12 +109,10 @@ void BFGS::PrepareStep(std::vector<std::vector<double>>& force,std::vector<std::
98109
}
99110
std::cout<<std::endl;
100111
}*/
112+
101113
//call dysev
102-
//std::cout<<size<<std::endl;
103114
std::vector<double> omega(3*size);
104-
//std::cout<<"enter prepareStep3"<<std::endl;
105115
std::vector<double> work(3*size*3*size);
106-
//std::cout<<"enter prepareStep4"<<std::endl;
107116
int lwork=3*size*3*size;
108117
int info;
109118
std::vector<double> H_flat;
@@ -115,7 +124,6 @@ void BFGS::PrepareStep(std::vector<std::vector<double>>& force,std::vector<std::
115124
int* ptr=&value;
116125
dsyev_("V","U",ptr,H_flat.data(),ptr,omega.data(),work.data(),&lwork,&info);
117126
std::vector<std::vector<double>> V(3*size, std::vector<double>(3*size, 0.0));
118-
//std::cout<<"enter prepareStep5"<<std::endl;
119127
for(int i = 0; i < 3*size; i++)
120128
{
121129
for(int j = 0; j < 3*size; j++)
@@ -124,38 +132,27 @@ void BFGS::PrepareStep(std::vector<std::vector<double>>& force,std::vector<std::
124132
}
125133
}
126134

127-
128135
/*for(int i=0;i<3*size;i++)
129136
{
130137
std::cout<<omega[i]<<' ';
131138
}
132139
std::cout<<std::endl;*/
140+
133141
std::vector<double> a=this->DotInMAndV2(V, changedforce);
134142
for(int i = 0; i < a.size(); i++)
135143
{
136-
if(omega[i]>0)
144+
a[i]/=std::abs(omega[i]);
145+
/*if(omega[i]>0)
137146
{
138147
a[i] /= omega[i];
139148
}
140149
else if(omega[i]<0)
141150
{
142151
a[i] /= (-omega[i]);
143-
}
144-
152+
}*/
145153
}
146-
147154
std::vector<double> tmpdpos = this->DotInMAndV1(V, a);
148-
149-
150155
dpos = this->ReshapeVToM(tmpdpos);
151-
/*for(int i=0;i<size;i++)
152-
{
153-
for(int j=0;j<3;j++)
154-
{
155-
std::cout<<dpos[i][j]<<' ';
156-
}
157-
std::cout<<std::endl;
158-
}*/
159156
for(int i = 0; i < size; i++)
160157
{
161158
double k = 0;
@@ -166,29 +163,77 @@ void BFGS::PrepareStep(std::vector<std::vector<double>>& force,std::vector<std::
166163
steplength[i] = sqrt(k);
167164
}
168165
pos0 = this->ReshapeMToV(pos);
166+
pos_taud0=this->ReshapeMToV(pos_taud);
169167
force0 = this->ReshapeMToV(force);
170168
}
171169

172-
void BFGS::Update(std::vector<double> pos, std::vector<double> force,std::vector<std::vector<double>>& H)
170+
void BFGS::Update(std::vector<double> pos, std::vector<double> force,std::vector<std::vector<double>>& H,UnitCell& ucell)
173171
{
174172
if(sign)
175173
{
176174
sign=false;
177175
return;
178176
}
179-
/*std::cout<<"PrintPos"<<std::endl;
177+
std::vector<double> dpos = this->VSubV(ReshapeMToV(pos_taud), pos_taud0);
180178
for(int i=0;i<3*size;i++)
181179
{
182-
std::cout<<pos[i]<<' ';
180+
double shortest_move = dpos[i];
181+
//dpos[i]/=ModuleBase::BOHR_TO_A;
182+
//dpos[i]/=ucell.lat0;
183+
for (int cell = -1; cell <= 1; ++cell)
184+
{
185+
const double now_move = dpos[i] + cell;
186+
if (std::abs(now_move) < std::abs(shortest_move))
187+
{
188+
shortest_move = now_move;
189+
}
190+
}
191+
//shortest_move=shortest_move*ModuleBase::BOHR_TO_A*ucell.lat0;
192+
dpos[i]=shortest_move;
183193
}
184-
std::cout<<std::endl;
185-
std::cout<<"PrintPos0"<<std::endl;
194+
std::vector<std::vector<double>> c=ReshapeVToM(dpos);
195+
for(int iat=0; iat<size; iat++)
196+
{
197+
//Cartesian coordinate
198+
//convert from Angstrom to unit of latvec (Bohr)
199+
200+
//convert unit
201+
ModuleBase::Vector3<double> move_ion_cart;
202+
move_ion_cart.x = c[iat][0] *ModuleBase::BOHR_TO_A * ucell.lat0;
203+
move_ion_cart.y = c[iat][1] * ModuleBase::BOHR_TO_A * ucell.lat0;
204+
move_ion_cart.z = c[iat][2] * ModuleBase::BOHR_TO_A * ucell.lat0;
205+
/*std::cout<<"moveioncart"<<std::endl;
206+
std::cout<<move_ion_cart.x<<' ';
207+
std::cout<<move_ion_cart.y<<' ';
208+
std::cout<<move_ion_cart.z<<' ';
209+
std::cout<<std::endl;*/
210+
211+
212+
//convert pos
213+
ModuleBase::Vector3<double> move_ion_dr = move_ion_cart* ucell.latvec;
214+
int it = ucell.iat2it[iat];
215+
int ia = ucell.iat2ia[iat];
216+
Atom* atom = &ucell.atoms[it];
217+
218+
if(atom->mbl[ia].x == 1)
219+
{
220+
dpos[iat * 3] = move_ion_dr.x;
221+
}
222+
if(atom->mbl[ia].y == 1)
223+
{
224+
dpos[iat * 3 + 1] = move_ion_dr.y ;
225+
}
226+
if(atom->mbl[ia].z == 1)
227+
{
228+
dpos[iat * 3 + 2] = move_ion_dr.z ;
229+
}
230+
}
231+
/*std::cout<<"PrintDpos"<<std::endl;
186232
for(int i=0;i<3*size;i++)
187233
{
188-
std::cout<<pos0[i]<<' ';
234+
std::cout<<dpos[i]<<' ';
189235
}
190236
std::cout<<std::endl;*/
191-
std::vector<double> dpos = this->VSubV(pos, pos0);
192237
if(*max_element(dpos.begin(), dpos.end()) < 1e-7)
193238
{
194239
return;
@@ -220,7 +265,18 @@ void BFGS::DetermineStep(std::vector<double> steplength,std::vector<std::vector<
220265

221266
void BFGS::UpdatePos(UnitCell& ucell)
222267
{
223-
double move_ion[3*size];
268+
double a[3*size];
269+
for(int i=0;i<size;i++)
270+
{
271+
for(int j=0;j<3;j++)
272+
{
273+
a[i*3+j]=pos[i][j]+dpos[i][j];
274+
a[i*3+j]/=ModuleBase::BOHR_TO_A;
275+
}
276+
}
277+
ucell.update_pos_tau(a);
278+
279+
/*double move_ion[3*size];
224280
ModuleBase::zeros(move_ion, size*3);
225281
226282
for(int iat=0; iat<size; iat++)
@@ -233,12 +289,17 @@ void BFGS::UpdatePos(UnitCell& ucell)
233289
move_ion_cart.x = dpos[iat][0] / ModuleBase::BOHR_TO_A / ucell.lat0;
234290
move_ion_cart.y = dpos[iat][1] / ModuleBase::BOHR_TO_A / ucell.lat0;
235291
move_ion_cart.z = dpos[iat][2] / ModuleBase::BOHR_TO_A / ucell.lat0;
292+
std::cout<<"moveioncart"<<std::endl;
293+
std::cout<<move_ion_cart.x<<' ';
294+
std::cout<<move_ion_cart.y<<' ';
295+
std::cout<<move_ion_cart.z<<' ';
296+
std::cout<<std::endl;
236297
237298
//convert to Direct coordinate
238299
//note here the old GT is used
239300
240301
//convert pos
241-
ModuleBase::Vector3<double> move_ion_dr = move_ion_cart * ucell.GT;
302+
ModuleBase::Vector3<double> move_ion_dr = move_ion_cart* ucell.GT;
242303
243304
244305
int it = ucell.iat2it[iat];
@@ -258,8 +319,43 @@ void BFGS::UpdatePos(UnitCell& ucell)
258319
move_ion[iat * 3 + 2] = move_ion_dr.z ;
259320
}
260321
}
261-
ucell.update_pos_taud(move_ion);
262-
pos = this->MAddM(pos, dpos);
322+
std::cout<<"move_ion_dr"<<std::endl;
323+
for(int i=0;i<3*size;i++)
324+
{
325+
std::cout<<move_ion[i]<<' ';
326+
}
327+
std::cout<<std::endl;
328+
std::cout<<"printtau"<<std::endl;
329+
for(int i=0;i<ucell.ntype;i++)
330+
{
331+
for(int j=0;j<ucell.atoms[i].na;j++)
332+
{
333+
std::cout<<ucell.atoms[i].tau[j].x<<' ';
334+
std::cout<<ucell.atoms[i].tau[j].y<<' ';
335+
std::cout<<ucell.atoms[i].tau[j].z<<' ';
336+
std::cout<<ucell.atoms[i].taud[j].x<<' ';
337+
std::cout<<ucell.atoms[i].taud[j].y<<' ';
338+
std::cout<<ucell.atoms[i].taud[j].z<<' ';
339+
}
340+
std::cout<<std::endl;
341+
}
342+
//ucell.update_pos_taud(move_ion);
343+
344+
std::cout<<"printtau"<<std::endl;
345+
for(int i=0;i<ucell.ntype;i++)
346+
{
347+
for(int j=0;j<ucell.atoms[i].na;j++)
348+
{
349+
std::cout<<ucell.atoms[i].tau[j].x<<' ';
350+
std::cout<<ucell.atoms[i].tau[j].y<<' ';
351+
std::cout<<ucell.atoms[i].tau[j].z<<' ';
352+
std::cout<<ucell.atoms[i].taud[j].x<<' ';
353+
std::cout<<ucell.atoms[i].taud[j].y<<' ';
354+
std::cout<<ucell.atoms[i].taud[j].z<<' ';
355+
}
356+
std::cout<<std::endl;
357+
}
358+
//pos = this->MAddM(pos, dpos);*/
263359
}
264360

265361
void BFGS::IsRestrain(std::vector<std::vector<double>>& dpos)
@@ -284,9 +380,48 @@ void BFGS::IsRestrain(std::vector<std::vector<double>>& dpos)
284380
}
285381
}
286382
}
287-
//std::cout<<a<<std::endl;
383+
std::cout<<"max dpos"<<std::endl;
384+
std::cout<<a<<std::endl;
288385
Ions_Move_Basic::converged = a<0.00001;
289386
}
387+
388+
void BFGS::CalculateLargestGrad(ModuleBase::matrix _force,UnitCell& ucell)
389+
{
390+
std::vector<double> grad= std::vector<double>(3*size, 0.0);
391+
int iat = 0;
392+
for (int it = 0; it < ucell.ntype; it++)
393+
{
394+
Atom *atom = &ucell.atoms[it];
395+
for (int ia = 0; ia < ucell.atoms[it].na; ia++)
396+
{
397+
for (int ik = 0; ik < 3; ++ik)
398+
{
399+
if (atom->mbl[ia][ik])
400+
{
401+
grad[3 * iat + ik] = -_force(iat, ik) * ucell.lat0;
402+
}
403+
}
404+
++iat;
405+
}
406+
}
407+
Ions_Move_Basic::largest_grad = 0.0;
408+
for (int i = 0; i < 3*size; i++)
409+
{
410+
if (Ions_Move_Basic::largest_grad < std::abs(grad[i]))
411+
{
412+
Ions_Move_Basic::largest_grad = std::abs(grad[i]);
413+
}
414+
}
415+
Ions_Move_Basic::largest_grad /= ucell.lat0;
416+
if (PARAM.inp.out_level == "ie")
417+
{
418+
std::cout << " LARGEST GRAD (eV/A) : " << Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / 0.529177
419+
<< std::endl;
420+
}
421+
422+
}
423+
424+
290425
// matrix methods
291426

292427
std::vector<double> BFGS::ReshapeMToV(std::vector<std::vector<double>> matrix)

0 commit comments

Comments
 (0)