bfgs_trad

Author: 19hello

source/module_relax/relax_new/bfgs.cpp

@@ -6,7 +6,11 @@
 void BFGS::init_relax(const int _size,UnitCell& ucell,double _maxstep) // initialize H0、H、pos0、force0、force
 {
     alpha=70;//relax_scale_force
-    maxstep=_maxstep; 
+    maxstep=_maxstep;
+    if(maxstep==0)
+    {
+        maxstep=0.8;
+    } 
     size=_size;
     sign =true;
     H = std::vector<std::vector<double>>(3*size, std::vector<double>(3*size, 0.0));
@@ -40,18 +44,19 @@ bool BFGS::relax_step(ModuleBase::matrix _force,UnitCell& ucell)
     //std::cout<<"enter Step0"<<std::endl;
     //std::cout<<size<<std::endl;
     //GetPos(ucell);
-    this->PrepareStep(force,pos,H,pos0,force0,steplength);
+    this->PrepareStep(force,pos,H,pos0,force0,steplength,dpos);
     //std::cout<<"enter Step1"<<std::endl;
     this->DetermineStep(steplength,dpos,maxstep);
+
     std::cout<<"enter Step2"<<std::endl;
-    for(int i=0;i<size;i++)
+    /*for(int i=0;i<size;i++)
     {
         for(int j=0;j<3;j++)
         {
             std::cout<<dpos[i][j]<<' ';
         }
         std::cout<<std::endl;
-    }
+    }*/
     this->UpdatePos(ucell);
     //std::cout<<"enter Step3"<<std::endl;
     return this->IsRestrain(dpos);
@@ -73,7 +78,7 @@ void BFGS::GetPos(UnitCell& ucell,std::vector<std::vector<double>>& pos)
 }
 
 
-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)
+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)
 {
     //std::cout<<"enter prepareStep0"<<std::endl;
     std::vector<double> changedforce = this->ReshapeMToV(force);
@@ -124,13 +129,29 @@ void BFGS::PrepareStep(std::vector<std::vector<double>>& force,std::vector<std::
     std::vector<double> a=this->DotInMAndV2(V, changedforce);
     for(int i = 0; i < a.size(); i++)
     {
-        a[i] /= abs(omega[i]);
+        if(omega[i]>0)
+        {
+            a[i] /= omega[i];
+        }
+        else if(omega[i]<0)
+        {
+           a[i] /= (-omega[i]);
+        }
+        
     }
 
     std::vector<double> tmpdpos = this->DotInMAndV1(V, a);
 
 
     dpos = this->ReshapeVToM(tmpdpos);
+    /*for(int i=0;i<size;i++)
+    {
+        for(int j=0;j<3;j++)
+        {
+            std::cout<<dpos[i][j]<<' ';
+        }
+        std::cout<<std::endl;
+    }*/
     for(int i = 0; i < size; i++)
     {
         double k = 0;
@@ -176,7 +197,7 @@ void BFGS::Update(std::vector<double> pos, std::vector<double> force,std::vector
     H = this->MSubM(H, this->MPlus(this->OuterVAndV(dg, dg), b));
 }
 
-void BFGS::DetermineStep(std::vector<double> steplength,std::vector<std::vector<double>>& dpos,int maxstep)
+void BFGS::DetermineStep(std::vector<double> steplength,std::vector<std::vector<double>>& dpos,double maxstep)
 {
     auto maxsteplength = max_element(steplength.begin(), steplength.end());
     double a = *maxsteplength;

source/module_relax/relax_new/bfgs.h

@@ -30,7 +30,7 @@ class BFGS
 
     void init_relax(const int _size,UnitCell& ucell,double _maxstep);//initialize parameters
     bool relax_step(ModuleBase::matrix _force,UnitCell& ucell);//
-    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);
+    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);
     bool IsRestrain(std::vector<std::vector<double>>& dpos);
 
 private:
@@ -39,7 +39,7 @@ class BFGS
 
     void GetPos(UnitCell& ucell,std::vector<std::vector<double>>& pos);
     void Update(std::vector<double> pos, std::vector<double> force,std::vector<std::vector<double>>& H);
-    void DetermineStep(std::vector<double> steplength,std::vector<std::vector<double>>& dpos,int maxstep);
+    void DetermineStep(std::vector<double> steplength,std::vector<std::vector<double>>& dpos,double maxstep);
     void UpdatePos(UnitCell& ucell);
 
 

tests/integrate/109_PW_CR/INPUT

@@ -4,12 +4,12 @@ suffix              autotest
 pseudo_dir		../../PP_ORB
 
 nbands			8
-calculation             cell_relax
+calculation             cell-relax
 #Parameters (Accuracy)
 ecutwfc			20
 scf_nmax			20
 kpar  2
-nspin 1
+nspin 2
 
 basis_type		pw
 relax_nmax                   2
@@ -20,7 +20,6 @@ stress_thr              1e-6
 cal_force                   1
 force_thr_ev            1.0e-3
 
-
 #ks_solver               scalapack_gvx
 mixing_type         broyden
 mixing_beta         0.7

tests/integrate/109_PW_CR/KPT

@@ -1,4 +1,4 @@
 K_POINTS
 0
 Gamma
-2 2 2 0 0 0
+3 1 2 3 1 2