diff --git a/source/module_cell/module_neighbor/sltk_grid.cpp b/source/module_cell/module_neighbor/sltk_grid.cpp
index d20f9bb625..6fd6575c23 100644
--- a/source/module_cell/module_neighbor/sltk_grid.cpp
+++ b/source/module_cell/module_neighbor/sltk_grid.cpp
@@ -152,12 +152,19 @@ void Grid::setMemberVariables(std::ofstream& ofs_in, //  output data to ofs
     this->clear_atoms();
 
     // random selection, in order to estimate again.
-    this->x_min = ucell.atoms[0].tau[0].x;
-    this->y_min = ucell.atoms[0].tau[0].y;
-    this->z_min = ucell.atoms[0].tau[0].z;
-    this->x_max = ucell.atoms[0].tau[0].x;
-    this->y_max = ucell.atoms[0].tau[0].y;
-    this->z_max = ucell.atoms[0].tau[0].z;
+    for (int it = 0; it < ucell.ntype; it++)
+    {
+        if (ucell.atoms[it].na > 0)
+        {
+            this->x_min = ucell.atoms[it].tau[0].x;
+            this->y_min = ucell.atoms[it].tau[0].y;
+            this->z_min = ucell.atoms[it].tau[0].z;
+            this->x_max = ucell.atoms[it].tau[0].x;
+            this->y_max = ucell.atoms[it].tau[0].y;
+            this->z_max = ucell.atoms[it].tau[0].z;
+            break;
+        }
+    }
 
     ModuleBase::Vector3<double> vec1(ucell.latvec.e11, ucell.latvec.e12, ucell.latvec.e13);
     ModuleBase::Vector3<double> vec2(ucell.latvec.e21, ucell.latvec.e22, ucell.latvec.e23);
diff --git a/source/module_cell/update_cell.cpp b/source/module_cell/update_cell.cpp
index 1ea761afba..c4c324ce45 100644
--- a/source/module_cell/update_cell.cpp
+++ b/source/module_cell/update_cell.cpp
@@ -384,41 +384,21 @@ void periodic_boundary_adjustment(Atom* atoms,
     // first adjust direct coordinates,
     // then update them into cartesian coordinates,
     //----------------------------------------------
-    for (int i=0;i<ntype;i++) {
-        Atom* atom = &atoms[i];
-        for (int j=0;j<atom->na;j++) {
-            printf("the taud is %f %f %f\n",atom->taud[j].x,atom->taud[j].y,atom->taud[j].z);
-        }
-    }
     for (int it = 0; it < ntype; it++) {
         Atom* atom = &atoms[it];
         for (int ia = 0; ia < atom->na; ia++) {
             // mohan update 2011-03-21
-            if (atom->taud[ia].x < 0)
+            for (int ik = 0; ik < 3; ik++) 
             {
-                atom->taud[ia].x += 1.0;
-            }
-            if (atom->taud[ia].y < 0)
-            {
-                atom->taud[ia].y += 1.0;
-            }
-            if (atom->taud[ia].z < 0)
-            {
-                atom->taud[ia].z += 1.0;
-            }
-            if (atom->taud[ia].x >= 1.0)
-            {
-                atom->taud[ia].x -= 1.0;
-            }
-            if (atom->taud[ia].y >= 1.0)
-            {
-                atom->taud[ia].y -= 1.0;
-            }
-            if (atom->taud[ia].z >= 1.0)
-            {
-                atom->taud[ia].z -= 1.0;
+                if (atom->taud[ia][ik] < 0) 
+                {
+                    atom->taud[ia][ik] += 1.0;
+                }
+                if (atom->taud[ia][ik] >= 1.0) 
+                {
+                    atom->taud[ia][ik] -= 1.0;
+                }
             }
-
             if (atom->taud[ia].x < 0 
                 || atom->taud[ia].y < 0
                 || atom->taud[ia].z < 0 
diff --git a/source/module_elecstate/module_charge/charge_init.cpp b/source/module_elecstate/module_charge/charge_init.cpp
index 6a8e59bee8..9d104312b7 100644
--- a/source/module_elecstate/module_charge/charge_init.cpp
+++ b/source/module_elecstate/module_charge/charge_init.cpp
@@ -37,6 +37,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
     this->pgrid = &pgrid;
 
     bool read_error = false;
+    bool read_kin_error = false;
     if (PARAM.inp.init_chg == "file" || PARAM.inp.init_chg == "auto")
     {
         GlobalV::ofs_running << " try to read charge from file" << std::endl;
@@ -101,72 +102,100 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
             }
         }
 
-        if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
+        if (read_error)
         {
-            GlobalV::ofs_running << " try to read kinetic energy density from file" << std::endl;
-            // try to read charge from binary file first, which is the same as QE
-            std::vector<std::complex<double>> kin_g_space(PARAM.inp.nspin * this->ngmc, {0.0, 0.0});
-            std::vector<std::complex<double>*> kin_g;
-            for (int is = 0; is < PARAM.inp.nspin; is++)
+            const std::string warn_msg
+                = " WARNING: \"init_chg\" is enabled but ABACUS failed to read charge density from file.\n"
+                  " Please check if there is SPINX_CHG.cube (X=1,...) or {suffix}-CHARGE-DENSITY.restart in the "
+                  "directory.\n";
+            std::cout << std::endl << warn_msg;
+            if (PARAM.inp.init_chg == "file")
             {
-                kin_g.push_back(kin_g_space.data() + is * this->ngmc);
+                ModuleBase::WARNING_QUIT("Charge::init_rho",
+                                         "Failed to read in charge density from file.\nIf you want to use atomic "
+                                         "charge initialization, \nplease set init_chg to atomic in INPUT.");
             }
+        }
 
-            std::stringstream binary;
-            binary << PARAM.globalv.global_readin_dir << PARAM.inp.suffix + "-TAU-DENSITY.restart";
-            if (ModuleIO::read_rhog(binary.str(), rhopw, kin_g.data()))
+        if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
+        {
+            // If the charge density is not read in, then the kinetic energy density is not read in either
+            if (!read_error)
             {
-                GlobalV::ofs_running << " Read in the kinetic energy density: " << binary.str() << std::endl;
-                for (int is = 0; is < PARAM.inp.nspin; ++is)
+                GlobalV::ofs_running << " try to read kinetic energy density from file" << std::endl;
+                // try to read charge from binary file first, which is the same as QE
+                std::vector<std::complex<double>> kin_g_space(PARAM.inp.nspin * this->ngmc, {0.0, 0.0});
+                std::vector<std::complex<double>*> kin_g;
+                for (int is = 0; is < PARAM.inp.nspin; is++)
                 {
-                    rhopw->recip2real(kin_g[is], this->kin_r[is]);
+                    kin_g.push_back(kin_g_space.data() + is * this->ngmc);
                 }
-            }
-            else
-            {
-                for (int is = 0; is < PARAM.inp.nspin; is++)
+
+                std::stringstream binary;
+                binary << PARAM.globalv.global_readin_dir << PARAM.inp.suffix + "-TAU-DENSITY.restart";
+                if (ModuleIO::read_rhog(binary.str(), rhopw, kin_g.data()))
                 {
-                    std::stringstream ssc;
-                    ssc << PARAM.globalv.global_readin_dir << "SPIN" << is + 1 << "_TAU.cube";
-                    // mohan update 2012-02-10, sunliang update 2023-03-09
-                    if (ModuleIO::read_vdata_palgrid(
-                            pgrid,
-                            (PARAM.inp.esolver_type == "sdft" ? GlobalV::RANK_IN_STOGROUP : GlobalV::MY_RANK),
-                            GlobalV::ofs_running,
-                            ssc.str(),
-                            this->kin_r[is],
-                            ucell.nat))
+                    GlobalV::ofs_running << " Read in the kinetic energy density: " << binary.str() << std::endl;
+                    for (int is = 0; is < PARAM.inp.nspin; ++is)
                     {
-                        GlobalV::ofs_running << " Read in the kinetic energy density: " << ssc.str() << std::endl;
+                        rhopw->recip2real(kin_g[is], this->kin_r[is]);
                     }
                 }
+                else
+                {
+                    for (int is = 0; is < PARAM.inp.nspin; is++)
+                    {
+                        std::stringstream ssc;
+                        ssc << PARAM.globalv.global_readin_dir << "SPIN" << is + 1 << "_TAU.cube";
+                        // mohan update 2012-02-10, sunliang update 2023-03-09
+                        if (ModuleIO::read_vdata_palgrid(
+                                pgrid,
+                                (PARAM.inp.esolver_type == "sdft" ? GlobalV::RANK_IN_STOGROUP : GlobalV::MY_RANK),
+                                GlobalV::ofs_running,
+                                ssc.str(),
+                                this->kin_r[is],
+                                ucell.nat))
+                        {
+                            GlobalV::ofs_running << " Read in the kinetic energy density: " << ssc.str() << std::endl;
+                        }
+                        else
+                        {
+                            read_kin_error = true;
+                            std::cout << " WARNING: \"init_chg\" is enabled but ABACUS failed to read kinetic energy "
+                                         "density from file.\n"
+                                         " Please check if there is SPINX_TAU.cube (X=1,...) or "
+                                         "{suffix}-TAU-DENSITY.restart in the directory.\n"
+                                      << std::endl;
+                            break;
+                        }
+                    }
+                }
+            }
+            else
+            {
+                read_kin_error = true;
             }
         }
     }
-    if (read_error)
+
+    if (PARAM.inp.init_chg == "atomic" || read_error) // mohan add 2007-10-17
     {
-        const std::string warn_msg = " WARNING: \"init_chg\" is enabled but ABACUS failed to read charge density from file.\n"
-                                     " Please check if there is SPINX_CHG.cube (X=1,...) or {suffix}-CHARGE-DENSITY.restart in the directory.\n";
-        std::cout << std::endl << warn_msg;
-        if (PARAM.inp.init_chg == "auto")
+        if (read_error)
         {
-            std::cout << " Charge::init_rho: use atomic initialization instead." << std::endl << std::endl;
-        }
-        else if (PARAM.inp.init_chg == "file")
-        {
-            ModuleBase::WARNING_QUIT("Charge::init_rho", "Failed to read in charge density from file.\nIf you want to use atomic charge initialization, \nplease set init_chg to atomic in INPUT.");
+            std::cout << " Charge::init_rho: use atomic initialization instead." << std::endl;
         }
+        this->atomic_rho(PARAM.inp.nspin, ucell.omega, rho, strucFac, ucell);
     }
 
-    if (PARAM.inp.init_chg == "atomic" || 
-        (PARAM.inp.init_chg == "auto" && read_error)) // mohan add 2007-10-17
+    // wenfei 2021-7-29 : initial tau = 3/5 rho^2/3, Thomas-Fermi
+    if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
     {
-        this->atomic_rho(PARAM.inp.nspin, ucell.omega, rho, strucFac, ucell);
-
-        // liuyu 2023-06-29 : move here from atomic_rho(), which will be called several times in charge extrapolation
-        // wenfei 2021-7-29 : initial tau = 3/5 rho^2/3, Thomas-Fermi
-        if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
+        if (PARAM.inp.init_chg == "atomic" || read_kin_error)
         {
+            if (read_kin_error)
+            {
+                std::cout << " Charge::init_rho: init kinetic energy density from rho." << std::endl;
+            }
             const double pi = 3.141592653589790;
             const double fact = (3.0 / 5.0) * pow(3.0 * pi * pi, 2.0 / 3.0);
             for (int is = 0; is < PARAM.inp.nspin; ++is)
diff --git a/source/module_esolver/esolver_gets.cpp b/source/module_esolver/esolver_gets.cpp
index cf51f6bbbe..842b0340e1 100644
--- a/source/module_esolver/esolver_gets.cpp
+++ b/source/module_esolver/esolver_gets.cpp
@@ -147,7 +147,7 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep)
 }
 
 void ESolver_GetS::after_all_runners(UnitCell& ucell) {};
-double ESolver_GetS::cal_energy() {};
+double ESolver_GetS::cal_energy() { return 0.0; };
 void ESolver_GetS::cal_force(UnitCell& ucell, ModuleBase::matrix& force) {};
 void ESolver_GetS::cal_stress(UnitCell& ucell, ModuleBase::matrix& stress) {};
 
diff --git a/source/module_esolver/esolver_ks.h b/source/module_esolver/esolver_ks.h
index c20a53f555..16932bcc6a 100644
--- a/source/module_esolver/esolver_ks.h
+++ b/source/module_esolver/esolver_ks.h
@@ -43,7 +43,7 @@ class ESolver_KS : public ESolver_FP
     virtual void iter_init(UnitCell& ucell, const int istep, const int iter);
 
     //! Something to do after hamilt2density function in each iter loop.
-    virtual void iter_finish(UnitCell& ucell, const int istep, int& iter);
+    virtual void iter_finish(UnitCell& ucell, const int istep, int& iter) override;
 
     // calculate electron density from a specific Hamiltonian with ethr
     virtual void hamilt2density_single(UnitCell& ucell, const int istep, const int iter, const double ethr);
diff --git a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.cpp b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.cpp
index 7207a64a27..8f2c4d8779 100644
--- a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.cpp
+++ b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.cpp
@@ -69,7 +69,7 @@ void LCAO_Deepks_Interface<TK, TR>::out_deepks_labels(const double& etot,
 
         if (PARAM.inp.deepks_bandgap)
         {
-            const int nocc = PARAM.inp.nelec / 2;
+            const int nocc = (PARAM.inp.nelec+1) / 2;
             std::vector<double> o_tot(nks);
             for (int iks = 0; iks < nks; ++iks)
             {
diff --git a/source/module_hsolver/diago_iter_assist.cpp b/source/module_hsolver/diago_iter_assist.cpp
index 5ec443ab4e..89aa94d1a9 100644
--- a/source/module_hsolver/diago_iter_assist.cpp
+++ b/source/module_hsolver/diago_iter_assist.cpp
@@ -118,6 +118,9 @@ void DiagoIterAssist<T, Device>::diagH_subspace(const hamilt::Hamilt<T, Device>*
     // after generation of H and S matrix, diag them
     DiagoIterAssist::diagH_LAPACK(nstart, n_band, hcc, scc, nstart, en, vcc);
 
+
+    const int ld_temp = in_place ? dmax : dmin;
+
     { // code block to calculate evc
         gemm_op<T, Device>()(ctx,
                              'N',
@@ -132,12 +135,12 @@ void DiagoIterAssist<T, Device>::diagH_subspace(const hamilt::Hamilt<T, Device>*
                              nstart,
                              &zero,
                              temp,
-                             dmin);
+                             ld_temp);
     }
 
     if (!in_place)
     {
-        matrixSetToAnother<T, Device>()(ctx, n_band, temp, dmin, evc.get_pointer(), dmax);
+        matrixSetToAnother<T, Device>()(ctx, n_band, temp, ld_temp, evc.get_pointer(), dmax);
         delmem_complex_op()(ctx, temp);
     }
     delmem_complex_op()(ctx, hcc);
diff --git a/source/module_io/get_pchg_lcao.cpp b/source/module_io/get_pchg_lcao.cpp
index 6e069fd017..5a1b759866 100644
--- a/source/module_io/get_pchg_lcao.cpp
+++ b/source/module_io/get_pchg_lcao.cpp
@@ -200,7 +200,9 @@ void IState_Charge::begin(Gint_k& gk,
         mode = 3;
     }
 
-    const int fermi_band = static_cast<int>((nelec + 1) / 2 + 1.0e-8);
+    int fermi_band = static_cast<int>((nelec + 1) / 2 + 1.0e-8);
+    // for nspin = 4 case, fermi_band is the number of electrons because of each band has 1 electron
+    if(PARAM.inp.nspin == 4) fermi_band = nelec;
     std::cout << " number of electrons = " << nelec << std::endl;
     std::cout << " number of occupied bands = " << fermi_band << std::endl;
 
@@ -216,7 +218,7 @@ void IState_Charge::begin(Gint_k& gk,
             elecstate::DensityMatrix<std::complex<double>, double> DM(this->ParaV, nspin_dm, kv.kvec_d, kv.get_nks() / nspin_dm);
 
 #ifdef __MPI
-            this->idmatrix(ib, nspin, nelec, nlocal, wg, DM, kv, if_separate_k);
+            this->idmatrix(ib, nspin_dm, nelec, nlocal, wg, DM, kv, if_separate_k);
 #else
             ModuleBase::WARNING_QUIT("IState_Charge::begin", "The `pchg` calculation is only available for MPI now!");
 #endif
@@ -224,7 +226,7 @@ void IState_Charge::begin(Gint_k& gk,
             if (if_separate_k)
             {
                 // For multi-k, loop over all real k-points
-                for (int ik = 0; ik < kv.get_nks() / nspin; ++ik)
+                for (int ik = 0; ik < kv.get_nks() / nspin_dm; ++ik)
                 {
                     for (int is = 0; is < nspin; ++is)
                     {
@@ -506,16 +508,8 @@ void IState_Charge::idmatrix(const int& ib,
     ModuleBase::TITLE("IState_Charge", "idmatrix");
     assert(wg.nr == kv.get_nks());
 
-    const int fermi_band = static_cast<int>((nelec + 1) / 2 + 1.0e-8);
-
     // To ensure the normalization of charge density in multi-k calculation (if if_separate_k is true)
-    double wg_sum_k = 0;
-    double wg_sum_k_homo = 0;
-    for (int ik = 0; ik < kv.get_nks() / nspin; ++ik)
-    {
-        wg_sum_k += wg(ik, ib);
-        wg_sum_k_homo += wg(ik, fermi_band - 1);
-    }
+    double wg_sum_k = 1.0;
 
     for (int ik = 0; ik < kv.get_nks(); ++ik)
     {
@@ -530,11 +524,11 @@ void IState_Charge::idmatrix(const int& ib,
             double wg_value;
             if (if_separate_k)
             {
-                wg_value = (ib < fermi_band) ? wg_sum_k : wg_sum_k_homo;
+                wg_value = wg_sum_k;
             }
             else
             {
-                wg_value = (ib < fermi_band) ? wg(ik, ib) : wg(ik, fermi_band - 1);
+                wg_value = wg(ik, 0);
             }
             wg_local[ib_local] = wg_value;
         }
diff --git a/source/module_io/input_conv.h b/source/module_io/input_conv.h
index a756cc5499..e14c10dac7 100644
--- a/source/module_io/input_conv.h
+++ b/source/module_io/input_conv.h
@@ -12,8 +12,8 @@
 #include <iomanip>
 #include <iostream>
 #include <regex.h>
-#include <stdio.h>
-#include <string.h>
+#include <cstdio>
+#include <cstring>
 #include <string>
 #include <vector>
 #include <algorithm>
@@ -31,7 +31,7 @@ void tmp_convert();
  * @brief Pass the data members from the INPUT instance(defined in
  * module_io/input.cpp) to GlobalV and GlobalC.
  */
-void Convert(void);
+void Convert();
 
 /**
  * @brief To parse input parameters as expressions into vectors
@@ -47,7 +47,7 @@ void parse_expression(const std::string& fn, std::vector<T>& vec)
 {
     ModuleBase::TITLE("Input_Conv", "parse_expression");
     int count = 0;
-    std::string pattern("([0-9]+\\*[0-9.]+|[0-9,.]+)");
+    std::string pattern("([-+]?[0-9]+\\*[-+]?[0-9.]+|[-+]?[0-9,.]+)");
     std::vector<std::string> str;
     std::stringstream ss(fn);
     std::string section;
@@ -103,6 +103,7 @@ void parse_expression(const std::string& fn, std::vector<T>& vec)
         {
             size_t pos = sub_str.find("*");
             int num = stoi(sub_str.substr(0, pos));
+            assert(num>=0);
             T occ = stof(sub_str.substr(pos + 1, sub_str.size()));
             // std::vector<double> ocp_temp(num, occ);
             // const std::vector<double>::iterator dest = vec.begin() + count;
diff --git a/source/module_io/output_mulliken.h b/source/module_io/output_mulliken.h
index a5c0c21bbc..4611c124d1 100644
--- a/source/module_io/output_mulliken.h
+++ b/source/module_io/output_mulliken.h
@@ -138,7 +138,7 @@ void cal_mag(Parallel_Orbitals* pv,
             auto sc_lambda
                 = new hamilt::DeltaSpin<hamilt::OperatorLCAO<TK, double>>(nullptr,
                                                                           kv.kvec_d,
-                                                                          nullptr,
+                                                                          dynamic_cast<hamilt::HamiltLCAO<TK, double>*>(p_ham)->getHR(),
                                                                           ucell,
                                                                           &gd,
                                                                           two_center_bundle.overlap_orb_onsite.get(),
@@ -164,7 +164,7 @@ void cal_mag(Parallel_Orbitals* pv,
             auto sc_lambda = new hamilt::DeltaSpin<hamilt::OperatorLCAO<std::complex<double>, std::complex<double>>>(
                 nullptr,
                 kv.kvec_d,
-                nullptr,
+                dynamic_cast<hamilt::HamiltLCAO<std::complex<double>, std::complex<double>>*>(p_ham)->getHR(),
                 ucell,
                 &gd,
                 two_center_bundle.overlap_orb_onsite.get(),
diff --git a/source/module_lr/AX/AX.h b/source/module_lr/AX/AX.h
index cb5e9f18f6..c1cdac3532 100644
--- a/source/module_lr/AX/AX.h
+++ b/source/module_lr/AX/AX.h
@@ -13,13 +13,13 @@ namespace LR
         const psi::Psi<double>& c,
         const int& nocc,
         const int& nvirt,
-        double* AX_istate);
+        double* const AX_istate);
     void cal_AX_blas(
         const std::vector<container::Tensor>& V_istate,
         const psi::Psi<double>& c,
         const int& nocc,
         const int& nvirt,
-        double* AX_istate,
+        double* const AX_istate,
         const bool add_on = true);
 #ifdef __MPI
     void cal_AX_pblas(
@@ -31,7 +31,7 @@ namespace LR
         const int& nocc,
         const int& nvirt,
         const Parallel_2D& pX,
-        double* AX_istate,
+        double* const AX_istate,
         const bool add_on=true);
 #endif
     // complex
@@ -40,13 +40,13 @@ namespace LR
         const psi::Psi<std::complex<double>>& c,
         const int& nocc,
         const int& nvirt,
-        std::complex<double>* AX_istate);
+        std::complex<double>* const AX_istate);
     void cal_AX_blas(
         const std::vector<container::Tensor>& V_istate,
         const psi::Psi<std::complex<double>>& c,
         const int& nocc,
         const int& nvirt,
-        std::complex<double>* AX_istate,
+        std::complex<double>* const AX_istate,
         const bool add_on = true);
 
 #ifdef __MPI
@@ -59,7 +59,7 @@ namespace LR
         const int& nocc,
         const int& nvirt,
         const Parallel_2D& pX,
-        std::complex<double>* AX_istate,
+        std::complex<double>* const AX_istate,
         const bool add_on = true);
 #endif
 }
\ No newline at end of file
diff --git a/source/module_lr/AX/AX_parallel.cpp b/source/module_lr/AX/AX_parallel.cpp
index 221962416e..c5b956bac6 100644
--- a/source/module_lr/AX/AX_parallel.cpp
+++ b/source/module_lr/AX/AX_parallel.cpp
@@ -71,7 +71,7 @@ namespace LR
         const int& nocc,
         const int& nvirt,
         const Parallel_2D& pX,
-        std::complex<double>* AX_istate,
+        std::complex<double>* const AX_istate,
         const bool add_on)
     {
         ModuleBase::TITLE("hamilt_lrtd", "cal_AX_plas");
diff --git a/source/module_lr/AX/AX_serial.cpp b/source/module_lr/AX/AX_serial.cpp
index abe990e5ff..234d4412dd 100644
--- a/source/module_lr/AX/AX_serial.cpp
+++ b/source/module_lr/AX/AX_serial.cpp
@@ -40,7 +40,7 @@ namespace LR
         const psi::Psi<std::complex<double>>& c,
         const int& nocc,
         const int& nvirt,
-        std::complex<double>* AX_istate)
+        std::complex<double>* const AX_istate)
     {
         ModuleBase::TITLE("hamilt_lrtd", "cal_AX_forloop");
         const int nks = V_istate.size();
@@ -107,7 +107,7 @@ namespace LR
         const psi::Psi<std::complex<double>>& c,
         const int& nocc,
         const int& nvirt,
-        std::complex<double>* AX_istate,
+        std::complex<double>* const AX_istate,
         const bool add_on)
     {
         ModuleBase::TITLE("hamilt_lrtd", "cal_AX_blas");
diff --git a/source/module_lr/dm_trans/dm_trans.h b/source/module_lr/dm_trans/dm_trans.h
index 385d33defb..77cb0870f1 100644
--- a/source/module_lr/dm_trans/dm_trans.h
+++ b/source/module_lr/dm_trans/dm_trans.h
@@ -14,7 +14,7 @@ namespace LR
 /// \f[ \tilde{\rho}_{\mu_j\mu_b}=\sum_{jb}c_{j,\mu_j}X_{jb}c^*_{b,\mu_b} \f]
     template<typename T>
     std::vector<container::Tensor> cal_dm_trans_pblas(
-        const T* X_istate,
+        const T* const X_istate,
         const Parallel_2D& px,
         const psi::Psi<T>& c,
         const Parallel_2D& pc,
@@ -29,7 +29,7 @@ namespace LR
     /// @brief calculate the 2d-block transition density matrix in AO basis using ?gemm
     template<typename T>
     std::vector<container::Tensor> cal_dm_trans_blas(
-        const T* X_istate,
+        const T* const X_istate,
         const psi::Psi<T>& c,
         const int& nocc, const int& nvirt,
         const bool renorm_k = true,
@@ -39,7 +39,7 @@ namespace LR
     /// @brief calculate the 2d-block transition density matrix in AO basis using for loop (for test)
     template<typename T>
     std::vector<container::Tensor> cal_dm_trans_forloop_serial(
-        const T* X_istate,
+        const T* const X_istate,
         const psi::Psi<T>& c,
         const int& nocc, const int& nvirt,
         const bool renorm_k = true,
diff --git a/source/module_lr/dm_trans/dm_trans_parallel.cpp b/source/module_lr/dm_trans/dm_trans_parallel.cpp
index e230ddbc00..701bd87621 100644
--- a/source/module_lr/dm_trans/dm_trans_parallel.cpp
+++ b/source/module_lr/dm_trans/dm_trans_parallel.cpp
@@ -10,7 +10,7 @@ namespace LR
     // c: nao*nbands in para2d, nbands*nao in psi  (row-para and constructed: nao)
     // X: nvirt*nocc in para2d, nocc*nvirt in psi (row-para and constructed: nvirt)
 template <>
-std::vector<container::Tensor> cal_dm_trans_pblas(const double* X_istate,
+std::vector<container::Tensor> cal_dm_trans_pblas(const double* const X_istate,
     const Parallel_2D& px,
     const psi::Psi<double>& c,
     const Parallel_2D& pc,
@@ -62,7 +62,7 @@ std::vector<container::Tensor> cal_dm_trans_pblas(const double* X_istate,
     return dm_trans;
 }
 template <>
-std::vector<container::Tensor> cal_dm_trans_pblas(const std::complex<double>* X_istate,
+std::vector<container::Tensor> cal_dm_trans_pblas(const std::complex<double>* const X_istate,
     const Parallel_2D& px,
     const psi::Psi<std::complex<double>>& c,
     const Parallel_2D& pc,
diff --git a/source/module_lr/dm_trans/dm_trans_serial.cpp b/source/module_lr/dm_trans/dm_trans_serial.cpp
index 70dc6e6542..2a098c08c4 100644
--- a/source/module_lr/dm_trans/dm_trans_serial.cpp
+++ b/source/module_lr/dm_trans/dm_trans_serial.cpp
@@ -6,7 +6,7 @@
 namespace LR
 {
     template<> std::vector<container::Tensor> cal_dm_trans_forloop_serial(
-        const double* X_istate,
+        const double* const X_istate,
         const psi::Psi<double>& c,
         const int& nocc,
         const int& nvirt,
@@ -41,7 +41,7 @@ namespace LR
     }
 
     template<> std::vector<container::Tensor> cal_dm_trans_forloop_serial(
-        const std::complex<double>* X_istate,
+        const std::complex<double>* const X_istate,
         const psi::Psi<std::complex<double>>& c,
         const int& nocc,
         const int& nvirt,
@@ -78,7 +78,7 @@ namespace LR
 
 
     template<> std::vector<container::Tensor> cal_dm_trans_blas(
-        const double* X_istate,
+        const double* const X_istate,
         const psi::Psi<double>& c,
         const int& nocc,
         const int& nvirt,
@@ -112,7 +112,7 @@ namespace LR
 
 
     template<> std::vector<container::Tensor> cal_dm_trans_blas(
-        const std::complex<double>* X_istate,
+        const std::complex<double>* const X_istate,
         const psi::Psi<std::complex<double>>& c,
         const int& nocc,
         const int& nvirt,
diff --git a/source/module_lr/esolver_lrtd_lcao.cpp b/source/module_lr/esolver_lrtd_lcao.cpp
index 8c00ef2fdd..97842897b9 100644
--- a/source/module_lr/esolver_lrtd_lcao.cpp
+++ b/source/module_lr/esolver_lrtd_lcao.cpp
@@ -568,7 +568,7 @@ void LR::ESolver_LR<T, TR>::after_all_runners(UnitCell& ucell)
         {
             spectrum.optical_absorption_method1(freq, input.abs_broadening);
             // =============================================== for test ====================================================
-            // spectrum.optical_absorption_method2(freq, input.abs_broadening, spin_types[is]);
+            // spectrum.optical_absorption_method2(freq, input.abs_broadening);
             // spectrum.test_transition_dipoles_velocity_ks(eig_ks.c);
             // spectrum.write_transition_dipole(PARAM.globalv.global_out_dir + "dipole_velocity_ks.dat");
             // =============================================== for test ====================================================
diff --git a/source/module_lr/hamilt_casida.h b/source/module_lr/hamilt_casida.h
index 47eeb05a7a..b28ad9fde8 100644
--- a/source/module_lr/hamilt_casida.h
+++ b/source/module_lr/hamilt_casida.h
@@ -115,7 +115,7 @@ namespace LR
             }
 #endif
 
-            this->cal_dm_trans = [&, this](const int& is, const T* X)->void
+            this->cal_dm_trans = [&, this](const int& is, const T* const X)->void
                 {
                     const auto psi_ks_is = LR_Util::get_psi_spin(psi_ks_in, is, nk);
 #ifdef __MPI
@@ -134,7 +134,7 @@ namespace LR
 
         std::vector<T> matrix()const;
 
-        void hPsi(const T* psi_in, T* hpsi, const int ld_psi, const int& nband) const
+        void hPsi(const T* const psi_in, T* const hpsi, const int ld_psi, const int& nband) const
         {
             assert(ld_psi == nk * pX[0].get_local_size());
             for (int ib = 0;ib < nband;++ib)
@@ -189,6 +189,6 @@ namespace LR
         /// first node operator, add operations from each operators
         hamilt::Operator<T, base_device::DEVICE_CPU>* ops = nullptr;
 
-        std::function<void(const int&, const T*)> cal_dm_trans;
+        std::function<void(const int&, const T* const)> cal_dm_trans;
     };
 }
diff --git a/source/module_lr/hamilt_ulr.hpp b/source/module_lr/hamilt_ulr.hpp
index 77ee62ce73..c058d14803 100644
--- a/source/module_lr/hamilt_ulr.hpp
+++ b/source/module_lr/hamilt_ulr.hpp
@@ -68,7 +68,7 @@ namespace LR
             }
 #endif
 
-            this->cal_dm_trans = [&, this](const int& is, const T* X)->void
+            this->cal_dm_trans = [&, this](const int& is, const T* const X)->void
                 {
                     const auto psi_ks_is = LR_Util::get_psi_spin(psi_ks_in, is, nk);
                     // LR_Util::print_value(X, pX_in[is].get_local_size());
@@ -88,7 +88,7 @@ namespace LR
         {
             for (auto& op : ops) { delete op; }
         }
-        void hPsi(const T* psi_in, T* hpsi, const int ld_psi, const int& nband) const
+        void hPsi(const T* const psi_in, T* const hpsi, const int ld_psi, const int& nband) const
         {
             ModuleBase::TITLE("HamiltULR", "hPsi");
             assert(ld_psi == this->ldim);
@@ -224,7 +224,7 @@ namespace LR
         /// Hxc+Exx: size=nbands, store the result of each bands for common use
         std::unique_ptr<elecstate::DensityMatrix<T, T>> DM_trans;
 
-        std::function<void(const int&, const T*)> cal_dm_trans;
+        std::function<void(const int&, const T* const)> cal_dm_trans;
         const bool tdm_sym = false;     ///< whether to symmetrize the transition density matrix
     };
 }
\ No newline at end of file
diff --git a/source/module_lr/lr_spectrum.cpp b/source/module_lr/lr_spectrum.cpp
index c6db92c67b..0a2c5735f6 100644
--- a/source/module_lr/lr_spectrum.cpp
+++ b/source/module_lr/lr_spectrum.cpp
@@ -148,7 +148,8 @@ template<> double LR::LR_Spectrum<double>::cal_mean_squared_dipole(ModuleBase::V
 }
 template<> double LR::LR_Spectrum<std::complex<double>>::cal_mean_squared_dipole(ModuleBase::Vector3<std::complex<double>> dipole)
 {
-    return dipole.norm2().real() / 3.;
+    // return dipole.norm2().real() / 3.;       // ModuleBase::Vector3::norm2 calculates x*x + y*y + z*z, but here we need x*x.conj() + y*y.conj() + z*z.conj()
+    return (std::norm(dipole.x) + std::norm(dipole.y) + std::norm(dipole.z)) / 3.;
 }
 
 template<typename T>
@@ -217,12 +218,12 @@ void LR::LR_Spectrum<T>::transition_analysis(const std::string& spintype)
     ofs << std::setw(40) << spintype << std::endl;
     ofs << "==================================================================== " << std::endl;
     ofs << std::setw(8) << "State" << std::setw(30) << "Excitation Energy (Ry, eV)" <<
-        std::setw(45) << "Transition dipole x, y, z (a.u.)" << std::setw(30) << "Oscillator strength(a.u.)" << std::endl;
+        std::setw(90) << "Transition dipole x, y, z (a.u.)" << std::setw(30) << "Oscillator strength(a.u.)" << std::endl;
     ofs << "------------------------------------------------------------------------------------ " << std::endl;
     for (int istate = 0;istate < nstate;++istate)
         ofs << std::setw(8) << istate << std::setw(15) << std::setprecision(6) << eig[istate] << std::setw(15) << eig[istate] * ModuleBase::Ry_to_eV
-        << std::setw(15) << transition_dipole_[istate].x << std::setw(15) << transition_dipole_[istate].y << std::setw(15) << transition_dipole_[istate].z
-        << std::setw(30) << oscillator_strength_[istate] << std::endl;
+        << std::setprecision(4) << std::setw(30) << transition_dipole_[istate].x << std::setw(30) << transition_dipole_[istate].y << std::setw(30) << transition_dipole_[istate].z
+        << std::setprecision(6) << std::setw(30) << oscillator_strength_[istate] << std::endl;
     ofs << "------------------------------------------------------------------------------------ " << std::endl;
     ofs << std::setw(8) << "State" << std::setw(20) << "Occupied orbital"
         << std::setw(20) << "Virtual orbital" << std::setw(30) << "Excitation amplitude"
diff --git a/source/module_lr/lr_spectrum_velocity.cpp b/source/module_lr/lr_spectrum_velocity.cpp
index 5a6927f74d..1c3778f973 100644
--- a/source/module_lr/lr_spectrum_velocity.cpp
+++ b/source/module_lr/lr_spectrum_velocity.cpp
@@ -23,9 +23,14 @@ namespace LR
         return vR;
     }
 
-    inline double lorentz_delta(const double freq_diff, const double eta)
+    inline double lorentz_delta(const double dfreq_au, const double eta_au)
     {
-        return eta / (freq_diff * freq_diff + eta * eta) / M_PI;
+        return eta_au / (dfreq_au * dfreq_au + eta_au * eta_au) / M_PI;
+    }
+    inline double gauss_delta(const double dfreq_au, const double eta_au)
+    {
+        const double c = eta_au / std::sqrt(2. * std::log(2.));
+        return std::exp(-dfreq_au * dfreq_au / (2 * c * c)) / (std::sqrt(2 * M_PI) * c);
     }
 
     template<typename T> inline ModuleBase::Vector3<T> convert_vector_to_vector3(const std::vector<std::complex<double>>& vec);
@@ -109,13 +114,13 @@ namespace LR
         // 4*pi^2/V * mean_squared_dipole *delta(w-Omega_S)
         std::ofstream ofs(PARAM.globalv.global_out_dir + "absorption.dat");
         if (GlobalV::MY_RANK == 0) { ofs << "Frequency (eV) | wave length(nm) | Absorption (a.u.)" << std::endl; }
-        const double fac = 4 * M_PI * M_PI / ucell.omega * ModuleBase::e2 / this->nk;  // e2: Ry to Hartree in the denominator
+        const double fac = 4 * M_PI * M_PI / ucell.omega / this->nk;
         for (int f = 0;f < freq.size();++f)
         {
             double abs_value = 0.0;
             for (int i = 0;i < nstate;++i)
             {
-                abs_value += this->mean_squared_transition_dipole_[i] * lorentz_delta((freq[f] - eig[i]), eta);
+                abs_value += this->mean_squared_transition_dipole_[i] * lorentz_delta((freq[f] - eig[i]) / ModuleBase::e2, eta / ModuleBase::e2); // e2: Ry to Hartree 
             }
             abs_value *= fac;
             if (GlobalV::MY_RANK == 0) { ofs << freq[f] * ModuleBase::Ry_to_eV << "\t" << 91.126664 / freq[f] << "\t" << abs_value << std::endl; }
diff --git a/source/module_ri/module_exx_symmetry/irreducible_sector.h b/source/module_ri/module_exx_symmetry/irreducible_sector.h
index 6f9ab3cb8c..7a12b9a232 100644
--- a/source/module_ri/module_exx_symmetry/irreducible_sector.h
+++ b/source/module_ri/module_exx_symmetry/irreducible_sector.h
@@ -127,8 +127,6 @@ namespace ModuleSymmetry
 
         /// symmetry info for BvK supercell
         std::vector<int> isymbvk_to_isym_;
-        std::vector<ModuleBase::Matrix3> bvk_gmatrix_;
-        std::vector<ModuleBase::Vector3<double>> bvk_gtrans_;
         int bvk_nsym_;
 
         friend class Symmetry_rotation;
diff --git a/source/module_ri/module_exx_symmetry/irreducible_sector_bvk.cpp b/source/module_ri/module_exx_symmetry/irreducible_sector_bvk.cpp
index 18c9954bce..5efd046d73 100644
--- a/source/module_ri/module_exx_symmetry/irreducible_sector_bvk.cpp
+++ b/source/module_ri/module_exx_symmetry/irreducible_sector_bvk.cpp
@@ -40,13 +40,15 @@ namespace ModuleSymmetry
         ModuleBase::TITLE("Irreducible_Sector", "gen_symmetry_BvK");
         auto set_matrix3 = [](const ModuleBase::Vector3<double>& a1, const ModuleBase::Vector3<double>& a2, const ModuleBase::Vector3<double>& a3)
             -> ModuleBase::Matrix3 {return ModuleBase::Matrix3(a1.x, a1.y, a1.z, a2.x, a2.y, a2.z, a3.x, a3.y, a3.z);};
-
+        auto set_bvk_same_as_ucell = [&symm, this]()->void
+            {
+                this->bvk_nsym_ = symm.nrotk;
+                this->isymbvk_to_isym_.resize(symm.nrotk);
+                for (int isym = 0;isym < symm.nrotk;++isym) { this->isymbvk_to_isym_[isym] = isym; }
+            };
         if (bvk_period[0] == bvk_period[1] && bvk_period[0] == bvk_period[2])
         {   //the BvK supercell has the same symmetry as the original cell
-            this->bvk_nsym_ = symm.nrotk;
-            this->isymbvk_to_isym_.resize(symm.nrotk);
-            for (int isym = 0;isym < symm.nrotk;++isym) { this->isymbvk_to_isym_[isym] = isym;
-}
+            set_bvk_same_as_ucell();
             return;
         }
 
@@ -129,23 +131,30 @@ namespace ModuleSymmetry
         bvk_op.resize(bvk_nop);
         int bvk_npg, bvk_nsg, bvk_pgnum, bvk_sgnum;
         std::string bvk_pgname, bvk_sgname;
-        this->bvk_gmatrix_.resize(48);
-        this->bvk_gtrans_.resize(48);
+        std::vector<ModuleBase::Matrix3> bvk_gmatrix(48);
+        std::vector<ModuleBase::Vector3<double>> bvk_gtrans(48);
         symm.getgroup(bvk_npg, bvk_nsg, GlobalV::ofs_running, bvk_nop,
-            bvk_op.data(), this->bvk_gmatrix_.data(), this->bvk_gtrans_.data(),
+            bvk_op.data(), bvk_gmatrix.data(), bvk_gtrans.data(),
             bvk_dpos.data(), bvk_rot_dpos.data(), order_index.data(),
             bvk_itmin_type, bvk_itmin_start, bvk_istart.data(), bvk_na.data());
-        this->bvk_gmatrix_.resize(bvk_nsg);
-        this->bvk_gtrans_.resize(bvk_nsg);
+        bvk_gmatrix.resize(bvk_nsg);
+        bvk_gtrans.resize(bvk_nsg);
         this->bvk_nsym_ = bvk_nsg;
-        symm.pointgroup(bvk_npg, bvk_pgnum, bvk_pgname, this->bvk_gmatrix_.data(), GlobalV::ofs_running);
+        // bvk suppercell cannot have higher symmetry than the original cell
+        if (this->bvk_nsym_ > symm.nrotk)
+        {
+            std::cout << "reset bvk symmetry to the same as the original cell" << std::endl;
+            set_bvk_same_as_ucell();
+            return;
+        }
+        symm.pointgroup(bvk_npg, bvk_pgnum, bvk_pgname, bvk_gmatrix.data(), GlobalV::ofs_running);
         ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "POINT GROUP OF BvK SCELL", bvk_pgname);
-        symm.pointgroup(bvk_nsg, bvk_sgnum, bvk_sgname, this->bvk_gmatrix_.data(), GlobalV::ofs_running);
+        symm.pointgroup(bvk_nsg, bvk_sgnum, bvk_sgname, bvk_gmatrix.data(), GlobalV::ofs_running);
         ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "POINT GROUP IN SPACE GROUP OF BvK SCELL", bvk_sgname);
-        symm.gmatrix_convert_int(this->bvk_gmatrix_.data(), this->bvk_gmatrix_.data(), bvk_nsg, bvk_min_optlat, lat.latvec);
-        symm.gtrans_convert(this->bvk_gtrans_.data(), this->bvk_gtrans_.data(), bvk_nsg, bvk_min_optlat, lat.latvec);
+        symm.gmatrix_convert_int(bvk_gmatrix.data(), bvk_gmatrix.data(), bvk_nsg, bvk_min_optlat, lat.latvec);
+        symm.gtrans_convert(bvk_gtrans.data(), bvk_gtrans.data(), bvk_nsg, bvk_min_optlat, lat.latvec);
         // get map from bvk-op to original op
-        this->isymbvk_to_isym_ = get_isymbvk_to_isym_map(this->bvk_gmatrix_, symm);
+        this->isymbvk_to_isym_ = get_isymbvk_to_isym_map(bvk_gmatrix, symm);
         return;
     }
 
diff --git a/source/module_ri/module_exx_symmetry/symmetry_rotation.cpp b/source/module_ri/module_exx_symmetry/symmetry_rotation.cpp
index c73d587ff2..27cfc8e723 100644
--- a/source/module_ri/module_exx_symmetry/symmetry_rotation.cpp
+++ b/source/module_ri/module_exx_symmetry/symmetry_rotation.cpp
@@ -210,7 +210,7 @@ namespace ModuleSymmetry
     // gmatc should be a rotation matrix, i.e. det(gmatc)=1
     TCdouble Symmetry_rotation::get_euler_angle(const ModuleBase::Matrix3& gmatc) const
     {
-        double threshold = 1e-8;
+        double threshold = this->eps_;
         double alpha, beta, gamma;
         if (std::fabs(gmatc.e32) > threshold || std::fabs(gmatc.e31) > threshold) // sin(beta) is not zero
         {
diff --git a/tests/integrate/223_NO_zero_atom/INPUT b/tests/integrate/223_NO_zero_atom/INPUT
new file mode 100644
index 0000000000..abe15ed4a0
--- /dev/null
+++ b/tests/integrate/223_NO_zero_atom/INPUT
@@ -0,0 +1,27 @@
+INPUT_PARAMETERS
+#Parameters (1.General)
+suffix			autotest
+calculation     	scf
+
+nbands			6
+symmetry		0
+pseudo_dir      	../../PP_ORB/
+orbital_dir     	../../PP_ORB/
+
+#Parameters (2.Iteration)
+ecutwfc			20
+
+#Parameters (3.Basis)
+basis_type		lcao
+
+#Parameters (4.Smearing)
+smearing_method		gauss
+smearing_sigma		0.0002
+
+#Parameters (5.Mixing)
+mixing_type		broyden
+mixing_beta		0.7
+cal_force		1
+test_force		1
+cal_stress		1
+test_stress		1
diff --git a/tests/integrate/223_NO_zero_atom/KPT b/tests/integrate/223_NO_zero_atom/KPT
new file mode 100644
index 0000000000..f5f7f4ec34
--- /dev/null
+++ b/tests/integrate/223_NO_zero_atom/KPT
@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+2 2 2 0 0 0
diff --git a/tests/integrate/223_NO_zero_atom/STRU b/tests/integrate/223_NO_zero_atom/STRU
new file mode 100644
index 0000000000..a8d7fba106
--- /dev/null
+++ b/tests/integrate/223_NO_zero_atom/STRU
@@ -0,0 +1,27 @@
+ATOMIC_SPECIES
+Al 13 Al_ONCV_PBE-1.0.upf upf201
+Si 14 Si_ONCV_PBE-1.0.upf upf201
+
+NUMERICAL_ORBITAL
+Al_gga_8au_100Ry_4s4p1d.orb
+Si_gga_8au_100Ry_2s2p1d.orb
+
+LATTICE_CONSTANT
+10.2  // add lattice constant
+
+LATTICE_VECTORS
+0.5 0.5 0.0
+0.5 0.0 0.5
+0.0 0.5 0.5
+
+ATOMIC_POSITIONS
+Direct 
+
+Al
+0.0
+0
+Si      // Element type
+0.0     // magnetism
+2
+0.00 0.00 0.00 1 1 1
+0.25 0.25 0.25 1 1 1
diff --git a/tests/integrate/223_NO_zero_atom/jd b/tests/integrate/223_NO_zero_atom/jd
new file mode 100644
index 0000000000..e8fb4e42b6
--- /dev/null
+++ b/tests/integrate/223_NO_zero_atom/jd
@@ -0,0 +1 @@
+test calculation when atom number is zero with LCAO.
diff --git a/tests/integrate/223_NO_zero_atom/result.ref b/tests/integrate/223_NO_zero_atom/result.ref
new file mode 100644
index 0000000000..af95ced2af
--- /dev/null
+++ b/tests/integrate/223_NO_zero_atom/result.ref
@@ -0,0 +1,5 @@
+etotref -211.4490868003635171
+etotperatomref -105.7245434002
+totalforceref 0.000000
+totalstressref 366.791415
+totaltimeref 2.69
diff --git a/tests/integrate/CASES_CPU.txt b/tests/integrate/CASES_CPU.txt
index 6e6093523a..134be1fbae 100644
--- a/tests/integrate/CASES_CPU.txt
+++ b/tests/integrate/CASES_CPU.txt
@@ -188,6 +188,7 @@
 #220_NO_KP_MD_MSST_level2
 220_NO_KP_MD_NVT
 220_NO_KP_MD_wfc_out
+223_NO_zero_atom
 #230_NO_KP_MD_TD
 240_NO_KP_15_SO
 240_NO_KP_15_SO_average