Modify BFGS test code

source/source_relax/bfgs.cpp

@@ -10,10 +10,12 @@
 //! initialize H0、H、pos0、force0、force
 void BFGS::allocate(const int _size) 
 {
+    assert(_size > 0);
     alpha=70;//default value in ase is 70
     maxstep=PARAM.inp.relax_bfgs_rmax;
     size=_size;
-    sign =true;
+    largest_grad=0.0;
+    sign=true;
     H = std::vector<std::vector<double>>(3*size, std::vector<double>(3*size, 0.0));
 
     for (int i = 0; i < 3*size; ++i) 
@@ -37,6 +39,7 @@ void BFGS::relax_step(const ModuleBase::matrix& _force,UnitCell& ucell)
     GetPos(ucell,pos);  
     GetPostaud(ucell,pos_taud);
     ucell.ionic_position_updated = true;
+    assert(_force.nr == force.size() && _force.nc == force[0].size());
     for(int i = 0; i < _force.nr; i++)
     {
         for(int j=0;j<_force.nc;j++)
@@ -65,7 +68,7 @@ void BFGS::relax_step(const ModuleBase::matrix& _force,UnitCell& ucell)
         k+=ucell.atoms[i].na;
     }
 
-    this->PrepareStep(force,pos,H,pos0,force0,steplength,dpos,ucell);
+    this->PrepareStep(force,pos,H,pos0,force0,steplength,dpos,size,ucell);
     this->DetermineStep(steplength,dpos,maxstep);
     this->UpdatePos(ucell);
     this->CalculateLargestGrad(_force,ucell);
@@ -76,6 +79,7 @@ void BFGS::relax_step(const ModuleBase::matrix& _force,UnitCell& ucell)
 
 void BFGS::GetPos(UnitCell& ucell,std::vector<std::vector<double>>& pos)
 {
+    assert(pos.size() == ucell.nat);
     int k=0;
     for(int i=0;i<ucell.ntype;i++)
     {
@@ -92,6 +96,7 @@ void BFGS::GetPos(UnitCell& ucell,std::vector<std::vector<double>>& pos)
 void BFGS::GetPostaud(UnitCell& ucell,
                       std::vector<std::vector<double>>& pos_taud)
 {
+    assert(pos_taud.size() == ucell.nat);
     int k=0;
     for(int i=0;i<ucell.ntype;i++)
     {
@@ -112,6 +117,7 @@ void BFGS::PrepareStep(std::vector<std::vector<double>>& force,
                        std::vector<double>& force0,
                        std::vector<double>& steplength,
                        std::vector<std::vector<double>>& dpos,
+                       int& size,
                        UnitCell& ucell)
 {
     std::vector<double> changedforce = ReshapeMToV(force);
@@ -142,8 +148,10 @@ void BFGS::PrepareStep(std::vector<std::vector<double>>& force,
         }
     }
     std::vector<double> a=DotInMAndV2(V, changedforce);
+    double threshold=1e-8;
     for(int i = 0; i < a.size(); i++)
     {
+        assert(std::abs(omega[i]) > threshold);
         a[i]/=std::abs(omega[i]);    
     }
     std::vector<double> tmpdpos = DotInMAndV1(V, a);
@@ -221,7 +229,8 @@ void BFGS::Update(std::vector<double>& pos,
             dpos[iat * 3 + 2] = move_ion_dr.z ;
         }
     }
-    if(*max_element(dpos.begin(), dpos.end()) < 1e-7)
+    double threshold=1e-7;
+    if(*max_element(dpos.begin(), dpos.end()) < threshold)
     {
         return;
     } 
@@ -243,6 +252,8 @@ void BFGS::DetermineStep(std::vector<double>& steplength,
 {
     std::vector<double>::iterator maxsteplength = max_element(steplength.begin(), steplength.end());
     double a = *maxsteplength;
+    double threshold=1e-10;
+    assert(a > threshold);
     if(a >= maxstep)
     {
         double scale = maxstep / a;
@@ -311,8 +322,8 @@ void BFGS::UpdatePos(UnitCell& ucell)
 
 void BFGS::IsRestrain(std::vector<std::vector<double>>& dpos)
 {
-    Ions_Move_Basic::converged = Ions_Move_Basic::largest_grad 
-        * ModuleBase::Ry_to_eV / 0.529177<PARAM.inp.force_thr_ev;
+    Ions_Move_Basic::converged = largest_grad 
+        * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A<PARAM.inp.force_thr_ev;
 }
 
 void BFGS::CalculateLargestGrad(const ModuleBase::matrix& _force,UnitCell& ucell)
@@ -334,20 +345,20 @@ void BFGS::CalculateLargestGrad(const ModuleBase::matrix& _force,UnitCell& ucell
             ++iat;
         }
     }
-    Ions_Move_Basic::largest_grad = 0.0;
+    largest_grad = 0.0;
     for (int i = 0; i < 3*size; i++)
     {
-        if (Ions_Move_Basic::largest_grad < std::abs(grad[i]))
+        if (largest_grad < std::abs(grad[i]))
         {
-            Ions_Move_Basic::largest_grad = std::abs(grad[i]);
+            largest_grad = std::abs(grad[i]);
         }
     }
+    assert(ucell.lat0 != 0);
     Ions_Move_Basic::largest_grad /= ucell.lat0;
     if (PARAM.inp.out_level == "ie")
     {
-        std::cout << " LARGEST GRAD (eV/Angstrom)  : " << Ions_Move_Basic::largest_grad 
-            * ModuleBase::Ry_to_eV / 0.5291772109
+        std::cout << " LARGEST GRAD (eV/Angstrom)  : " << largest_grad 
+            * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A
                   << std::endl;
     }
-
 }

source/source_relax/bfgs.h

@@ -13,6 +13,17 @@
 class BFGS
 {
 public:
+    void allocate(const int _size);//initialize parameters
+    void relax_step(const ModuleBase::matrix& _force,UnitCell& ucell);//a full iteration step
+
+
+private:
+    bool sign;//check if this is the first iteration
+    double alpha;//initialize H,diagonal element is alpha
+    double maxstep;//every movement smaller than maxstep
+    double largest_grad;
+    int size;//number of atoms
+
     std::vector<double> steplength;//the length of atoms displacement 
     std::vector<std::vector<double>> H;//Hessian matrix
     std::vector<double> force0;//force in previous step
@@ -22,17 +33,8 @@ class BFGS
     std::vector<double> pos_taud0;//atom pos in previous step(relative coordinates)
     std::vector<std::vector<double>> pos_taud;
     std::vector<std::vector<double>> dpos;
-
-    void allocate(const int _size);//initialize parameters
-    void relax_step(const ModuleBase::matrix& _force,UnitCell& ucell);//a full iteration step
-    void 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);//calculate the atomic displacement in one iteration step
-
-private:
-    bool sign;//check if this is the first iteration
-    double alpha;//initialize H,diagonal element is alpha
-    double maxstep;//every movement smaller than maxstep
-    int size;//number of atoms
 
+    void 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,int& size,UnitCell& ucell);//calculate the atomic displacement in one iteration step
     void IsRestrain(std::vector<std::vector<double>>& dpos);//check if converged
     void CalculateLargestGrad(const ModuleBase::matrix& _force,UnitCell& ucell);
     void GetPos(UnitCell& ucell,std::vector<std::vector<double>>& pos);

source/source_relax/matrix_methods.cpp

@@ -4,6 +4,8 @@
 
 std::vector<double> ReshapeMToV(std::vector<std::vector<double>>& matrix) 
 {
+    assert(!matrix.empty());
+    assert(matrix[0].size() == 3);
     int size = matrix.size();
     std::vector<double> result;
     result.reserve(3*size);
@@ -16,6 +18,8 @@ std::vector<double> ReshapeMToV(std::vector<std::vector<double>>& matrix)
 std::vector<std::vector<double>> MAddM(std::vector<std::vector<double>>& a, 
                                              std::vector<std::vector<double>>& b) 
 {
+    assert(!a.empty() && !b.empty());
+    assert(a.size() == b.size() && a[0].size() == b[0].size());
     std::vector<std::vector<double>> result = std::vector<std::vector<double>>(a.size(), std::vector<double>(a[0].size(), 0.0));
     for(int i = 0; i < a.size(); i++)
     {
@@ -29,6 +33,7 @@ std::vector<std::vector<double>> MAddM(std::vector<std::vector<double>>& a,
 
 std::vector<double> VSubV(std::vector<double>& a, std::vector<double>& b) 
 {
+    assert(a.size() == b.size());
     std::vector<double> result = std::vector<double>(a.size(), 0.0);
     for(int i = 0; i < a.size(); i++)
     {
@@ -39,6 +44,7 @@ std::vector<double> VSubV(std::vector<double>& a, std::vector<double>& b)
 
 std::vector<std::vector<double>> ReshapeVToM(std::vector<double>& matrix) 
 {
+    assert(matrix.size() % 3 == 0);
     std::vector<std::vector<double>> result = std::vector<std::vector<double>>(matrix.size() / 3, std::vector<double>(3));
     for(int i = 0; i < result.size(); i++)
     {
@@ -52,6 +58,8 @@ std::vector<std::vector<double>> ReshapeVToM(std::vector<double>& matrix)
 
 std::vector<double> DotInMAndV1(std::vector<std::vector<double>>& matrix, std::vector<double>& vec) 
 {
+    assert(!matrix.empty());
+    assert(matrix[0].size() == vec.size());
     std::vector<double> result(matrix.size(), 0.0);
     for(int i = 0; i < result.size(); i++)
     {
@@ -64,6 +72,8 @@ std::vector<double> DotInMAndV1(std::vector<std::vector<double>>& matrix, std::v
 }
 std::vector<double> DotInMAndV2(std::vector<std::vector<double>>& matrix, std::vector<double>& vec) 
 {
+    assert(!matrix.empty());
+    assert(matrix.size() == vec.size());
     std::vector<double> result(matrix.size(), 0.0);
     for(int i = 0; i < result.size(); i++)
     {
@@ -77,6 +87,7 @@ std::vector<double> DotInMAndV2(std::vector<std::vector<double>>& matrix, std::v
 
 double DotInVAndV(std::vector<double>& vec1, std::vector<double>& vec2) 
 {
+    assert(vec1.size() == vec2.size());
     double result = 0.0;
     for(int i = 0; i < vec1.size(); i++)
     {
@@ -87,6 +98,7 @@ double DotInVAndV(std::vector<double>& vec1, std::vector<double>& vec2)
 
 std::vector<std::vector<double>> OuterVAndV(std::vector<double>& a, std::vector<double>& b) 
 {
+    assert(a.size() == b.size());
     std::vector<std::vector<double>> result = std::vector<std::vector<double>>(a.size(), std::vector<double>(b.size(), 0.0));
     for(int i = 0; i < a.size(); i++)
     {
@@ -100,6 +112,8 @@ std::vector<std::vector<double>> OuterVAndV(std::vector<double>& a, std::vector<
 
 std::vector<std::vector<double>> MPlus(std::vector<std::vector<double>>& a, double b)
 {
+    assert(!a.empty());
+    assert(b != 0);
     std::vector<std::vector<double>> result = std::vector<std::vector<double>>(a.size(), std::vector<double>(a[0].size(), 0.0));
     for(int i = 0; i < a.size(); i++)
     {
@@ -113,6 +127,8 @@ std::vector<std::vector<double>> MPlus(std::vector<std::vector<double>>& a, doub
 
 std::vector<std::vector<double>> MSubM(std::vector<std::vector<double>>& a, std::vector<std::vector<double>>& b)
 {
+    assert(!a.empty() && !b.empty());
+    assert(a.size() == b.size() && a[0].size() == b[0].size());
     std::vector<std::vector<double>> result = std::vector<std::vector<double>>(a.size(), std::vector<double>(a[0].size(), 0.0));
     for(int i = 0; i < a.size(); i++)
     {
@@ -126,6 +142,7 @@ std::vector<std::vector<double>> MSubM(std::vector<std::vector<double>>& a, std:
 
 std::vector<double> DotInVAndFloat(std::vector<double>& vec, double b) 
 {
+    assert(b != 0);
     std::vector<double> result(vec.size(), 0.0);
     for(int i = 0; i < vec.size(); i++)
     {
@@ -136,6 +153,7 @@ std::vector<double> DotInVAndFloat(std::vector<double>& vec, double b)
 
 std::vector<double> VAddV(std::vector<double>& a, std::vector<double>& b) 
 {
+    assert(a.size() == b.size());
     std::vector<double> result = std::vector<double>(a.size(), 0.0);
     for(int i = 0; i < a.size(); i++)
     {

source/source_relax/matrix_methods.h

@@ -2,6 +2,7 @@
 #define MATRIX_METHODS
 
 #include <vector>
+#include <cassert>