remove vectorvector in bfgs and add parameter cp2k

docs/advanced/input_files/input-main.md

@@ -113,6 +113,7 @@
     - [npart\_sto](#npart_sto)
   - [Geometry Relaxation](#geometry-relaxation)
     - [relax\_method](#relax_method)
+    - [cp2k](#cp2k)
     - [relax\_new](#relax_new)
     - [relax\_scale\_force](#relax_scale_force)
     - [relax\_nmax](#relax_nmax)

source/source_io/input_conv.cpp

@@ -190,6 +190,7 @@ void Input_Conv::Convert()
     Ions_Move_Basic::relax_bfgs_init = PARAM.inp.relax_bfgs_init;
     Ions_Move_Basic::out_stru = PARAM.inp.out_stru; // mohan add 2012-03-23
     Ions_Move_Basic::relax_method = PARAM.inp.relax_method;
+    Ions_Move_Basic::cp2k = PARAM.inp.cp2k;
     Lattice_Change_Basic::fixed_axes = PARAM.inp.fixed_axes;
 
 

source/source_io/module_parameter/input_parameter.h

@@ -152,6 +152,7 @@ struct Input_para
 
     // ==============   #Parameters (4.Relaxation) ===========================
     std::vector<std::string> relax_method = {"cg","1"}; ///< methods to move_ion: sd, bfgs, cg...
+    std::string cp2k = "none";                    
     bool relax_new = true;
     bool relax = false; ///< allow relaxation along the specific direction
     double relax_scale_force = 0.5;

source/source_io/read_input_item_relax.cpp

@@ -48,6 +48,20 @@ void ReadInput::item_relax()
         // };
         // this->add_item(item);
     }
+    {
+        Input_Item item("cp2k");
+        item.annotation = "x; y; z; none;";
+        read_sync_string(input.cp2k);
+        item.check_value = [](const Input_Item& item, const Parameter& para) {
+            const std::vector<std::string> cp2ks = {"x", "y", "z", "none"};
+            if (std::find(cp2ks.begin(),cp2ks.end(), para.input.cp2k)==cp2ks.end())
+            {
+                const std::string warningstr = nofound_str(cp2ks, "cp2k");
+                ModuleBase::WARNING_QUIT("ReadInput", warningstr);
+            }
+        };
+        this->add_item(item);
+    }
     {
         Input_Item item("relax_new");
         item.annotation = "whether to use the new relaxation method";

source/source_relax/bfgs.cpp

@@ -15,11 +15,11 @@ void BFGS::allocate(const int _size)
     size=_size;
     largest_grad=0.0;
     sign=true;
-    H = std::vector<std::vector<double>>(3*size, std::vector<double>(3*size, 0.0));
+    H = ModuleBase::matrix(3*size, 3*size);
 
     for (int i = 0; i < 3*size; ++i) 
     {
-        H[i][i] = alpha;  
+        H(i,i) = alpha;  
     }
 
     pos = std::vector<ModuleBase::Vector3<double>> (size, ModuleBase::Vector3<double>(0.0, 0.0, 0.0)); 
@@ -114,7 +114,7 @@ void BFGS::GetPostaud(UnitCell& ucell,
 
 void BFGS::PrepareStep(std::vector<ModuleBase::Vector3<double>>& force,
                        std::vector<ModuleBase::Vector3<double>>& pos,
-                       std::vector<std::vector<double>>& H,
+                       ModuleBase::matrix& H,
                        std::vector<double>& pos0,
                        std::vector<double>& force0,
                        std::vector<double>& steplength,
@@ -133,20 +133,23 @@ void BFGS::PrepareStep(std::vector<ModuleBase::Vector3<double>>& force,
     int info=0;
     std::vector<double> H_flat;
 
-    for(const auto& row : H)
+    for(int i=0;i<H.nr;i++)
     {
-        H_flat.insert(H_flat.end(), row.begin(), row.end());
-    }
+        for(int j=0;j<H.nc;j++)
+        {
+            H_flat.push_back(H(i,j));
+        }
+    } 
 
     int value=3*size;
     int* ptr=&value;
     dsyev_("V","U",ptr,H_flat.data(),ptr,omega.data(),work.data(),&lwork,&info);
-    std::vector<std::vector<double>> V(3*size, std::vector<double>(3*size, 0.0));
+    ModuleBase::matrix V(3*size, 3*size);
     for(int i = 0; i < 3*size; i++)
     {
         for(int j = 0; j < 3*size; j++)
         {
-            V[j][i] = H_flat[3*size*i + j];
+            V(j,i) = H_flat[3*size*i + j];
         }
     }
     std::vector<double> a=DotInMAndV2(V, changedforce);
@@ -174,7 +177,7 @@ void BFGS::PrepareStep(std::vector<ModuleBase::Vector3<double>>& force,
 
 void BFGS::Update(std::vector<double>& pos, 
                   std::vector<double>& force,
-                  std::vector<std::vector<double>>& H,
+                  ModuleBase::matrix& H,
                   UnitCell& ucell)
 {
     if(sign)
@@ -240,10 +243,10 @@ void BFGS::Update(std::vector<double>& pos,
     double a = DotInVAndV(dpos, dforce);
     std::vector<double> dg = DotInMAndV1(H, dpos);
     double b = DotInVAndV(dpos, dg);
-    std::vector<std::vector<double>> term1=OuterVAndV(dforce, dforce);
-    std::vector<std::vector<double>> term2=OuterVAndV(dg, dg);
-    std::vector<std::vector<double>> term3=MPlus(term1, a);
-    std::vector<std::vector<double>> term4=MPlus(term2, b);
+    ModuleBase::matrix term1=OuterVAndV(dforce, dforce);
+    ModuleBase::matrix term2=OuterVAndV(dg, dg);
+    ModuleBase::matrix term3=MPlus(term1, a);
+    ModuleBase::matrix term4=MPlus(term2, b);
     H = MSubM(H, term3);
     H = MSubM(H, term4);
 }

source/source_relax/bfgs.h

@@ -26,7 +26,7 @@ class BFGS
     int size;//number of atoms
 
     std::vector<double> steplength;//the length of atoms displacement 
-    std::vector<std::vector<double>> H;//Hessian matrix
+    ModuleBase::matrix H;//Hessian matrix
     std::vector<double> force0;//force in previous step
     std::vector<ModuleBase::Vector3<double>> force;
     std::vector<double> pos0;//atom pos in previous step(cartesian coordinates)
@@ -35,12 +35,12 @@ class BFGS
     std::vector<ModuleBase::Vector3<double>> pos_taud;
     std::vector<ModuleBase::Vector3<double>> dpos;
 
-    void PrepareStep(std::vector<ModuleBase::Vector3<double>>& force,std::vector<ModuleBase::Vector3<double>>& pos,std::vector<std::vector<double>>& H,std::vector<double>& pos0,std::vector<double>& force0,std::vector<double>& steplength,std::vector<ModuleBase::Vector3<double>>& dpos,int& size,UnitCell& ucell);//calculate the atomic displacement in one iteration step
+    void PrepareStep(std::vector<ModuleBase::Vector3<double>>& force,std::vector<ModuleBase::Vector3<double>>& pos,ModuleBase::matrix& H,std::vector<double>& pos0,std::vector<double>& force0,std::vector<double>& steplength,std::vector<ModuleBase::Vector3<double>>& dpos,int& size,UnitCell& ucell);//calculate the atomic displacement in one iteration step
     void IsRestrain();//check if converged
     void CalculateLargestGrad(const ModuleBase::matrix& _force,UnitCell& ucell);
     void GetPos(UnitCell& ucell,std::vector<ModuleBase::Vector3<double>>& pos);
     void GetPostaud(UnitCell& ucell,std::vector<ModuleBase::Vector3<double>>& pos_taud);
-    void Update(std::vector<double>& pos, std::vector<double>& force,std::vector<std::vector<double>>& H,UnitCell& ucell);//update hessian matrix
+    void Update(std::vector<double>& pos, std::vector<double>& force,ModuleBase::matrix& H,UnitCell& ucell);//update hessian matrix
     void DetermineStep(std::vector<double>& steplength,std::vector<ModuleBase::Vector3<double>>& dpos,double& maxstep);//normalize large atomic displacements based on maxstep
     void UpdatePos(UnitCell& ucell);//update ucell with the new coordinates
 

source/source_relax/ions_move_basic.h

@@ -22,6 +22,7 @@ extern double relax_bfgs_rmin;  // min value of trust radius,
 extern double relax_bfgs_init;  // initial value of trust radius,
 extern double best_xxx;         // the last step length of cg , we use it as  bfgs`s initial step length
 extern std::vector<std::string> relax_method; // relaxation method,
+extern std::string cp2k;
 extern int out_stru; // output the structure or not
 // funny way to pass this parameter, but nevertheless
 

source/source_relax/lattice_change_basic.cpp

@@ -16,6 +16,8 @@ double Lattice_Change_Basic::ediff = 0.0;
 double Lattice_Change_Basic::etot = 0.0;
 double Lattice_Change_Basic::etot_p = 0.0;
 
+std::string Ions_Move_Basic::cp2k = "none";
+
 // double Lattice_Change_Basic::lattice_change_ini = 0.5; // default is 0.5
 double Lattice_Change_Basic::lattice_change_ini = 0.01; // default is 0.5
 std::string Lattice_Change_Basic::fixed_axes = "None";
@@ -83,6 +85,18 @@ void Lattice_Change_Basic::change_lattice(UnitCell &ucell, double *move, double
 			for (int j = 0;j < 3;++j) 
 			{
 				move_mat_t(j, i) = move[i * 3 + j] / ucell.lat0;    //transpose
+                if(j==0&&Ions_Move_Basic::cp2k=="x")
+                {
+                    move_mat_t(j, i) = 0.0;
+                }
+                if(j==1&&Ions_Move_Basic::cp2k=="y")
+                {
+                    move_mat_t(j, i) = 0.0; 
+                }
+                if (j==2&&Ions_Move_Basic::cp2k=="z")
+                {
+                    move_mat_t(j, i) = 0.0;
+                }              
 			}
 		}
 		ModuleBase::matrix symm_move_mat_t = (move_mat_t * ucell.G.to_matrix());//symmetrize (latvec^{-1} * move_mat)^T
@@ -97,7 +111,6 @@ void Lattice_Change_Basic::change_lattice(UnitCell &ucell, double *move, double
 			}
 		}
     }
-
     if (ucell.lc[0] != 0)
     {
         ucell.latvec.e11 = (move[0] + lat[0]) / ucell.lat0;

source/source_relax/lattice_change_basic.h

@@ -3,6 +3,7 @@
 
 #include "source_base/matrix.h"
 #include "source_cell/unitcell.h"
+#include "ions_move_basic.h"
 
 namespace Lattice_Change_Basic
 {

source/source_relax/lbfgs.cpp

@@ -12,7 +12,6 @@ void LBFGS::allocate(const int _size) // initialize H0、H、pos0、force0、for
     size=_size;
     memory=100;
     iteration=0;
-    H = std::vector<std::vector<double>>(3*size, std::vector<double>(3*size, 0.0));
     H0=1/alpha;
     pos = std::vector<ModuleBase::Vector3<double>> (size, ModuleBase::Vector3<double>(0.0, 0.0, 0.0)); 
     pos0 = std::vector<double>(3*size, 0.0);
@@ -59,7 +58,7 @@ void LBFGS::relax_step(const ModuleBase::matrix _force,UnitCell& ucell,const dou
         }
         k+=ucell.atoms[i].na;
     }
-    this->prepare_step(force,pos,H,pos0,force0,dpos,ucell,etot);
+    this->prepare_step(force,pos,pos0,force0,dpos,ucell,etot);
     this->determine_step(steplength,dpos,maxstep);
     this->update_pos(ucell);
     this->calculate_largest_grad(_force,ucell);
@@ -100,7 +99,6 @@ void LBFGS::get_pos_taud(UnitCell& ucell,std::vector<ModuleBase::Vector3<double>
 
 void LBFGS::prepare_step(std::vector<ModuleBase::Vector3<double>>& force,
                          std::vector<ModuleBase::Vector3<double>>& pos,
-                         std::vector<std::vector<double>>& H,
                          std::vector<double>& pos0,
                          std::vector<double>& force0,
                          std::vector<ModuleBase::Vector3<double>>& dpos,

source/source_relax/lbfgs.h

@@ -50,7 +50,6 @@ class LBFGS
 
     ModuleESolver::ESolver* solver;         ///< Structure solver
     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
     std::vector<ModuleBase::Vector3<double>> force;
     std::vector<double> pos0;//atom pos in previous step(cartesian coordinates)
@@ -67,7 +66,6 @@ class LBFGS
      */
     void prepare_step(std::vector<ModuleBase::Vector3<double>>& force,
                       std::vector<ModuleBase::Vector3<double>>& pos,
-                      std::vector<std::vector<double>>& H,
                       std::vector<double>& pos0,
                       std::vector<double>& force0,
                       std::vector<ModuleBase::Vector3<double>>& dpos,

source/source_relax/matrix_methods.cpp

@@ -18,17 +18,17 @@ std::vector<double> ReshapeMToV(std::vector<ModuleBase::Vector3<double>>& matrix
     return result;
 }
 
-std::vector<std::vector<double>> MAddM(std::vector<std::vector<double>>& a, 
-                                             std::vector<std::vector<double>>& b) 
+ModuleBase::matrix MAddM(ModuleBase::matrix& a, 
+                                             ModuleBase::matrix& 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++)
+    assert(a.nr!=0 && a.nc!=0 && b.nr!=0 && b.nc!=0);
+    assert(a.nr == b.nr && a.nc == b.nc);
+    ModuleBase::matrix result = ModuleBase::matrix(a.nr, a.nc);
+    for(int i = 0; i < a.nr; i++)
     {
-        for(int j = 0; j < a[0].size(); j++)
+        for(int j = 0; j < a.nc; j++)
         {
-            result[i][j] = a[i][j] + b[i][j];
+            result(i,j) = a(i,j) + b(i,j);
         }
     }
     return result;
@@ -58,30 +58,30 @@ std::vector<ModuleBase::Vector3<double>> ReshapeVToM(std::vector<double>& matrix
     return result;
 }
 
-std::vector<double> DotInMAndV1(std::vector<std::vector<double>>& matrix, std::vector<double>& vec) 
+std::vector<double> DotInMAndV1(ModuleBase::matrix& matrix, std::vector<double>& vec) 
 {
-    assert(!matrix.empty());
-    assert(matrix[0].size() == vec.size());
-    std::vector<double> result(matrix.size(), 0.0);
+    assert(matrix.nr!=0 && matrix.nc!=0);
+    assert(matrix.nc == vec.size());
+    std::vector<double> result(matrix.nr, 0.0);
     for(int i = 0; i < result.size(); i++)
     {
         for(int j = 0; j < vec.size(); j++)
         {
-            result[i] += matrix[i][j] * vec[j];
+            result[i] += matrix(i,j) * vec[j];
         }
     }
     return result;
 }
-std::vector<double> DotInMAndV2(std::vector<std::vector<double>>& matrix, std::vector<double>& vec) 
+std::vector<double> DotInMAndV2(ModuleBase::matrix& matrix, std::vector<double>& vec) 
 {
-    assert(!matrix.empty());
-    assert(matrix.size() == vec.size());
-    std::vector<double> result(matrix.size(), 0.0);
+    assert(matrix.nr!=0 && matrix.nc!=0);
+    assert(matrix.nr == vec.size());
+    std::vector<double> result(matrix.nc, 0.0);
     for(int i = 0; i < result.size(); i++)
     {
         for(int j = 0; j < vec.size(); j++)
         {
-            result[i] += matrix[j][i] * vec[j];
+            result[i] += matrix(j,i) * vec[j];
         }
     }
     return result;
@@ -98,53 +98,52 @@ double DotInVAndV(std::vector<double>& vec1, std::vector<double>& vec2)
     return result;
 }
 
-std::vector<std::vector<double>> OuterVAndV(std::vector<double>& a, std::vector<double>& b) 
+ModuleBase::matrix 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));
+    ModuleBase::matrix result = ModuleBase::matrix(a.size(), b.size());
     for(int i = 0; i < a.size(); i++)
     {
         for(int j = 0; j < b.size(); j++)
         {
-            result[i][j] = a[i] * b[j];
+            result(i,j) = a[i] * b[j];
         }
     }
     return result;
 }
 
-std::vector<std::vector<double>> MPlus(std::vector<std::vector<double>>& a, double b)
+ModuleBase::matrix MPlus(ModuleBase::matrix& a, double b)
 {
-    assert(!a.empty());
+    assert(a.nr!=0 && a.nc!=0);
     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++)
+    ModuleBase::matrix result = ModuleBase::matrix(a.nr, a.nc);
+    for(int i = 0; i < a.nr; i++)
     {
-        for(int j = 0; j < a[0].size(); j++)
+        for(int j = 0; j < a.nc; j++)
         {
-            result[i][j] = a[i][j] / b;
+            result(i,j) = a(i,j) / b;
         }
     }
     return result;
 }
 
-std::vector<std::vector<double>> MSubM(std::vector<std::vector<double>>& a, std::vector<std::vector<double>>& b)
+ModuleBase::matrix MSubM(ModuleBase::matrix& a, ModuleBase::matrix& 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++)
+    assert(a.nr!=0 && a.nc!=0 && b.nr!=0 && b.nc!=0);
+    assert(a.nr == b.nr && a.nc == b.nc);
+    ModuleBase::matrix result = ModuleBase::matrix(a.nr, a.nc);
+    for(int i = 0; i < a.nr; i++)
     {
-        for(int j = 0; j < a[0].size(); j++)
+        for(int j = 0; j < a.nc; j++)
         {
-            result[i][j] = a[i][j] - b[i][j];
+            result(i,j) = a(i,j) - b(i,j);
         }
     }
     return result;
 }
 
 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++)
     {