22#include " module_hamilt_pw/hamilt_pwdft/global.h"
33#include " module_base/matrix3.h"
44#include " module_relax/relax_old/ions_move_basic.h"
5+ #include " module_parameter/parameter.h"
6+ #include " ions_move_basic.h"
57
6- void BFGS::init_relax (const int _size,UnitCell& ucell,double _maxstep) // initialize H0、H、pos0、force0、force
8+
9+ void BFGS::allocate (const int _size) // initialize H0、H、pos0、force0、force
710{
811 alpha=70 ;// relax_scale_force
9- maxstep=_maxstep ;
12+ maxstep=PARAM. inp . relax_bfgs_rmax ;
1013 if (maxstep==0 )
1114 {
1215 maxstep=0.8 ;
@@ -23,13 +26,19 @@ void BFGS::init_relax(const int _size,UnitCell& ucell,double _maxstep) // initia
2326 force0 = std::vector<double >(3 *size, 0.0 );
2427 force = std::vector<std::vector<double >>(size, std::vector<double >(3 , 0.0 ));
2528 steplength = std::vector<double >(size, 0.0 );
26- this ->GetPos (ucell,pos);
29+ this ->GetPos (GlobalC::ucell,pos);
30+
2731
2832}
2933
30- bool BFGS::relax_step (ModuleBase::matrix _force,UnitCell& ucell)
34+ void BFGS::relax_step (ModuleBase::matrix _force,UnitCell& ucell)
3135{
3236 std::cout<<" enter Step" <<std::endl;
37+ if (sign)
38+ {
39+
40+ }
41+ std::cout<<" enter Step1" <<std::endl;
3342 ucell.ionic_position_updated = true ;
3443 for (int i = 0 ; i < _force.nr ; i++)
3544 {
@@ -59,21 +68,21 @@ bool BFGS::relax_step(ModuleBase::matrix _force,UnitCell& ucell)
5968 }*/
6069 this ->UpdatePos (ucell);
6170 // std::cout<<"enter Step3"<<std::endl;
62- return this ->IsRestrain (dpos);
71+ this ->IsRestrain (dpos);
6372}
6473
6574void BFGS::GetPos (UnitCell& ucell,std::vector<std::vector<double >>& pos)
6675{
6776 int k=0 ;
68- for (int i=0 ;i<GlobalC:: ucell.ntype ;i++)
77+ for (int i=0 ;i<ucell.ntype ;i++)
6978 {
70- for (int j=0 ;j<GlobalC:: ucell.atoms [i].na ;j++)
79+ for (int j=0 ;j<ucell.atoms [i].na ;j++)
7180 {
72- pos[k+j][0 ]=GlobalC:: ucell.atoms [i].tau [j].x *ModuleBase::BOHR_TO_A*GlobalC:: ucell.lat0 ;
73- pos[k+j][1 ]=GlobalC:: ucell.atoms [i].tau [j].y *ModuleBase::BOHR_TO_A*GlobalC:: ucell.lat0 ;
74- pos[k+j][2 ]=GlobalC:: ucell.atoms [i].tau [j].z *ModuleBase::BOHR_TO_A*GlobalC:: ucell.lat0 ;
81+ pos[k+j][0 ]=ucell.atoms [i].tau [j].x *ModuleBase::BOHR_TO_A*ucell.lat0 ;
82+ pos[k+j][1 ]=ucell.atoms [i].tau [j].y *ModuleBase::BOHR_TO_A*ucell.lat0 ;
83+ pos[k+j][2 ]=ucell.atoms [i].tau [j].z *ModuleBase::BOHR_TO_A*ucell.lat0 ;
7584 }
76- k+=GlobalC:: ucell.atoms [i].na ;
85+ k+=ucell.atoms [i].na ;
7786 }
7887}
7988
@@ -86,14 +95,14 @@ void BFGS::PrepareStep(std::vector<std::vector<double>>& force,std::vector<std::
8695 // std::cout<<"enter prepareStep1"<<std::endl;
8796 this ->Update (changedpos, changedforce,H);
8897 // std::cout<<"enter prepareStep2"<<std::endl;
89- for (int i = 0 ; i < 3 *size; i++)
98+ /* for(int i = 0; i < 3*size; i++)
9099 {
91100 for(int j = 0; j < 3*size; j++)
92101 {
93102 std::cout<<H[i][j]<<' ';
94103 }
95104 std::cout<<std::endl;
96- }
105+ }*/
97106 // call dysev
98107 // std::cout<<size<<std::endl;
99108 std::vector<double > omega (3 *size);
@@ -226,20 +235,20 @@ void BFGS::UpdatePos(UnitCell& ucell)
226235
227236 // convert unit
228237 ModuleBase::Vector3<double > move_ion_cart;
229- move_ion_cart.x = dpos[iat][0 ] / ModuleBase::BOHR_TO_A / GlobalC:: ucell.lat0 ;
230- move_ion_cart.y = dpos[iat][1 ] / ModuleBase::BOHR_TO_A / GlobalC:: ucell.lat0 ;
231- move_ion_cart.z = dpos[iat][2 ] / ModuleBase::BOHR_TO_A / GlobalC:: ucell.lat0 ;
238+ move_ion_cart.x = dpos[iat][0 ] / ModuleBase::BOHR_TO_A / ucell.lat0 ;
239+ move_ion_cart.y = dpos[iat][1 ] / ModuleBase::BOHR_TO_A / ucell.lat0 ;
240+ move_ion_cart.z = dpos[iat][2 ] / ModuleBase::BOHR_TO_A / ucell.lat0 ;
232241
233242 // convert to Direct coordinate
234243 // note here the old GT is used
235244
236245 // convert pos
237- ModuleBase::Vector3<double > move_ion_dr = move_ion_cart * GlobalC:: ucell.GT ;
246+ ModuleBase::Vector3<double > move_ion_dr = move_ion_cart * ucell.GT ;
238247
239248
240- int it = GlobalC:: ucell.iat2it [iat];
241- int ia = GlobalC:: ucell.iat2ia [iat];
242- Atom* atom = &GlobalC:: ucell.atoms [it];
249+ int it = ucell.iat2it [iat];
250+ int ia = ucell.iat2ia [iat];
251+ Atom* atom = &ucell.atoms [it];
243252
244253 if (atom->mbl [ia].x == 1 )
245254 {
@@ -254,11 +263,11 @@ void BFGS::UpdatePos(UnitCell& ucell)
254263 move_ion[iat * 3 + 2 ] = move_ion_dr.z ;
255264 }
256265 }
257- GlobalC:: ucell.update_pos_taud (move_ion);
266+ ucell.update_pos_taud (move_ion);
258267 pos = this ->MAddM (pos, dpos);
259268}
260269
261- bool BFGS::IsRestrain (std::vector<std::vector<double >>& dpos)
270+ void BFGS::IsRestrain (std::vector<std::vector<double >>& dpos)
262271{
263272 double a=0 ;
264273 for (int i=0 ;i<size;i++)
@@ -281,7 +290,7 @@ bool BFGS::IsRestrain(std::vector<std::vector<double>>& dpos)
281290 }
282291 }
283292 std::cout<<a<<std::endl;
284- return a<0.00001 ;
293+ Ions_Move_Basic::converged = a<0.00001 ;
285294}
286295// matrix methods
287296
0 commit comments