@@ -13,7 +13,6 @@ void bfgs::initialize(int _size) // initialize H0、H、pos0、force0、force
1313 for (int i = 0 ; i < 3 *size; ++i) {
1414 H[i][i] = alpha;
1515 }
16-
1716 pos = std::vector<std::vector<double >> (size, std::vector<double >(3 , 0.0 ));
1817 pos0 = std::vector<double >(3 *size, 0.0 );
1918 dpos = std::vector<std::vector<double >>(size, std::vector<double >(3 , 0.0 ));
@@ -29,10 +28,13 @@ bool bfgs::Step(std::vector<std::vector<double>> _force,UnitCell& ucell)
2928 std::cout<<" enter Step" <<std::endl;
3029 GlobalC::ucell.ionic_position_updated = true ;
3130 force = _force;
32-
31+ // std::cout<<"enter Step0"<<std::endl;
32+ // std::cout<<size<<std::endl;
3333 // GetPos(ucell);
3434 PrepareStep ();
35+ // std::cout<<"enter Step1"<<std::endl;
3536 DetermineStep ();
37+ // std::cout<<"enter Step2"<<std::endl;
3638 /* for(int i=0;i<size;i++)
3739 {
3840 for(int j=0;j<3;j++)
@@ -42,6 +44,7 @@ bool bfgs::Step(std::vector<std::vector<double>> _force,UnitCell& ucell)
4244 std::cout<<std::endl;
4345 }*/
4446 UpdatePos ();
47+ // std::cout<<"enter Step3"<<std::endl;
4548 return IsRestrain ();
4649}
4750
@@ -63,9 +66,12 @@ void bfgs::GetPos()
6366
6467void bfgs::PrepareStep ()
6568{
69+ // std::cout<<"enter prepareStep0"<<std::endl;
6670 std::vector<double > changedforce = ReshapeMToV (force);
6771 std::vector<double > changedpos = ReshapeMToV (pos);
72+ // std::cout<<"enter prepareStep1"<<std::endl;
6873 Update (changedpos, changedforce);
74+ // std::cout<<"enter prepareStep2"<<std::endl;
6975 /* for(int i = 0; i < 3*size; i++)
7076 {
7177 for(int j = 0; j < 3*size; j++)
@@ -75,8 +81,11 @@ void bfgs::PrepareStep()
7581 std::cout<<std::endl;
7682 }*/
7783 // call dysev
84+ // std::cout<<size<<std::endl;
7885 std::vector<double > omega (3 *size);
86+ // std::cout<<"enter prepareStep3"<<std::endl;
7987 std::vector<double > work (3 *size*3 *size);
88+ // std::cout<<"enter prepareStep4"<<std::endl;
8089 int lwork=3 *size*3 *size;
8190 int info;
8291 std::vector<double > H_flat;
@@ -88,13 +97,15 @@ void bfgs::PrepareStep()
8897 int * ptr=&value;
8998 dsyev_ (" V" ," U" ,ptr,H_flat.data (),ptr,omega.data (),work.data (),&lwork,&info);
9099 std::vector<std::vector<double >> V (3 *size, std::vector<double >(3 *size, 0.0 ));
100+ // std::cout<<"enter prepareStep5"<<std::endl;
91101 for (int i = 0 ; i < 3 *size; i++)
92102 {
93103 for (int j = 0 ; j < 3 *size; j++)
94104 {
95105 V[j][i] = H_flat[3 *size*i + j];
96106 }
97107 }
108+
98109
99110 /* for(int i=0;i<3*size;i++)
100111 {
0 commit comments