diff --git a/source/module_relax/relax_new/line_search.cpp b/source/module_relax/relax_new/line_search.cpp
index 27772b4d10..138f36ca51 100644
--- a/source/module_relax/relax_new/line_search.cpp
+++ b/source/module_relax/relax_new/line_search.cpp
@@ -1,101 +1,96 @@
 #include "line_search.h"
-#include <cmath>
-#include <algorithm>
 
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
 
-bool Line_Search::line_search(
-    const bool restart,
-    const double x, //current point
-    const double y, //function value at current point
-    const double f, //derivative at current point
-    double & xnew,  //the next point that we want to try
-    const double conv_thr)
+#include <algorithm>
+#include <cmath>
+
+bool Line_Search::line_search(const bool restart,
+                              const double x, // current point
+                              const double y, // function value at current point
+                              const double f, // derivative at current point
+                              double& xnew,   // the next point that we want to try
+                              const double conv_thr)
 {
-    if(restart) ls_step = 0;
+    if (restart) {
+        ls_step = 0;
+}
 
-    if(ls_step == 0) //first point: make a trial step into trial direction
+    if (ls_step == 0) // first point: make a trial step into trial direction
     {
-        return this->first_order(x,y,f,xnew);
+        return this->first_order(x, y, f, xnew);
     }
-    if(ls_step == 1) //second point: third order extrapolation/interpolation
+    if (ls_step == 1) // second point: third order extrapolation/interpolation
     {
-        return this->third_order(x,y,f,xnew,conv_thr);
+        return this->third_order(x, y, f, xnew, conv_thr);
     }
 
-    //if still not converging after 2 steps: start Brent
-    if(ls_step == 2)
+    // if still not converging after 2 steps: start Brent
+    if (ls_step == 2)
     {
-        this->init_brent(x,y,f);
+        this->init_brent(x, y, f);
     }
-    if(ls_step >= 3)
+    if (ls_step >= 3)
     {
-        this->update_brent(x,y,f);
+        this->update_brent(x, y, f);
     }
 
-    if(ls_step >= 2)
+    if (ls_step >= 2)
     {
-        return this->brent(x,y,f,xnew,conv_thr);
+        return this->brent(x, y, f, xnew, conv_thr);
     }
-    ModuleBase::WARNING_QUIT("line_search","ls_step <0");
+    ModuleBase::WARNING_QUIT("line_search", "ls_step <0");
     __builtin_unreachable();
 }
 
-bool Line_Search::first_order(
-    const double x,
-    const double y,
-    const double f,
-    double & xnew)
+bool Line_Search::first_order(const double x, const double y, const double f, double& xnew)
 {
     xa = x; // set the first endpoint
     ya = y;
     fa = f;
     xnew = x + 1; // move one step towards trial direction
-    ls_step ++;
+    ls_step++;
     return false;
 }
 
-bool Line_Search::third_order(
-    const double x,
-    const double y,
-    const double f,
-    double & xnew,
-    const double conv_thr)
+bool Line_Search::third_order(const double x, const double y, const double f, double& xnew, const double conv_thr)
 {
-    double dmove, dmoveh, dmove1, dmove2;
+    double dmove = 0.0;
+    double dmoveh = 0.0;
+    double dmove1 = 0.0;
+    double dmove2 = 0.0;
 
     xb = x; // set the second endpoint
     yb = y;
     fb = f;
 
     double ffin = yb - ya;
-    double fab = (fa + fb)/2.0;
+    double fab = (fa + fb) / 2.0;
 
-    double k3 = 3.0*(fb    +fa-2.0*ffin);
-    double k2 =      fb+2.0*fa-3.0*ffin;
-    double k1 =      fa;
+    double k3 = 3.0 * (fb + fa - 2.0 * ffin);
+    double k2 = fb + 2.0 * fa - 3.0 * ffin;
+    double k1 = fa;
 
-    double tmp = k1*k3/k2/k2;
+    double tmp = k1 * k3 / k2 / k2;
 
-    //cubic extrapolation
-    if( std::abs(k3/k1)<1.0e-2 || std::abs(k3) < 2.4e-5 || tmp>1.0 ) //harmonic case
+    // cubic extrapolation
+    if (std::abs(k3 / k1) < 1.0e-2 || std::abs(k3) < 2.4e-5 || tmp > 1.0) // harmonic case
     {
-        dmove =  -fa/(fab-fa)/2.0;
-        if(dmove<0)
+        dmove = -fa / (fab - fa) / 2.0;
+        if (dmove < 0)
         {
-            dmove  = 4.0;
+            dmove = 4.0;
         }
     }
-    else //anharmonic case
+    else // anharmonic case
     {
-        dmove1 = k2/k3*(1.0-std::sqrt(1.0-tmp));
-        dmove2 = k2/k3*(1.0+std::sqrt(1.0-tmp));
+        dmove1 = k2 / k3 * (1.0 - std::sqrt(1.0 - tmp));
+        dmove2 = k2 / k3 * (1.0 + std::sqrt(1.0 - tmp));
 
-        double dy1, dy2;
-        dy1 = -(k1-(k2-k3*dmove1/3.0)*dmove1)*dmove1;
-        dy2 = -(k1-(k2-k3*dmove2/3.0)*dmove2)*dmove2;
+        double dy1 = -(k1 - (k2 - k3 * dmove1 / 3.0) * dmove1) * dmove1;
+        double dy2 = -(k1 - (k2 - k3 * dmove2 / 3.0) * dmove2) * dmove2;
 
-        if(dy1>dy2)
+        if (dy1 > dy2)
         {
             dmove = dmove1;
         }
@@ -104,35 +99,37 @@ bool Line_Search::third_order(
             dmove = dmove2;
         }
 
-        dmoveh = -fa/(fab-fa)/2.0;
-        if(dmoveh<0) dmoveh = 4.0;
+        dmoveh = -fa / (fab - fa) / 2.0;
+        if (dmoveh < 0) {
+            dmoveh = 4.0;
+}
 
-        if(dmove > 2.0*dmoveh || dmoveh > 2.0*dmove ||
-            (fa*fb > 0 && dmove < 1.0) ||
-            (fa*fb < 0 && dmove > 1.0) )
+        if (dmove > 2.0 * dmoveh || dmoveh > 2.0 * dmove || (fa * fb > 0 && dmove < 1.0)
+            || (fa * fb < 0 && dmove > 1.0))
         {
             dmove = dmoveh;
         }
-    } //end anharmonic case
+    } // end anharmonic case
 
-    if(dmove > 4.0) dmove = 4.0;
-    xnew  = dmove +xa;
+    if (dmove > 4.0) {
+        dmove = 4.0;
+}
+    xnew = dmove + xa;
 
-    double dy = (fb+(fab-fb)/(xa-xb)*(dmove-xb))*(dmove-xb);
-    if(std::abs(dy)<conv_thr) return true;
+    double dy = (fb + (fab - fb) / (xa - xb) * (dmove - xb)) * (dmove - xb);
+    if (std::abs(dy) < conv_thr) {
+        return true;
+}
 
-    ls_step ++;
+    ls_step++;
     return false;
 }
 
-void Line_Search::init_brent(
-    const double x,
-    const double y,
-    const double f)
+void Line_Search::init_brent(const double x, const double y, const double f)
 {
     bracked = true;
 
-    if(x>xb) // x > b, start interval [b,x]
+    if (x > xb) // x > b, start interval [b,x]
     {
         xa = xb;
         ya = yb;
@@ -140,18 +137,20 @@ void Line_Search::init_brent(
         xb = x;
         yb = y;
         fb = f;
-        if(fa*fb>0) bracked = false;
+        if (fa * fb > 0) {
+            bracked = false;
+}
         fstart = fa;
     }
     else // x < b
     {
-        if(fa*f<=0) // minimum between [a,x]
+        if (fa * f <= 0) // minimum between [a,x]
         {
             xb = x;
             yb = y;
             fb = f;
         }
-        else if (fb*f<=0) // minimum between [x,b]
+        else if (fb * f <= 0) // minimum between [x,b]
         {
             xa = xb;
             ya = yb;
@@ -160,7 +159,7 @@ void Line_Search::init_brent(
             yb = y;
             fb = f;
         }
-        else //problematic case, no minimum between [a,b]
+        else // problematic case, no minimum between [a,b]
         {
             xa = xb;
             ya = yb;
@@ -176,15 +175,12 @@ void Line_Search::init_brent(
     fc = fb;
 }
 
-void Line_Search::update_brent(
-    const double x,
-    const double y,
-    const double f)
+void Line_Search::update_brent(const double x, const double y, const double f)
 {
     xb = x;
     yb = y;
     fb = f;
-    if(!bracked && fstart*f<0)
+    if (!bracked && fstart * f < 0)
     {
         bracked = true;
         xc = xb;
@@ -192,84 +188,84 @@ void Line_Search::update_brent(
     }
 }
 
-bool Line_Search::brent(
-    const double x,
-    const double y,
-    const double f,
-    double & xnew,
-    const double conv_thr)
+bool Line_Search::brent(const double x, const double y, const double f, double& xnew, const double conv_thr)
 {
-    ls_step ++;
+    ls_step++;
 
-    double xd,xe,xm;
+    double xd = 0.0;
+    double xe = 0.0;
+    double xm = 0.0;
     // if no zero is between xa and xb
-    if(!bracked)
+    if (!bracked)
     {
-        if(std::abs(fc)<=std::abs(fb) || (xa-xc)<(xb-xa))
+        if (std::abs(fc) <= std::abs(fb) || (xa - xc) < (xb - xa))
         {
             xc = xa;
             fc = fa;
-            xd = xb-xa;
-            xe = xb-xa;
+            xd = xb - xa;
+            xe = xb - xa;
         }
-        double tol1 = 2.0*e8*std::abs(xb)+0.5*e8;
-        xm = 0.5*(xc-xb);
+        double tol1 = 2.0 * e8 * std::abs(xb) + 0.5 * e8;
+        xm = 0.5 * (xc - xb);
 
-        if(!(xc<=xa && xa<=xb))
+        if (!(xc <= xa && xa <= xb))
         {
-            ModuleBase::WARNING_QUIT("Brent","something wrong with Brent line search!");
+            ModuleBase::WARNING_QUIT("Brent", "something wrong with Brent line search!");
         }
-        if(std::abs(xm)<=tol1 || fb==0.0)
+        if (std::abs(xm) <= tol1 || fb == 0.0)
         {
             return true;
         }
-        if(std::abs(xe)>=tol1 && std::abs(fa)>std::abs(fb))
+        if (std::abs(xe) >= tol1 && std::abs(fa) > std::abs(fb))
         {
-            double s=fb/fa;
-            double p,qq;
-            if(xa==xc)
+            double s = fb / fa;
+            double p = 0.0;
+            double qq = 0.0;
+            if (xa == xc)
             {
-                p = 2.0*xm*s;
+                p = 2.0 * xm * s;
                 qq = 1.0 - s;
             }
             else
             {
-                qq = fa/fc;
-                double r = fb/fc;
-                p=s*(2.0*xm*qq*(qq-r)-(xb-xa)*(r-1.0));
-                qq=(qq-1.0)*(r-1.0)*(s-1.0);
+                qq = fa / fc;
+                double r = fb / fc;
+                p = s * (2.0 * xm * qq * (qq - r) - (xb - xa) * (r - 1.0));
+                qq = (qq - 1.0) * (r - 1.0) * (s - 1.0);
             }
-            if(p>0.0) qq=-qq;
+            if (p > 0.0) {
+                qq = -qq;
+}
             p = std::abs(p);
 
-            if( p < std::min(2.0*(xb-xa)*qq-std::abs(tol1*qq) , std::abs(xe*qq)/2.0) )
+            if (p < std::min(2.0 * (xb - xa) * qq - std::abs(tol1 * qq), std::abs(xe * qq) / 2.0))
             {
-                xe=xd;
-                xd=p/qq;
+                xe = xd;
+                xd = p / qq;
             }
             else
             {
-                xd=2.0*(xb-xa);
-                xe=xd;
+                xd = 2.0 * (xb - xa);
+                xe = xd;
             }
         }
         else
         {
-            xd = 2.0*(xb-xa);
+            xd = 2.0 * (xb - xa);
             xe = xd;
         }
 
-        double fab = (fa+fb)/2.0;
-        double dy = (fb+(fab-fb)/(xa-xb)*xd)*xd;
+        double fab = (fa + fb) / 2.0;
+        double dy = (fb + (fab - fb) / (xa - xb) * xd) * xd;
 
         xc = xa;
         fc = fa;
         xa = xb;
         fa = fb;
 
-        if(std::abs(xd)>tol1)
+        if (std::abs(xd) > tol1)
         {
-            xb = xb+xd;
+            xb = xb + xd;
         }
         else
         {
@@ -277,32 +273,35 @@ bool Line_Search::brent(
         }
 
         xnew = xb;
-        if(std::abs(dy)<conv_thr) return true;
-        if(ls_step == 4) //I'm not sure if this is a good choice, but the idea is there should not be so many line search steps
-                         //I feel if the line search does not converge, we'd better change the direction and restart line search
+        if (std::abs(dy) < conv_thr) {
+            return true;
+}
+        if (ls_step == 4) // I'm not sure if this is a good choice, but the idea is there should not be so many line
+                          // search steps I feel if the line search does not converge, we'd better change the direction
+                          // and restart line search
         {
             GlobalV::ofs_running << "Too many Brent steps, let's do next CG step" << std::endl;
             return true;
-            //ModuleBase::WARNING_QUIT("Brent","too many steps in line search, something wrong");
+            // ModuleBase::WARNING_QUIT("Brent","too many steps in line search, something wrong");
         }
 
         return false;
-    }//end bracked
+    } // end bracked
     else
     {
-        if(!( (xa<=xb && xb<=xc) || (xc<=xb && xb<=xa) ))
+        if (!((xa <= xb && xb <= xc) || (xc <= xb && xb <= xa)))
         {
-            ModuleBase::WARNING_QUIT("Brent","something wrong with Brent line search!");
+            ModuleBase::WARNING_QUIT("Brent", "something wrong with Brent line search!");
         }
 
-        if((fb>0 && fc>0) || (fb<0 && fc<0))
+        if ((fb > 0 && fc > 0) || (fb < 0 && fc < 0))
         {
             xc = xa;
             fc = fa;
-            xd = xb-xa;
+            xd = xb - xa;
             xe = xd;
         }
-        if(std::abs(fc)<std::abs(fb)) //rearrange b,c so that abs(fc)>=abs(fb)
+        if (std::abs(fc) < std::abs(fb)) // rearrange b,c so that abs(fc)>=abs(fb)
         {
             xa = xb;
             xb = xc;
@@ -312,37 +311,39 @@ bool Line_Search::brent(
             fc = fa;
         }
 
-        double tol1 = 2.0*e8*std::abs(xb)+0.5*e8;
-        xm = 0.5*(xc-xb);
-        if(std::abs(xm)<=tol1 || fb==0.0)
+        double tol1 = 2.0 * e8 * std::abs(xb) + 0.5 * e8;
+        xm = 0.5 * (xc - xb);
+        if (std::abs(xm) <= tol1 || fb == 0.0)
         {
             return true;
         }
 
-        if(std::abs(xe)>=tol1 && std::abs(fa)>std::abs(fb))
+        if (std::abs(xe) >= tol1 && std::abs(fa) > std::abs(fb))
         {
-            double s=fb/fa;
-            double p,
-            qq;
-            if(xa==xc)
+            double s = fb / fa;
+            double p = 0.0;
+            double qq = 0.0;
+            if (xa == xc)
             {
-                p = 2.0*xm*s;
+                p = 2.0 * xm * s;
                 qq = 1.0 - s;
             }
             else
             {
-                qq = fa/fc;
-                double r = fb/fc;
-                p=s*(2.0*xm*qq*(qq-r)-(xb-xa)*(r-1.0));
-                qq=(qq-1.0)*(r-1.0)*(s-1.0);
+                qq = fa / fc;
+                double r = fb / fc;
+                p = s * (2.0 * xm * qq * (qq - r) - (xb - xa) * (r - 1.0));
+                qq = (qq - 1.0) * (r - 1.0) * (s - 1.0);
             }
-            if(p>0.0) qq=-qq;
+            if (p > 0.0) {
+                qq = -qq;
+}
             p = std::abs(p);
 
-            if(2.0*p < std::min( 3.0*xm*qq-std::abs(tol1*qq), std::abs(xe*qq) ) )
+            if (2.0 * p < std::min(3.0 * xm * qq - std::abs(tol1 * qq), std::abs(xe * qq)))
             {
                 xe = xd;
-                xd = p/qq;
+                xd = p / qq;
             }
             else
             {
@@ -356,20 +357,20 @@ bool Line_Search::brent(
             xe = xd;
         }
 
-        double fab = (fa+fb)/2.0;
-        double dy = (fb+(fab-fb)/(xa-xb)*xd)*xd;
+        double fab = (fa + fb) / 2.0;
+        double dy = (fb + (fab - fb) / (xa - xb) * xd) * xd;
 
         xa = xb;
         fa = fb;
         ya = yb;
 
-        if(std::abs(xc)>tol1)
+        if (std::abs(xc) > tol1)
         {
             xb = xb + xd;
         }
         else
         {
-            if(xm > 0)
+            if (xm > 0)
             {
                 xb = xb + tol1;
             }
@@ -380,13 +381,15 @@ bool Line_Search::brent(
         }
 
         xnew = xb;
-        if(std::abs(dy)<conv_thr) return true;
-        if(ls_step == 4)
+        if (std::abs(dy) < conv_thr) {
+            return true;
+}
+        if (ls_step == 4)
         {
             GlobalV::ofs_running << "Too many Brent steps, let's do next CG step" << std::endl;
             return true;
-            //ModuleBase::WARNING_QUIT("Brent","too many steps in line search, something wrong");
+            // ModuleBase::WARNING_QUIT("Brent","too many steps in line search, something wrong");
         }
         return false;
-    }//end ibrack
+    } // end ibrack
 }
\ No newline at end of file
diff --git a/source/module_relax/relax_new/line_search.h b/source/module_relax/relax_new/line_search.h
index 3cfa27197f..8277c31923 100644
--- a/source/module_relax/relax_new/line_search.h
+++ b/source/module_relax/relax_new/line_search.h
@@ -3,55 +3,41 @@
 
 class Line_Search
 {
-    public:
-
-    Line_Search(){};
-    ~Line_Search(){};
-
-    bool line_search(
-        const bool restart, //whether to restart line search (when cg direction has changed)
-        const double x, //position for function at x
-        const double y, //value of function at x
-        const double f, //gradient of function at x
-        double & xnew, //postion where function has to be evaluated
-        const double conv_thr);  //predicted change of function value of function at xnew
-
-    bool first_order(
-        const double x,
-        const double y,
-        const double f,
-        double & xnew);
-
-    bool third_order(
-        const double x,
-        const double y,
-        const double f,
-        double & xnew,
-        const double conv_thr);
-
-    void init_brent(
-        const double x,
-        const double y,
-        const double f);
-
-    void update_brent(
-        const double x,
-        const double y,
-        const double f);
-
-    bool brent(
-        const double x,
-        const double y,
-        const double f,
-        double & xnew,
-        const double conv_thr);
-
-    private:
-
-    int ls_step = 0; // step of line search
-    bool bracked; //whether the minima is bracked by [a,c]
-    double xa,xb,xc,fa,fb,fc,ya,yb; //keep record of some points
-    double fstart; //keep record of initial gradient for brent method
+  public:
+    Line_Search() {};
+    ~Line_Search() {};
+
+    bool line_search(const bool restart,     // whether to restart line search (when cg direction has changed)
+                     const double x,         // position for function at x
+                     const double y,         // value of function at x
+                     const double f,         // gradient of function at x
+                     double& xnew,           // postion where function has to be evaluated
+                     const double conv_thr); // predicted change of function value of function at xnew
+
+    bool first_order(const double x, const double y, const double f, double& xnew);
+
+    bool third_order(const double x, const double y, const double f, double& xnew, const double conv_thr);
+
+    void init_brent(const double x, const double y, const double f);
+
+    void update_brent(const double x, const double y, const double f);
+
+    bool brent(const double x, const double y, const double f, double& xnew, const double conv_thr);
+
+  private:
+    int ls_step = 0;      // step of line search
+    bool bracked = false; // whether the minima is bracked by [a,c]
+    /// these variables used to keep record of some points
+    double xa = 0.0;
+    double xb = 0.0;
+    double xc = 0.0;
+    double fa = 0.0;
+    double fb = 0.0;
+    double fc = 0.0;
+    double ya = 0.0;
+    double yb = 0.0;
+    /// ------------------------------------
+    double fstart = 0.0; // keep record of initial gradient for brent method
     const double e8 = 1.0e-8;
 };
 
diff --git a/source/module_relax/relax_new/relax.cpp b/source/module_relax/relax_new/relax.cpp
index 19c7be712d..e873a7aa50 100644
--- a/source/module_relax/relax_new/relax.cpp
+++ b/source/module_relax/relax_new/relax.cpp
@@ -1,19 +1,19 @@
 #include "relax.h"
 
-#include "module_parameter/parameter.h"
-#include <cmath>
-
-#include "module_relax/relax_old/ions_move_basic.h"
 #include "module_base/matrix3.h"
 #include "module_base/parallel_common.h"
 #include "module_base/tool_title.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
+#include "module_parameter/parameter.h"
+#include "module_relax/relax_old/ions_move_basic.h"
+
+#include <cmath>
 
 void Relax::init_relax(const int nat_in)
 {
-    ModuleBase::TITLE("Relax","init_relax");
+    ModuleBase::TITLE("Relax", "init_relax");
 
-    //set some initial conditions / constants
+    // set some initial conditions / constants
     nat = nat_in;
     istep = 0;
     cg_step = 0;
@@ -24,48 +24,50 @@ void Relax::init_relax(const int nat_in)
     etot = 0;
     etot_p = 0;
 
-    force_thr_eva = PARAM.inp.force_thr  * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A; //convert to eV/A
-    fac_force  = PARAM.inp.relax_scale_force * 0.1;
+    force_thr_eva = PARAM.inp.force_thr * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A; // convert to eV/A
+    fac_force = PARAM.inp.relax_scale_force * 0.1;
     fac_stress = fac_force / nat;
 
-    //allocate some data structures
+    // allocate some data structures
 
-    //gradients in ion and cell; current and previous
-    grad_ion.create(nat,3);
-    grad_cell.create(3,3);
+    // gradients in ion and cell; current and previous
+    grad_ion.create(nat, 3);
+    grad_cell.create(3, 3);
 
-    grad_ion_p.create(nat,3);
-    grad_cell_p.create(3,3);
+    grad_ion_p.create(nat, 3);
+    grad_cell_p.create(3, 3);
 
-    //search directions for ion and cell; current and previous
-    search_dr_ion.create(nat,3);
-    search_dr_cell.create(3,3);
+    // search directions for ion and cell; current and previous
+    search_dr_ion.create(nat, 3);
+    search_dr_cell.create(3, 3);
 
-    search_dr_ion_p.create(nat,3);
-    search_dr_cell_p.create(3,3);
+    search_dr_ion_p.create(nat, 3);
+    search_dr_cell_p.create(3, 3);
 
-    //set if we are allowing lattice vectors to move
+    // set if we are allowing lattice vectors to move
     if_cell_moves = false;
-    if(PARAM.inp.calculation == "cell-relax") { if_cell_moves = true;
-}
+    if (PARAM.inp.calculation == "cell-relax")
+    {
+        if_cell_moves = true;
+    }
 }
 
 bool Relax::relax_step(UnitCell& ucell,
-                       const ModuleBase::matrix& force, 
-                       const ModuleBase::matrix &stress, 
+                       const ModuleBase::matrix& force,
+                       const ModuleBase::matrix& stress,
                        const double etot_in)
 {
-    ModuleBase::TITLE("Relax","relax_step");
+    ModuleBase::TITLE("Relax", "relax_step");
 
-    etot = etot_in * ModuleBase::Ry_to_eV; //convert to eV
-    if( istep == 0 ) 
-    { 
+    etot = etot_in * ModuleBase::Ry_to_eV; // convert to eV
+    if (istep == 0)
+    {
         etot_p = etot;
     }
 
-    bool relax_done = this->setup_gradient(ucell,force, stress);
-    if(relax_done) 
-    { 
+    bool relax_done = this->setup_gradient(ucell, force, stress);
+    if (relax_done)
+    {
         return relax_done;
     }
 
@@ -73,154 +75,170 @@ bool Relax::relax_step(UnitCell& ucell,
 
     bool ls_done = this->check_line_search();
 
-    if(ls_done)
+    if (ls_done)
     {
         this->new_direction();
-        this->move_cell_ions(ucell,true);
+        this->move_cell_ions(ucell, true);
     }
     else
     {
         this->perform_line_search();
-        this->move_cell_ions(ucell,false);
+        this->move_cell_ions(ucell, false);
         dmovel = dmove;
     }
 
-    istep ++;
+    istep++;
 
     return relax_done;
 }
 
-bool Relax::setup_gradient(const UnitCell &ucell,
-                           const ModuleBase::matrix& force, 
-                           const ModuleBase::matrix &stress)
+bool Relax::setup_gradient(const UnitCell& ucell, const ModuleBase::matrix& force, const ModuleBase::matrix& stress)
 {
-    ModuleBase::TITLE("Relax","setup_gradient");
+    ModuleBase::TITLE("Relax", "setup_gradient");
 
-    //if not relax, then return converged
-    if( !( PARAM.inp.calculation == "relax" || PARAM.inp.calculation == "cell-relax" ) ) 
-    { 
+    // if not relax, then return converged
+    if (!(PARAM.inp.calculation == "relax" || PARAM.inp.calculation == "cell-relax"))
+    {
         return true;
     }
 
-    //indicating whether force & stress are converged
+    // indicating whether force & stress are converged
     bool force_converged = true;
     double max_grad = 0.0;
 
     //=========================================
-    //set gradient for ions degrees of freedom
+    // set gradient for ions degrees of freedom
     //=========================================
 
     grad_ion.zero_out();
-    ModuleBase::matrix force_eva = force * (ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A); //convert to eV/A
-
-	int iat=0;
-	for(int it = 0;it < ucell.ntype;it++)
-	{
-		Atom* atom = &ucell.atoms[it];
-		for(int ia =0;ia< ucell.atoms[it].na;ia++)
-		{	
+    ModuleBase::matrix force_eva = force * (ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A); // convert to eV/A
+
+    int iat = 0;
+    for (int it = 0; it < ucell.ntype; it++)
+    {
+        Atom* atom = &ucell.atoms[it];
+        for (int ia = 0; ia < ucell.atoms[it].na; ia++)
+        {
             double force2 = 0.0;
-			if(atom->mbl[ia].x == 1)
-			{
-				grad_ion(iat, 0) = force_eva(iat, 0);
-                if( std::abs(force_eva(iat,0)) > max_grad) 
-                { 
-                    max_grad = std::abs(force_eva(iat,0));
+            if (atom->mbl[ia].x == 1)
+            {
+                grad_ion(iat, 0) = force_eva(iat, 0);
+                if (std::abs(force_eva(iat, 0)) > max_grad)
+                {
+                    max_grad = std::abs(force_eva(iat, 0));
                 }
-			}
-			if(atom->mbl[ia].y == 1)
-			{
-				grad_ion(iat, 1) = force_eva(iat, 1);
-                if( std::abs(force_eva(iat,1)) > max_grad) 
-                { 
-                    max_grad = std::abs(force_eva(iat,1));
+            }
+            if (atom->mbl[ia].y == 1)
+            {
+                grad_ion(iat, 1) = force_eva(iat, 1);
+                if (std::abs(force_eva(iat, 1)) > max_grad)
+                {
+                    max_grad = std::abs(force_eva(iat, 1));
                 }
-			}
-			if(atom->mbl[ia].z == 1)
-			{
-				grad_ion(iat, 2) = force_eva(iat, 2);
-                if( std::abs(force_eva(iat,2)) > max_grad) 
-                { 
-                    max_grad = std::abs(force_eva(iat,2));
+            }
+            if (atom->mbl[ia].z == 1)
+            {
+                grad_ion(iat, 2) = force_eva(iat, 2);
+                if (std::abs(force_eva(iat, 2)) > max_grad)
+                {
+                    max_grad = std::abs(force_eva(iat, 2));
                 }
-			}
-			++iat;
-		}
-	}
-    assert(iat==nat);
-
-    if(max_grad > force_thr_eva) 
-    { 
+            }
+            ++iat;
+        }
+    }
+    assert(iat == nat);
+
+    if (max_grad > force_thr_eva)
+    {
         force_converged = false;
     }
-	if(PARAM.inp.out_level=="ie")
-	{
+    if (PARAM.inp.out_level == "ie")
+    {
         std::cout << " ETOT DIFF (eV)       : " << etot - etot_p << std::endl;
-		std::cout << " LARGEST GRAD (eV/A)  : " << max_grad << std::endl;
+        std::cout << " LARGEST GRAD (eV/A)  : " << max_grad << std::endl;
         etot_p = etot;
-	}
-    
+    }
+
     GlobalV::ofs_running << "\n Largest gradient in force is " << max_grad << " eV/A." << std::endl;
     GlobalV::ofs_running << " Threshold is " << PARAM.inp.force_thr_ev << " eV/A." << std::endl;
-//=========================================
-//set gradient for cell degrees of freedom
-//=========================================
+    //=========================================
+    // set gradient for cell degrees of freedom
+    //=========================================
 
-    if(if_cell_moves)
+    if (if_cell_moves)
     {
         grad_cell.zero_out();
         ModuleBase::matrix stress_ev = stress * (ucell.omega * ModuleBase::Ry_to_eV);
 
-        if(PARAM.inp.fixed_axes == "shape")
+        if (PARAM.inp.fixed_axes == "shape")
         {
-            double pressure = (stress_ev(0,0) + stress_ev(1,1) + stress_ev(2,2)) / 3.0;
+            double pressure = (stress_ev(0, 0) + stress_ev(1, 1) + stress_ev(2, 2)) / 3.0;
             stress_ev.zero_out();
-            stress_ev(0,0) = pressure; // apply constraints
-            stress_ev(1,1) = pressure;
-            stress_ev(2,2) = pressure;
+            stress_ev(0, 0) = pressure; // apply constraints
+            stress_ev(1, 1) = pressure;
+            stress_ev(2, 2) = pressure;
         }
-        else if(PARAM.inp.fixed_axes == "volume")
+        else if (PARAM.inp.fixed_axes == "volume")
         {
-            double pressure = (stress_ev(0,0) + stress_ev(1,1) + stress_ev(2,2)) / 3.0;
-            stress_ev(0,0) -= pressure;
-            stress_ev(1,1) -= pressure;
-            stress_ev(2,2) -= pressure;
+            double pressure = (stress_ev(0, 0) + stress_ev(1, 1) + stress_ev(2, 2)) / 3.0;
+            stress_ev(0, 0) -= pressure;
+            stress_ev(1, 1) -= pressure;
+            stress_ev(2, 2) -= pressure;
         }
         else if (PARAM.inp.fixed_axes != "None")
         {
-            //Note stress is given in the directions of lattice vectors
-            //So we need to first convert to Cartesian and then apply the constraint
+            // Note stress is given in the directions of lattice vectors
+            // So we need to first convert to Cartesian and then apply the constraint
             ModuleBase::Matrix3 stress_cart;
-            stress_cart.e11 = stress_ev(0,0); stress_cart.e12 = stress_ev(0,1); stress_cart.e13 = stress_ev(0,2);
-            stress_cart.e21 = stress_ev(1,0); stress_cart.e22 = stress_ev(1,1); stress_cart.e23 = stress_ev(1,2);
-            stress_cart.e31 = stress_ev(2,0); stress_cart.e32 = stress_ev(2,1); stress_cart.e33 = stress_ev(2,2);
+            stress_cart.e11 = stress_ev(0, 0);
+            stress_cart.e12 = stress_ev(0, 1);
+            stress_cart.e13 = stress_ev(0, 2);
+            stress_cart.e21 = stress_ev(1, 0);
+            stress_cart.e22 = stress_ev(1, 1);
+            stress_cart.e23 = stress_ev(1, 2);
+            stress_cart.e31 = stress_ev(2, 0);
+            stress_cart.e32 = stress_ev(2, 1);
+            stress_cart.e33 = stress_ev(2, 2);
 
             stress_cart = ucell.latvec * stress_cart;
 
-            if(ucell.lc[0] == 0)
+            if (ucell.lc[0] == 0)
             {
-                stress_cart.e11 = 0; stress_cart.e12 = 0; stress_cart.e13 = 0;
+                stress_cart.e11 = 0;
+                stress_cart.e12 = 0;
+                stress_cart.e13 = 0;
             }
-            if(ucell.lc[1] == 0)
+            if (ucell.lc[1] == 0)
             {
-                stress_cart.e21 = 0; stress_cart.e22 = 0; stress_cart.e23 = 0;
+                stress_cart.e21 = 0;
+                stress_cart.e22 = 0;
+                stress_cart.e23 = 0;
             }
-            if(ucell.lc[2] == 0)
+            if (ucell.lc[2] == 0)
             {
-                stress_cart.e31 = 0; stress_cart.e32 = 0; stress_cart.e33 = 0;
+                stress_cart.e31 = 0;
+                stress_cart.e32 = 0;
+                stress_cart.e33 = 0;
             }
 
             stress_cart = ucell.GT * stress_cart;
-            stress_ev(0,0) = stress_cart.e11; stress_ev(0,1) = stress_cart.e12; stress_ev(0,2) = stress_cart.e13;
-            stress_ev(1,0) = stress_cart.e21; stress_ev(1,1) = stress_cart.e22; stress_ev(1,2) = stress_cart.e23;
-            stress_ev(2,0) = stress_cart.e31; stress_ev(2,1) = stress_cart.e32; stress_ev(2,2) = stress_cart.e33;           
+            stress_ev(0, 0) = stress_cart.e11;
+            stress_ev(0, 1) = stress_cart.e12;
+            stress_ev(0, 2) = stress_cart.e13;
+            stress_ev(1, 0) = stress_cart.e21;
+            stress_ev(1, 1) = stress_cart.e22;
+            stress_ev(1, 2) = stress_cart.e23;
+            stress_ev(2, 0) = stress_cart.e31;
+            stress_ev(2, 1) = stress_cart.e32;
+            stress_ev(2, 2) = stress_cart.e33;
         }
 
-        for(int i=0;i<3;i++)
+        for (int i = 0; i < 3; i++)
         {
-            for(int j=0;j<3;j++)
+            for (int j = 0; j < 3; j++)
             {
-                grad_cell(i,j) = stress_ev(i,j); // apply constraints
+                grad_cell(i, j) = stress_ev(i, j); // apply constraints
             }
         }
 
@@ -231,9 +249,9 @@ bool Relax::setup_gradient(const UnitCell &ucell,
         {
             for (int j = 0; j < 3; j++)
             {
-                double grad = grad_cell(i,j) / (ucell.omega * ModuleBase::Ry_to_eV);
-                if ( largest_grad < std::abs(grad) ) 
-                { 
+                double grad = grad_cell(i, j) / (ucell.omega * ModuleBase::Ry_to_eV);
+                if (largest_grad < std::abs(grad))
+                {
                     largest_grad = std::abs(grad);
                 }
             }
@@ -247,227 +265,229 @@ bool Relax::setup_gradient(const UnitCell &ucell,
             force_converged = false;
         }
 
-        GlobalV::ofs_running << "\n Largest gradient in stress is " << largest_grad  << " kbar." << std::endl;
+        GlobalV::ofs_running << "\n Largest gradient in stress is " << largest_grad << " kbar." << std::endl;
         GlobalV::ofs_running << " Threshold is " << PARAM.inp.stress_thr << " kbar." << std::endl;
     }
 
-	if(force_converged)
-	{
-		GlobalV::ofs_running << "\n Relaxation is converged!" << std::endl;
-	}
-	else
-	{
-		GlobalV::ofs_running << "\n Relaxation is not converged yet!" << std::endl;
-	}
+    if (force_converged)
+    {
+        GlobalV::ofs_running << "\n Relaxation is converged!" << std::endl;
+    }
+    else
+    {
+        GlobalV::ofs_running << "\n Relaxation is not converged yet!" << std::endl;
+    }
 
     return force_converged;
 }
 
 void Relax::calculate_gamma()
 {
-    ModuleBase::TITLE("Relax","calculate_gamma");
+    ModuleBase::TITLE("Relax", "calculate_gamma");
 
-    //no need to calculate gamma if last step is trial
-    //since we won't update search direction
-    if(ltrial)
+    // no need to calculate gamma if last step is trial
+    // since we won't update search direction
+    if (ltrial)
     {
         return;
     }
 
-    grp_grp = 0.0; //grad_p*grad_p
-    gr_grp  = 0.0; //grad  *grad_p
-    gr_gr = 0.0;   //grad  *grad
-    gr_sr = 0.0;   //grad  *search_dir
+    grp_grp = 0.0; // grad_p*grad_p
+    gr_grp = 0.0;  // grad  *grad_p
+    gr_gr = 0.0;   // grad  *grad
+    gr_sr = 0.0;   // grad  *search_dir
 
-    for(int iat=0; iat<nat; iat++)
+    for (int iat = 0; iat < nat; iat++)
     {
-        for(int i=0;i<3;i++)
+        for (int i = 0; i < 3; i++)
         {
-            grp_grp += grad_ion_p(iat,i) *    grad_ion_p(iat,i);
-            gr_grp  += grad_ion_p(iat,i) *      grad_ion(iat,i);
-            gr_gr   += grad_ion(iat,i) *      grad_ion(iat,i);
-            gr_sr   += grad_ion(iat,i) * search_dr_ion(iat,i);
+            grp_grp += grad_ion_p(iat, i) * grad_ion_p(iat, i);
+            gr_grp += grad_ion_p(iat, i) * grad_ion(iat, i);
+            gr_gr += grad_ion(iat, i) * grad_ion(iat, i);
+            gr_sr += grad_ion(iat, i) * search_dr_ion(iat, i);
         }
     }
 
-    if(if_cell_moves)
+    if (if_cell_moves)
     {
-        for(int i=0;i<3;i++)
+        for (int i = 0; i < 3; i++)
         {
-            for(int j=0;j<3;j++)
+            for (int j = 0; j < 3; j++)
             {
-                grp_grp += grad_cell_p(i,j) *  grad_cell_p(i,j) / nat;
-                gr_grp  += grad_cell_p(i,j) *    grad_cell(i,j) / nat;
-                gr_gr   += grad_cell(i,j) *      grad_cell(i,j) / nat;
-                gr_sr   += grad_cell(i,j) * search_dr_cell(i,j) / nat;
+                grp_grp += grad_cell_p(i, j) * grad_cell_p(i, j) / nat;
+                gr_grp += grad_cell_p(i, j) * grad_cell(i, j) / nat;
+                gr_gr += grad_cell(i, j) * grad_cell(i, j) / nat;
+                gr_sr += grad_cell(i, j) * search_dr_cell(i, j) / nat;
             }
         }
     }
 
-    if(cg_step == 0)
+    if (cg_step == 0)
     {
         gamma = 0.0;
     }
     else
-    {   
-        gamma = (gr_gr - gr_grp) / grp_grp; //Polak-Riebere
+    {
+        gamma = (gr_gr - gr_grp) / grp_grp; // Polak-Riebere
     }
 }
 
 bool Relax::check_line_search()
 {
-    ModuleBase::TITLE("Relax","check_line_search");
-    
-    //if last step is trial step towards new direction
-    //then line search is not finished
-    //we will perform line search
-    if(ltrial)
+    ModuleBase::TITLE("Relax", "check_line_search");
+
+    // if last step is trial step towards new direction
+    // then line search is not finished
+    // we will perform line search
+    if (ltrial)
     {
         ltrial = false;
         return false;
     }
 
-    if(std::abs(gr_sr)*std::max(gamma,1.0) > std::abs(gr_gr)/5.0 && !brent_done) //last brent line search not finished
+    if (std::abs(gr_sr) * std::max(gamma, 1.0) > std::abs(gr_gr) / 5.0
+        && !brent_done) // last brent line search not finished
     {
         return false;
     }
-    
+
     return true;
 }
 
 void Relax::perform_line_search()
 {
-    ModuleBase::TITLE("Relax","line_search");
+    ModuleBase::TITLE("Relax", "line_search");
 
     double f = 0.0; // 1st order energy difference
-    for(int iat=0;iat<nat;iat++)
+    for (int iat = 0; iat < nat; iat++)
     {
-        for(int i=0;i<3;i++)
+        for (int i = 0; i < 3; i++)
         {
-            f -= step_size * fac_force * search_dr_ion(iat,i) *grad_ion(iat,i);
+            f -= step_size * fac_force * search_dr_ion(iat, i) * grad_ion(iat, i);
         }
     }
 
-    if(if_cell_moves)
+    if (if_cell_moves)
     {
-        for(int i=0;i<3;i++)
+        for (int i = 0; i < 3; i++)
         {
-            for(int j=0;j<3;j++)
+            for (int j = 0; j < 3; j++)
             {
-                f -= step_size * fac_stress * search_dr_cell(i,j) * grad_cell(i,j);
+                f -= step_size * fac_stress * search_dr_cell(i, j) * grad_cell(i, j);
             }
         }
     }
 
-    //perform line search
+    // perform line search
     bool restart_brent = false;
-    double x=dmovel, y=etot;
+    double x = dmovel, y = etot;
     double xnew, yd;
 
     brent_done = this->ls.line_search(restart_brent, x, y, f, xnew, force_thr_eva);
-    dmove  = xnew;
+    dmove = xnew;
 
     return;
 }
 
 void Relax::new_direction()
 {
-    ModuleBase::TITLE("Relax","new_direction");
-    if(cg_step != 0) { step_size += 0.2 * step_size * (dmovel - 1.0);
-}
+    ModuleBase::TITLE("Relax", "new_direction");
+    if (cg_step != 0)
+    {
+        step_size += 0.2 * step_size * (dmovel - 1.0);
+    }
 
-    //set GAMMA to zero if line minimization was not sufficient
-    if(5.0*std::abs(gr_sr)*gamma>std::abs(gr_gr))
+    // set GAMMA to zero if line minimization was not sufficient
+    if (5.0 * std::abs(gr_sr) * gamma > std::abs(gr_gr))
     {
         gamma = 0.0;
-        cg_step = 0; //reset cg
+        cg_step = 0; // reset cg
     }
 
-    //perform trial step
+    // perform trial step
 
-    //set search vectors
-    search_dr_ion_p  = search_dr_ion;
+    // set search vectors
+    search_dr_ion_p = search_dr_ion;
     search_dr_cell_p = search_dr_cell;
 
-    grad_ion_p  = grad_ion;
+    grad_ion_p = grad_ion;
     grad_cell_p = grad_cell;
 
     search_dr_ion = grad_ion + gamma * search_dr_ion;
     search_dr_cell = grad_cell + gamma * search_dr_cell;
 
-    //modify step if necessary
+    // modify step if necessary
     sr_sr = 1.0e-10;
-    for(int iat=0;iat<nat;iat++)
+    for (int iat = 0; iat < nat; iat++)
     {
-        for(int i=0;i<3;i++)
+        for (int i = 0; i < 3; i++)
         {
-            sr_sr += 1/fac_force * search_dr_ion(iat,i) * search_dr_ion(iat,i);
+            sr_sr += 1 / fac_force * search_dr_ion(iat, i) * search_dr_ion(iat, i);
         }
     }
 
-    if(if_cell_moves)
+    if (if_cell_moves)
     {
-        for(int i=0;i<3;i++)
+        for (int i = 0; i < 3; i++)
         {
-            for(int j=0;j<3;j++)
+            for (int j = 0; j < 3; j++)
             {
-                sr_sr += 1/fac_stress * search_dr_cell(i,j) * search_dr_cell(i,j);
+                sr_sr += 1 / fac_stress * search_dr_cell(i, j) * search_dr_cell(i, j);
             }
         }
     }
 
-    //if length of search vector increased, rescale step to avoid too large trial steps
-    if(sr_sr > srp_srp)
+    // if length of search vector increased, rescale step to avoid too large trial steps
+    if (sr_sr > srp_srp)
     {
         step_size *= srp_srp / sr_sr;
     }
     srp_srp = sr_sr;
-    
-    double f = 0.0; //first order change in energy (gradient in the search direction)
-    for(int iat=0;iat<nat;iat++)
+
+    double f = 0.0; // first order change in energy (gradient in the search direction)
+    for (int iat = 0; iat < nat; iat++)
     {
-        for(int i=0;i<3;i++)
+        for (int i = 0; i < 3; i++)
         {
-            f -= step_size * fac_force * search_dr_ion(iat,i) * grad_ion(iat,i);
+            f -= step_size * fac_force * search_dr_ion(iat, i) * grad_ion(iat, i);
         }
     }
 
-    if(if_cell_moves)
+    if (if_cell_moves)
     {
-        for(int i=0;i<3;i++)
+        for (int i = 0; i < 3; i++)
         {
-            for(int j=0;j<3;j++)
+            for (int j = 0; j < 3; j++)
             {
-                f -= step_size * fac_stress * search_dr_cell(i,j) * grad_cell(i,j);
+                f -= step_size * fac_stress * search_dr_cell(i, j) * grad_cell(i, j);
             }
         }
     }
 
-    //prepare for line search
+    // prepare for line search
     bool restart = true;
-    double x=0, y=etot;
-    //TODO: add a certain threshold for the progress.
+    double x = 0, y = etot;
+    // TODO: add a certain threshold for the progress.
     double xnew, yd = 1e-8;
 
     this->ls.line_search(restart, x, y, f, xnew, yd);
 
     dmovel = 1.0;
     ltrial = true;
-    cg_step ++;
+    cg_step++;
 
     return;
 }
 
-void Relax::move_cell_ions(UnitCell& ucell,
-                           const bool is_new_dir)
+void Relax::move_cell_ions(UnitCell& ucell, const bool is_new_dir)
 {
-    ModuleBase::TITLE("Relax","move_cell_ions");
+    ModuleBase::TITLE("Relax", "move_cell_ions");
 
     // I'm keeping this only because we have to
     // be compatible with old code
     ucell.ionic_position_updated = true;
-    if(if_cell_moves) 
-    { 
+    if (if_cell_moves)
+    {
         ucell.cell_parameter_updated = true;
     }
 
@@ -476,7 +496,7 @@ void Relax::move_cell_ions(UnitCell& ucell,
     // and the input variable is_new_dir is used to make the distinction
 
     double fac; // fac1 for force, fac2 for stress
-    if(is_new_dir)
+    if (is_new_dir)
     {
         fac = 1.0;
     }
@@ -493,27 +513,26 @@ void Relax::move_cell_ions(UnitCell& ucell,
     // Thirdly, in update_pos_taud, update Cartesian coordinates of atoms
     // in this step we are using the NEW latvec (lattice vectors)
     // Finally, update G, GT and other stuff, and print the new STRU, and update something for next SCF
-    
+
     // =================================================================
     // Step 1 : updating latvec
     // =================================================================
-    
-    if(if_cell_moves)
+
+    if (if_cell_moves)
     {
         // imo matrix3 class is not a very clever way to store 3*3 matrix ...
         ModuleBase::Matrix3 sr_dr_cell;
-        auto cp_mat_to_mat3 = [&sr_dr_cell, this]() -> void
-            {
-                sr_dr_cell.e11 = search_dr_cell(0, 0);
-                sr_dr_cell.e12 = search_dr_cell(0, 1);
-                sr_dr_cell.e13 = search_dr_cell(0, 2);
-                sr_dr_cell.e21 = search_dr_cell(1, 0);
-                sr_dr_cell.e22 = search_dr_cell(1, 1);
-                sr_dr_cell.e23 = search_dr_cell(1, 2);
-                sr_dr_cell.e31 = search_dr_cell(2, 0);
-                sr_dr_cell.e32 = search_dr_cell(2, 1);
-                sr_dr_cell.e33 = search_dr_cell(2, 2);
-            };
+        auto cp_mat_to_mat3 = [&sr_dr_cell, this]() -> void {
+            sr_dr_cell.e11 = search_dr_cell(0, 0);
+            sr_dr_cell.e12 = search_dr_cell(0, 1);
+            sr_dr_cell.e13 = search_dr_cell(0, 2);
+            sr_dr_cell.e21 = search_dr_cell(1, 0);
+            sr_dr_cell.e22 = search_dr_cell(1, 1);
+            sr_dr_cell.e23 = search_dr_cell(1, 2);
+            sr_dr_cell.e31 = search_dr_cell(2, 0);
+            sr_dr_cell.e32 = search_dr_cell(2, 1);
+            sr_dr_cell.e33 = search_dr_cell(2, 2);
+        };
         cp_mat_to_mat3();
 
         if (ModuleSymmetry::Symmetry::symm_flag && ucell.symm.nrotk > 0)
@@ -529,34 +548,40 @@ void Relax::move_cell_ions(UnitCell& ucell,
         // different from when the current CG step starts;
         // as a result, we need to save latvec at the beginning of
         // each CG step
-        if(is_new_dir) 
-        { 
+        if (is_new_dir)
+        {
             latvec_save = ucell.latvec;
         }
 
         ModuleBase::Matrix3 move_cell = latvec_save * sr_dr_cell;
 
-        //should be close to 0, but set again to avoid numerical issues
-        if(ucell.lc[0] == 0)
+        // should be close to 0, but set again to avoid numerical issues
+        if (ucell.lc[0] == 0)
         {
-            move_cell.e11 = 0; move_cell.e12 = 0; move_cell.e13 = 0;
+            move_cell.e11 = 0;
+            move_cell.e12 = 0;
+            move_cell.e13 = 0;
         }
-        if(ucell.lc[1] == 0)
+        if (ucell.lc[1] == 0)
         {
-            move_cell.e21 = 0; move_cell.e22 = 0; move_cell.e23 = 0;
+            move_cell.e21 = 0;
+            move_cell.e22 = 0;
+            move_cell.e23 = 0;
         }
-        if(ucell.lc[2] == 0)
+        if (ucell.lc[2] == 0)
         {
-            move_cell.e31 = 0; move_cell.e32 = 0; move_cell.e33 = 0;
+            move_cell.e31 = 0;
+            move_cell.e32 = 0;
+            move_cell.e33 = 0;
         }
         ucell.latvec += move_cell * (step_size * fac * fac_stress);
 
-        if(PARAM.inp.fixed_axes == "volume")
+        if (PARAM.inp.fixed_axes == "volume")
         {
             double omega_new = std::abs(ucell.latvec.Det()) * pow(ucell.lat0, 3);
-            ucell.latvec *= pow(ucell.omega / omega_new, 1.0/3.0);
+            ucell.latvec *= pow(ucell.omega / omega_new, 1.0 / 3.0);
         }
-        if(PARAM.inp.fixed_ibrav)
+        if (PARAM.inp.fixed_ibrav)
         {
             ucell.remake_cell();
         }
@@ -565,57 +590,57 @@ void Relax::move_cell_ions(UnitCell& ucell,
     // =================================================================
     // Step 2 & 3 : update direct & Cartesian atomic positions
     // =================================================================
-    
+
     // Calculating displacement in Cartesian coordinate (in Angstrom)
     double move_ion[nat * 3];
-    ModuleBase::zeros(move_ion, nat*3);
+    ModuleBase::zeros(move_ion, nat * 3);
 
-    for(int iat=0; iat<nat; iat++)
+    for (int iat = 0; iat < nat; iat++)
     {
-        //Cartesian coordinate
-        //convert from Angstrom to unit of latvec (Bohr)
+        // Cartesian coordinate
+        // convert from Angstrom to unit of latvec (Bohr)
         ModuleBase::Vector3<double> move_ion_cart;
-        move_ion_cart.x = step_size * fac_force * search_dr_ion(iat,0) / ModuleBase::BOHR_TO_A / ucell.lat0;
-        move_ion_cart.y = step_size * fac_force * search_dr_ion(iat,1) / ModuleBase::BOHR_TO_A / ucell.lat0;
-        move_ion_cart.z = step_size * fac_force * search_dr_ion(iat,2) / ModuleBase::BOHR_TO_A / ucell.lat0;
+        move_ion_cart.x = step_size * fac_force * search_dr_ion(iat, 0) / ModuleBase::BOHR_TO_A / ucell.lat0;
+        move_ion_cart.y = step_size * fac_force * search_dr_ion(iat, 1) / ModuleBase::BOHR_TO_A / ucell.lat0;
+        move_ion_cart.z = step_size * fac_force * search_dr_ion(iat, 2) / ModuleBase::BOHR_TO_A / ucell.lat0;
 
-        //convert to Direct coordinate
-        //note here the old GT is used
+        // convert to Direct coordinate
+        // note here the old GT is used
         ModuleBase::Vector3<double> move_ion_dr = move_ion_cart * ucell.GT;
 
         int it = ucell.iat2it[iat];
         int ia = ucell.iat2ia[iat];
         Atom* atom = &ucell.atoms[it];
 
-        if(atom->mbl[ia].x == 1)
+        if (atom->mbl[ia].x == 1)
         {
             move_ion[iat * 3] = move_ion_dr.x * fac;
         }
-        if(atom->mbl[ia].y == 1)
+        if (atom->mbl[ia].y == 1)
         {
             move_ion[iat * 3 + 1] = move_ion_dr.y * fac;
         }
-        if(atom->mbl[ia].z == 1)
+        if (atom->mbl[ia].z == 1)
         {
             move_ion[iat * 3 + 2] = move_ion_dr.z * fac;
         }
     }
 
-    if (ModuleSymmetry::Symmetry::symm_flag && ucell.symm.all_mbl && ucell.symm.nrotk > 0) 
+    if (ModuleSymmetry::Symmetry::symm_flag && ucell.symm.all_mbl && ucell.symm.nrotk > 0)
     {
         ucell.symm.symmetrize_vec3_nat(move_ion);
     }
 
-	ucell.update_pos_taud(move_ion);
+    ucell.update_pos_taud(move_ion);
 
-	// Print the structure file.
-	ucell.print_tau();
+    // Print the structure file.
+    ucell.print_tau();
 
     // =================================================================
     // Step 4 : update G,GT and other stuff
     // =================================================================
 
-    if(if_cell_moves)
+    if (if_cell_moves)
     {
         ucell.a1.x = ucell.latvec.e11;
         ucell.a1.y = ucell.latvec.e12;
@@ -651,8 +676,7 @@ void Relax::move_cell_ions(UnitCell& ucell,
         Parallel_Common::bcast_double(ucell.a3.z);
 #endif
 
-        ucell.omega
-            = std::abs(ucell.latvec.Det()) * ucell.lat0 * ucell.lat0 * ucell.lat0;
+        ucell.omega = std::abs(ucell.latvec.Det()) * ucell.lat0 * ucell.lat0 * ucell.lat0;
 
         ucell.GT = ucell.latvec.Inverse();
         ucell.G = ucell.GT.Transpose();
@@ -663,11 +687,11 @@ void Relax::move_cell_ions(UnitCell& ucell,
     // =================================================================
     // Step 6 : prepare something for next SCF
     // =================================================================
-    //I have a strong feeling that this part should be
-    //at the beginning of the next step (namely 'beforescf'),
-    //but before we have a better organized Esolver
-    //I do not want to change it
-    if(if_cell_moves)
+    // I have a strong feeling that this part should be
+    // at the beginning of the next step (namely 'beforescf'),
+    // but before we have a better organized Esolver
+    // I do not want to change it
+    if (if_cell_moves)
     {
         ucell.setup_cell_after_vc(GlobalV::ofs_running);
         ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SETUP UNITCELL");
diff --git a/source/module_relax/relax_new/relax.h b/source/module_relax/relax_new/relax.h
index bfa9a32bf1..11bd7fdf90 100644
--- a/source/module_relax/relax_new/relax.h
+++ b/source/module_relax/relax_new/relax.h
@@ -1,57 +1,52 @@
-//Wenfei Li, November 2022
-//A new implementation of CG relaxation
+// Wenfei Li, November 2022
+// A new implementation of CG relaxation
 #ifndef RELAX1_H
 #define RELAX1_H
 
+#include "line_search.h"
 #include "module_base/matrix.h"
 #include "module_base/matrix3.h"
-#include "line_search.h"
 #include "module_cell/unitcell.h"
 
 class Relax
 {
-    public:
-
-    Relax(){};
-    ~Relax(){};
+  public:
+    Relax() {};
+    ~Relax() {};
 
-    //prepare for relaxation
+    // prepare for relaxation
     void init_relax(const int nat_in);
-    //perform a single relaxation step
+    // perform a single relaxation step
     bool relax_step(UnitCell& ucell,
-                    const ModuleBase::matrix& force, 
-                    const ModuleBase::matrix &stress, 
+                    const ModuleBase::matrix& force,
+                    const ModuleBase::matrix& stress,
                     const double etot_in);
 
-    private:
+  private:
+    int istep = 0; // count ionic step
 
-    int istep; //count ionic step
+    // setup gradient based on force and stress
+    // constraints are considered here
+    // also check if relaxation has converged
+    // based on threshold in force & stress
+    bool setup_gradient(const UnitCell& ucell, const ModuleBase::matrix& force, const ModuleBase::matrix& stress);
 
-    //setup gradient based on force and stress
-    //constraints are considered here
-    //also check if relaxation has converged
-    //based on threshold in force & stress
-    bool setup_gradient(const UnitCell &ucell,
-                        const ModuleBase::matrix& force, 
-                        const ModuleBase::matrix &stress);
-
-    //check whether previous line search is done
+    // check whether previous line search is done
     bool check_line_search();
 
-    //if line search not done : perform line search
+    // if line search not done : perform line search
     void perform_line_search();
 
-    //if line search done: find new search direction and make a trial move
+    // if line search done: find new search direction and make a trial move
     void new_direction();
 
-    //move ions and lattice vectors
-    void move_cell_ions(UnitCell& ucell,
-                        const bool is_new_dir);
+    // move ions and lattice vectors
+    void move_cell_ions(UnitCell& ucell, const bool is_new_dir);
 
-    int nat; // number of atoms
-    bool ltrial; // if last step is trial step
+    int nat = 0;         // number of atoms
+    bool ltrial = false; // if last step is trial step
 
-    double step_size;
+    double step_size = 0.0;
 
     // Gradients; _p means previous step
     ModuleBase::matrix grad_ion;
@@ -66,32 +61,41 @@ class Relax
     ModuleBase::matrix search_dr_cell_p;
 
     // Used for applyting constraints
-    bool if_cell_moves;
+    bool if_cell_moves = false;
 
-    //Keeps track of how many CG trial steps have been performed,
-    //namely the number of CG directions followed
-    //Note : this should not be confused with number of ionic steps
-    //which includes both trial and line search steps
-    int cg_step;
+    // Keeps track of how many CG trial steps have been performed,
+    // namely the number of CG directions followed
+    // Note : this should not be confused with number of ionic steps
+    // which includes both trial and line search steps
+    int cg_step = 0;
 
-    //in CG, search_dr = search_dr_p + grad * gamma
-    double gamma;
+    // in CG, search_dr = search_dr_p + grad * gamma
+    double gamma = 0.0;
     void calculate_gamma();
 
-    //Intermediate variables
-    //I put them here because they are used across different subroutines
-    double sr_sr, srp_srp; //inner/cross products between search directions
-    double gr_gr, gr_grp, grp_grp; //inner/cross products between gradients
-    double gr_sr; //cross product between search direction and gradient
-    double e1ord1, e1ord2, e2ord, e2ord2;
-    double dmove,dmovel,dmoveh;
-    double etot, etot_p;
-    double force_thr_eva;
-
-    bool brent_done; //if brent line search is finished
-
-    double fac_force;
-    double fac_stress;
+    // Intermediate variables
+    // I put them here because they are used across different subroutines
+    double sr_sr = 0.0;
+    double srp_srp = 0.0; // inner/cross products between search directions
+    double gr_gr = 0.0;
+    double gr_grp = 0.0;
+    double grp_grp = 0.0; // inner/cross products between gradients
+    double gr_sr = 0.0;   // cross product between search direction and gradient
+    double e1ord1 = 0.0;
+    double e1ord2 = 0.0;
+    double e2ord = 0.0;
+    double e2ord2 = 0.0;
+    double dmove = 0.0;
+    double dmovel = 0.0;
+    double dmoveh = 0.0;
+    double etot = 0.0;
+    double etot_p = 0.0;
+    double force_thr_eva = 0.0;
+
+    bool brent_done = false; // if brent line search is finished
+
+    double fac_force = 0.0;
+    double fac_stress = 0.0;
 
     ModuleBase::Matrix3 latvec_save;
     Line_Search ls;